普段使いはPythonだが、用途によっては計算時間が開発のボトルネックになることがあって、今後も増えそうな予感。Pythonスクリプトの高速化の方法は何通りもあるが、どれもイマイチだなぁということで計算機のパワーでゴリ押しすることが多い。 動的型付け言語の様な書き方が出来て、数学ライブラリが使えて、適度に最適化されたコンパイル出来る言語が欲しいという需要が自分の中にあった。
JuliaというPythonやMatlabのようなスクリプト言語っぽい書き方の出来るコンパイル言語がイケてそうだったので、自分用に使う機能などを試すためにロケットの飛翔計算について初歩的な内容を書いている。
が試せればOK。
結果としては補間がダメだった。Interpolationというパッケージが所望の挙動をせず、PyCallによってScipyのinterp1d関数を呼び出してもエラーが消えなかった(自分の書き方が悪かっただけかも)。それ以外は使えて、早いので、使えるかも。。。というところ。
以下は試してみた結果。
using PyPlot
rc("font",family ="Osaka",size=16)
using ODE
ロケットの垂直方向に打ち上がって落ちるまでを考える。ロケットの運動は1自由度(1次元で垂直方向のみ)、3自由度(垂直と横方向の2次元平面上と機体の傾きの合わせて3自由度)、6自由度(空間3軸、姿勢3軸の6自由度)の運動を計算することがあるが、ここでは簡単のために1自由度で考える。
運動方程式は $$ m\frac{d^2y}{dt^2} = F - D - mg $$ ここで、
記号 | 単位 | 内容 |
---|---|---|
m | kg | 質量 |
y | m | 垂直方向位置 |
F | N | 推力 |
D | N | 抗力 |
g | kg/s^2 | 重力加速度 |
ロケットの比推力をIsp[s]とすると、推進剤を噴射することにより減るロケットの質量は $$\dot{m} = - \frac{F}{I_{SP} \cdot g}$$
1階の連立常微分方程式の形で記述できれば、ルンゲクッタ法などによって数値計算することが出来る。連立常微分方程式の形にすると
$$ \frac{d}{dt}\begin{bmatrix} m \\ \dot{y} \\ y \end{bmatrix} = \begin{bmatrix} - \frac{F}{I_{SP} \cdot g} \\ \frac{F-D}{m} -g \\ \dot{y} \end{bmatrix} $$それぞれの変数について考えていく。
ただし、g0:地表面での重力加速度[kg/s2]、R:地球半径[m]、h:地表からの高度[m]
function gravity(altitude)
# altitude: meter
g0 = 9.80665
radius_earth = 6378137
g = g0 * (radius_earth / (radius_earth + altitude)) ^2
end
gravity (generic function with 1 method)
# gravity関数のテスト
test_alt = 1:1000:10000000
test_g = zeros(length(test_alt))
index = 1
for i = test_alt
test_g[index] = gravity(i)
index += 1
end
plot(test_alt/1000, test_g)
ylim(0,10)
xlabel("高度 km");ylabel("重力加速度 kg/s2");grid()
抗力は、 $$D = Cd \cdot \frac{1}{2}\rho S V^2$$
抗力係数Cdはマッハ数の関数として、下記のようなグラフを利用 音速や空気密度はUS Standard Atmosphere 1976を使用。
function air(altitude)
# 標準大気モデルを用いた、高度による温度、音速、大気圧、空気密度の関数
# 高度は基準ジオポテンシャル高度を元にしている。
# 標準大気の各層ごとの気温減率から定義式を用いて計算している。
# Standard Atmosphere 1976 ISO 2533:1975
# 中間圏高度86kmまでの気温に対応している。それ以上は国際標準大気に当てはまらないので注意。
# cf. http://www.pdas.com/hydro.pdf
# @param altitude 高度[m]
# @return T 温度[K]
# @return a 音速[m/s]
# @return P 気圧[Pa]
# @return rho 空気密度[kg/m3]
h = altitude
g = 9.80655
gamma = 1.4
R = 287.0531
# height of atomospheric layer
HAL = [0, 11000, 20000, 32000, 47000, 51000, 71000, 84852]
# Lapse Rate Kelvin per meter
LR = [-0.0065, 0.0, 0.001, 0.0028, 0, -0.0028, -0.002, 0.0]
# Tempareture Kelvin
T0 = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.95]
# Pressure Pa
P0 = [101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564, 0.3734]
if ( h >= HAL[1] && h < HAL[2])
k = 1
elseif ( h >= HAL[2] && h < HAL[3])
k = 2
elseif ( h >= HAL[3] && h < HAL[4])
k = 3
elseif ( h >= HAL[4] && h < HAL[5])
k = 4
elseif ( h >= HAL[5] && h < HAL[6])
k = 5
elseif ( h >= HAL[6] && h < HAL[7])
k = 6
elseif ( h >= HAL[7] && h < HAL[8])
k = 7
elseif ( h >= HAL[8])
k = 8
else
k = 1
end
T = T0[k] + LR[k] * (h - HAL[k])
a = sqrt(T * gamma * R)
if LR[k] != 0
P = P0[k] * ((T0[k] + LR[k] * (h - HAL[k])) / T0[k]) ^ (g / -LR[k] / R)
else
P = P0[k] * exp(g / R * (HAL[k] - h) / T0[k])
end
rho = P / R / T
[T, a, P, rho]
end
air (generic function with 1 method)
# air関数のテスト
test_alt = 1:100:100000
test_air = zeros(length(test_alt))
index = 1
for i = test_alt
T, a, P, test_air[index] = air(i)
index += 1
end
plot(test_alt/1000, test_air)
xlabel("高度 km");ylabel("空気密度 kg/m3");grid()
# using PyCall
# @pyimport scipy.interpolate as interp
# 補間が使えなかった。。。
function Cd(mach)
# M-Vロケットの抗力係数を参考にした
# @param Mach マッハ数[-]
# @return self.Cd 抗力係数[-]
M = [0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0];
cd_row = [0.28, 0.28, 0.28, 0.29, 0.35, 0.64, 0.67, 0.69, 0.66, 0.62, 0.58, 0.55,
0.48, 0.42, 0.38, 0.355, 0.33]
# 線形補間のアルゴリズム、外挿はせずに
for i = 1:length(M)-1
if mach >= M[i] && mach < M[i+1]
alpha = (mach - M[i]) / (M[i+1] - M[i])
Cd = cd_row[i] + alpha * (cd_row[i+1] - cd_row[i])
break
end
end
if mach < M[1]
Cd = cd_row[1]
elseif mach >= M[end]
Cd = cd_row[end]
end
Cd
end
Cd (generic function with 1 method)
test_mach = 0:0.1:7
test_Cd = zeros(length(test_mach))
index = 1
for i = test_mach
test_Cd[index] = Cd(i)
index += 1
end
plot(test_mach, test_Cd)
ylim(0,1)
xlabel("マッハ数 -");ylabel("Cd -");grid()
# ロケットTypeを作る。ダイナミクスなどの関数にまとめてパラメータを入れるため
type Rocket
vy::Float64
y::Float64
Isp::Float64
thrust::Float64
burn_time::Float64
mass_init::Float64
area::Float64
end
rocket = Rocket(0.0,0.0,200.0,500.0,3.0,10.0,0.5)
Rocket(0.0,0.0,200.0,500.0,3.0,10.0,0.5)
function drag(rocket::Rocket)
T, a, P, rho = air(rocket.y)
mach = rocket.vy / a
drag = 0.5 * Cd(mach) * rho * rocket.area * rocket.vy ^2
end
drag (generic function with 1 method)
test_vel = 0:10:700
test_drag = zeros(length(test_vel))
index = 1
for i = test_vel
rocket = Rocket(i, 0, 200, 500, 3, 10, 0.12)
test_drag[index] = drag(rocket)
index += 1
end
plot(test_vel, test_drag, label="面積0.12m2のとき")
# ylim(0,1)
xlabel("速度 m/s");ylabel("抗力 N");grid()
legend(loc="best")
PyObject <matplotlib.legend.Legend object at 0x313ac3b10>
ここで常微分方程式をもう一度書いておく、JuliaのODEパッケージのodexx関数では連立常微分方程式の関数を記述するとサクッと数値積分してくれる。
$$ \frac{d}{dt}\begin{bmatrix} m \\ \dot{y} \\ y \end{bmatrix} = \begin{bmatrix} - \frac{F}{I_{SP} \cdot g} \\ \frac{F-D}{m} -g \\ \dot{y} \end{bmatrix} $$高度が0以下の時は計算を止めている。
function flight_dynamics(t, x, rocket::Rocket)
# x = [mass, vy, y]
# type Rocket : vy, y, Isp, thrust, burn_time, mass_init, area
dx = similar(x)
g0 = 9.80655
if t < rocket.burn_time
thrust = rocket.thrust
else
thrust = 0
end
m_dot = thrust / rocket.Isp / g0
rocket.vy = x[2]; rocket.y = x[3]
if x[2] > 0
D = drag(rocket)
else
D = - drag(rocket)
end
dx[1] = -m_dot
dx[2] = 1 / x[1] * (thrust - D) - gravity(x[3])
dx[3] = x[2]
if x[3] < 0.0
dx[1] = 0; dx[2] = 0; dx[3] = 0
end
dx
end
flight_dynamics (generic function with 1 method)
せっかくなので既に実用化されているロケットで計算、グラフ化をしてみる。ここではJAXAが鹿児島県内之浦から打上げている観測ロケットであるS-310を見てみる。
公式資料が見つからなかったから、S-510ロケットと同じような比推力の推進剤を使っていると仮定して、地上と真空中の比推力が変わらないというかなり大雑把な近似をしてしまって計算する。 (かなり大雑把な値の取り方と近似をしているので注意)
名称 | 値 |
---|---|
初期質量 | 722 kg |
直径 | 30 cm |
推力 | 31800 N |
比推力 | 260 s |
燃焼時間 | 28 秒 |
t = [0., 400]
Isp = 200; thrust = 31800; burn_time = 30; mass_init = 722; area = 0.07 # S-310 rocket
rocket = Rocket(0.0, 0.0, Isp, thrust, burn_time, mass_init, area)
@time t,y = ode23s((t, x) -> flight_dynamics(t,x,rocket), [rocket.mass_init, 0.0, 0.0], t; abstol=1e-5, reltol=1e-5,
maxstep=1e11/10, minstep=1e11/1e18)
Void
Void
0.028209 seconds (176.61 k allocations: 10.181 MB)
# ==== PLOT ====
mass = zeros(length(t))
velocity = zeros(length(t))
altitude = zeros(length(t))
for i = 1:length(t)
mass[i] = y[i][1]
velocity[i] = y[i][2]
altitude[i] = y[i][3]
end
time_S310 = t
mass_S310 = mass
vel_S310 = velocity
alt_S310 = altitude
subplot(311)
title("JAXA S-310")
plot(t,mass); ylim(0, rocket.mass_init*1.1)
ylabel("質量 kg");grid()
subplot(312)
plot(t, velocity)
ylabel("速度 m/s")
axhline(330, linewidth=1, linestyle = "--")
subplot(313)
plot(t, altitude/1000)
ylabel("高度 km")
xlabel("時間 s");grid()
上の計算では最大高度が150kmを少し超える程度になっている、Wikipediaによると190kmとか210kmなどと出ているので、実際にはもう少し性能が良いようだ。
Blue Originという会社はそもそもかなり秘密主義なので、こちらも公式に性能表が出ていない。Wikipediaなどの怪しい情報を元に考えてみる。 垂直上昇の能力のみを考えてみる。
youtude 307,000 Feet |Blue Origin
名称 | 値 |
---|---|
初期質量 | 29000 kg |
直径 | 2.5 m |
推力 | 420 kN |
比推力 | 360 s |
燃焼時間 | 120 秒 |
初期質量や直径については根拠の無い数字なので注意です。
t = [0., 400]
Isp = 360; thrust = 420000; burn_time = 120; mass_init = 29000; area = 4.9 # Blue Origin New Shepard
rocket = Rocket(0.0, 0.0, Isp, thrust, burn_time, mass_init, area)
@time t,y = ode23s((t, x) -> flight_dynamics(t,x,rocket), [rocket.mass_init, 0.0, 0.0], t; abstol=1e-4, reltol=1e-4,
maxstep=1e11/10, minstep=1e11/1e18)
Void
Void
0.030199 seconds (79.22 k allocations: 4.507 MB, 32.61% gc time)
# ==== PLOT ====
mass = zeros(length(t))
velocity = zeros(length(t))
altitude = zeros(length(t))
for i = 1:length(t)
mass[i] = y[i][1]
velocity[i] = y[i][2]
altitude[i] = y[i][3]
end
time_NS = t
mass_NS = mass
vel_NS = velocity
alt_NS = altitude
subplot(311)
title("Blue Origin New Shepard")
plot(t,mass); ylim(0, rocket.mass_init*1.1)
ylabel("質量 kg");grid()
subplot(312)
plot(t, velocity)
ylabel("速度 m/s");grid()
axhline(330, linewidth=1, linestyle = "--")
subplot(313)
plot(t, altitude/1000)
ylabel("高度 km")
xlabel("時間 s");grid()
2種類の高度100kmを超える垂直打上げロケットの飛翔プロファイルを計算したので、比較してみる。
subplot(311)
title("S-310 VS New Shepard")
plot(time_S310,mass_S310, label="S-310")
plot(time_NS,mass_NS, label="New Shepard")
ylim(0,30000);xlim(0,300)
legend(loc = "best", fontsize = 10)
# plot(time_NS,mass_NS)
ylabel("質量 kg")
subplot(312)
plot(time_S310, vel_S310, label="S-310")
plot(time_NS, vel_NS, label="New Shepard")
ylabel("速度 m/s")
axhline(330, linewidth=1, linestyle = "--")
ylim(-500, 2000);grid();xlim(0,300)
subplot(313)
plot(time_S310, alt_S310/1000, label="S-310")
plot(time_NS, alt_NS/1000, label="New Shepard")
ylabel("高度 km")
xlabel("時間 s")
grid();ylim(0,200);xlim(0,300)
(0,300)
S−310のような固体燃料ロケットは構造が簡単で大推力が出せる特性がある。一方、New Shepardに使われているBE-3エンジンは液体水素/液体酸素のロケットであり、高性能(高比推力)ではあるが、推力はそこそこという特性がある。
エンジンの違いによって、高度100km以上に到達するロケットにこれだけ違いが出てくる。 特にNew Shepardは上部に有人カプセルを載せており、人にかかるGを低減するためにあえて液体ロケットを使っている。