大気圏突入時の熱流束を知りたいー3. 落下してくる宇宙船の運動モデルー

blunt body(Appolloの宇宙船カプセルみたいな形の物体)が落下してくる時の、速度を計算する方法を調べた。

Apollo Reentry GuidanceのAppendix Aにある微分方程式を少し修正して、以下の式を得る。

f:id:imakov:20211226231147p:plain

元の式の船体角度γの微分方程式には、lifting parameterを含んだ項があるのだが、上の式では無視している。根拠として薄弱なのは承知だが、International Journal of Aerospace Engineeringの"High-Fidelity Aerothermal Engineering Analysis for Planetary Probes Using DOTNET Framework and OLAP Cubes Database"という論文で、lifting parameterを省略しているのを参考にさせてもらった。現実問題としては、Apolloの簡単な資料にlifting parameterが載っていないので、同パラメータの項を無視できるのは助かった。

上述の式で、同Apolloの資料にあるパラメータ(Apollo 6、突入速度10 km/s、突入角度-5.9°、ballistic係数395.8 kg/m2)を用いて計算した結果が次の図である。大気密度は前回のblog記事の式を用いた。縦軸が高度、横軸が速度で、各高度での速度が分かるようになっている。ちなみに、この図は上述のInternational Journal of Aerospace Engineeringにある図と整合している。 f:id:imakov:20211226232246p:plain

close all
clear variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Initial values
Ve_ms = 10e3; % Velocity at entrance to earth's atmosphere (km/s)
h0_m  = 76.2e3; % Initial altitude (250 000 feet)
gamma0 = -5.9/180*pi; % entry angle

% Solve differential equation
t_range = [1, 3e2];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,y] = ode45(@R415, t_range, [Ve_ms; h0_m; gamma0;], opts);
wh = y(:,2)>=0;
t_s = t(wh)'; % time (s)
V_ms = y(wh,1)'; % velosity (m/s)
h_m = y(wh,2)'; % altitude (m)
gamma_deg = y(wh,3)'/pi*180; % angle (degree)

%Plot results
figure();
subplot(3,1,1)
plot(t_s, h_m*1e-3)
xlabel('time (s)')
ylabel('Altitude (km)')
grid on

subplot(3,1,2)
plot(t_s, V_ms*1e-3)
xlabel('time (s)')
ylabel('Velocity (km/s)')
grid on

subplot(3,1,3)
plot(t_s, gamma_deg)
xlabel('time (s)')
ylabel('angle (deg)')
grid on

%Plot results
figure();
plot(V_ms*1e-3, h_m*1e-3)
xlabel('velocity (km/s)')
ylabel('Altitude (km)')
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = R415(t, y)
% y(1): velocity (m/s)
% y(2): altitude (m)
% y(3): angle (rad)

Cb_kgm2 = 395.8 ; % ballistic coefficient (kg/m2)
g= 9.8; % gravity acceleration(m/s2)
Re = 6.3781e6; % Earth radius (m)

dydt = [-0.5*exp_rho(y(2))*y(1)^2/Cb_kgm2 - g*sin(y(3)); 
    y(1)*sin(y(3));
    (y(1)/(Re+y(2)) - g/y(1))*cos(y(3))];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Allen-Eggersのexponential atmosphere model
function rho_kgm3 = exp_rho(h_m)
% h: altitude (m)
rho0_kgm3 = 1.752 ; % 0.0034 slugs/ft3
h0_m = 6705.6; % 22000 (1/ft)
rho_kgm3 = rho0_kgm3*exp(-h_m/h0_m);
end