大気圏突入時の熱流束を知りたいー3. 落下してくる宇宙船の運動モデルー
blunt body(Appolloの宇宙船カプセルみたいな形の物体)が落下してくる時の、速度を計算する方法を調べた。
Apollo Reentry GuidanceのAppendix Aにある微分方程式を少し修正して、以下の式を得る。
元の式の船体角度γの微分方程式には、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にある図と整合している。
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