大気圏突入時の熱流束を知りたいー4. Sutton-Gravesの熱伝達モデルー

前のblogに書いた大気密度のモデルと宇宙船の運動モデルに加えて、Sutton-Gravesの熱伝達モデルを使うと大気から宇宙船への熱流束を評価することができる。

J.FayとF. Riddellが考案（Journal of Aeronautical Sciences, vol. 25, no. 2, 1958）し、K. SuttonとR. A. Gravesが簡易化（NASA TR-376, 1971）した便利な式があって、大気から宇宙船への熱流束は次の式で表せる。

これまでblogで書いてきたことを総動員すると、Apollo 6に対する各高度の熱流束は次の図のように計算できる。高度が高いところでは、船体の速度が速いが空気が薄いので熱流束は低い。また、高度が低いところでは、空気が濃いが船体の速度が遅いので熱流束は低い。ということで、どこか途中のところでピークがあるような分布となっている。私の計算では最大熱流束は260 (W/cm²)となっている。前回のblogで参照したApolloの簡単な資料によると、最大熱流束は240 (W/cm²)とあるので、簡易的な計算でも十分使えるのが分かる。

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)

%  Allen-Eggersのexponential atmosphere model
rho_kgm3 = exp_rho(h_m);

%  Sutton-Gravesのconvective heating model
k = 1.7415e-4; % Earth
Rn_m = 3; % Nose radius (m)
qs = k*sqrt(rho_kgm3/Rn_m) .* (V_ms).^3 / 1e4; % W/cm2

%PLOT
subplot(3,1,1)
plot(rho_kgm3, h_m*1e-3)
xlabel('density (kg/m3)')
ylabel('Altitude (km)')
grid on

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

subplot(3,1,3)
plot(qs,  h_m*1e-3)
xlabel('Convective heating (W/cm2)')
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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