この文章の位置づけ

普段使いはPythonだが、用途によっては計算時間が開発のボトルネックになることがあって、今後も増えそうな予感。Pythonスクリプトの高速化の方法は何通りもあるが、どれもイマイチだなぁということで計算機のパワーでゴリ押しすることが多い。 動的型付け言語の様な書き方が出来て、数学ライブラリが使えて、適度に最適化されたコンパイル出来る言語が欲しいという需要が自分の中にあった。

JuliaというPythonやMatlabのようなスクリプト言語っぽい書き方の出来るコンパイル言語がイケてそうだったので、自分用に使う機能などを試すためにロケットの飛翔計算について初歩的な内容を書いている。

  • 微分方程式(ODE)
  • 補間
  • グラフ化

が試せればOK。

結果としては補間がダメだった。Interpolationというパッケージが所望の挙動をせず、PyCallによってScipyのinterp1d関数を呼び出してもエラーが消えなかった(自分の書き方が悪かっただけかも)。それ以外は使えて、早いので、使えるかも。。。というところ。

以下は試してみた結果。

In [1]:
using PyPlot
rc("font",family ="Osaka",size=16)
using ODE

ロケット1自由度(1次元)のダイナミクス

ロケットの垂直方向に打ち上がって落ちるまでを考える。ロケットの運動は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} $$ それぞれの変数について考えていく。

重力加速度

$$g(h) = g_0 \left( \frac{R}{R+h} \right) ^2$$ ただし、g0:地表面での重力加速度[kg/s2]、R:地球半径[m]、h:地表からの高度[m]

In [2]:
function gravity(altitude)
    #     altitude: meter
    g0 = 9.80665
    radius_earth = 6378137
    g = g0 * (radius_earth / (radius_earth + altitude)) ^2
end
Out[2]:
gravity (generic function with 1 method)
In [3]:
# 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を使用。

In [4]:
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
Out[4]:
air (generic function with 1 method)
In [5]:
# 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()
In [6]:
# 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
Out[6]:
Cd (generic function with 1 method)
In [7]:
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()
In [8]:
# ロケット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)
Out[8]:
Rocket(0.0,0.0,200.0,500.0,3.0,10.0,0.5)
In [9]:
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
Out[9]:
drag (generic function with 1 method)
In [10]:
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")
Out[10]:
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以下の時は計算を止めている。

In [11]:
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
Out[11]:
flight_dynamics (generic function with 1 method)

S-130ロケットの飛翔

せっかくなので既に実用化されているロケットで計算、グラフ化をしてみる。ここではJAXAが鹿児島県内之浦から打上げている観測ロケットであるS-310を見てみる。

公式資料が見つからなかったから、S-510ロケットと同じような比推力の推進剤を使っていると仮定して、地上と真空中の比推力が変わらないというかなり大雑把な近似をしてしまって計算する。 (かなり大雑把な値の取り方と近似をしているので注意)

名称
初期質量 722 kg
直径 30 cm
推力 31800 N
比推力 260 s
燃焼時間 28 秒
In [21]:
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
  
Out[21]:
Void
0.028209 seconds (176.61 k allocations: 10.181 MB)
In [22]:
# ==== 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()

S-310ロケットの結果の考察

上の計算では最大高度が150kmを少し超える程度になっている、Wikipediaによると190kmとか210kmなどと出ているので、実際にはもう少し性能が良いようだ。

Blue Originの再使用ロケットNew Shepardロケットの場合

Blue Originという会社はそもそもかなり秘密主義なので、こちらも公式に性能表が出ていない。Wikipediaなどの怪しい情報を元に考えてみる。 垂直上昇の能力のみを考えてみる。

youtude 307,000 Feet |Blue Origin

名称
初期質量 29000 kg
直径 2.5 m
推力 420 kN
比推力 360 s
燃焼時間 120 秒

初期質量や直径については根拠の無い数字なので注意です。

In [24]:
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
  
Out[24]:
Void
0.030199 seconds (79.22 k allocations: 4.507 MB, 32.61% gc time)
In [25]:
# ==== 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()

S-310とNew Shepardの比較

2種類の高度100kmを超える垂直打上げロケットの飛翔プロファイルを計算したので、比較してみる。

In [19]:
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)
Out[19]:
(0,300)

S-310とNew Shepardの違い

S−310のような固体燃料ロケットは構造が簡単で大推力が出せる特性がある。一方、New Shepardに使われているBE-3エンジンは液体水素/液体酸素のロケットであり、高性能(高比推力)ではあるが、推力はそこそこという特性がある。

エンジンの違いによって、高度100km以上に到達するロケットにこれだけ違いが出てくる。 特にNew Shepardは上部に有人カプセルを載せており、人にかかるGを低減するためにあえて液体ロケットを使っている。