ターミナル速度方程式

ターミナル速度方程式

ロケット工学に関連する方程式の一例として、ロケットの運動方程式であるターミナル速度方程式を挙げて、それをPythonで3Dグラフ化してみましょう。

ターミナル速度方程式：

$$
v_t = \sqrt{\frac{2 \cdot g \cdot I_{sp} \cdot \ln\left(\frac{m_0}{m_f}\right)}{A \cdot \rho_c \cdot C_d}}
$$

ここで、各変数の意味は以下の通りです：

  • (v_t)：ターミナル速度（Terminal Velocity）
  • (g)：重力加速度（通常は地球の場合は約9.8 m/s(^2)）
  • (I_{sp})：比推力（Specific Impulse）
  • (m_0)：初期質量
  • (m_f)：最終質量
  • (A)：断面積
  • (\rho_c)：大気密度
  • (C_d)：抗力係数

以下は、この方程式を3Dグラフ化するためのPythonコードの一例です。

ただし、実際の値やパラメータは具体的なロケットに依存します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# ターミナル速度方程式
def terminal_velocity(m0, mf, Isp, A, rho_c, Cd):
    g = 9.8  # 重力加速度 (m/s^2)
    return np.sqrt((2 * g * Isp * np.log(m0 / mf)) / (A * rho_c * Cd))

# パラメータの設定（例）
m0 = 1000  # 初期質量 (kg)
mf = 500   # 最終質量 (kg)
Isp = 300  # 比推力 (s)
A = 5      # 断面積 (m^2)
rho_c = 1  # 大気密度 (kg/m^3)
Cd = 0.5   # 抗力係数

# 初期質量と最終質量の範囲を設定
m0_values = np.linspace(800, 1200, 50)
mf_values = np.linspace(400, 600, 50)

# 3Dグリッドの生成
m0_grid, mf_grid = np.meshgrid(m0_values, mf_values)

# ターミナル速度の計算
vt_values = terminal_velocity(m0_grid, mf_grid, Isp, A, rho_c, Cd)

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(m0_grid, mf_grid, vt_values, cmap='viridis')

# 軸ラベルとタイトルの設定
ax.set_xlabel('Initial Mass (kg)')
ax.set_ylabel('Final Mass (kg)')
ax.set_zlabel('Terminal Velocity (m/s)')
ax.set_title('Terminal Velocity vs. Initial and Final Mass')

plt.show()

この例では、初期質量と最終質量を変化させたときのターミナル速度を3Dグラフで表示しています。

[実行結果]

パラメータや方程式を実際の状況に合わせて変更することが重要です。

ソースコード解説

以下に、コードの主な要素を章立てして説明します。

1. 関数 terminal_velocity の定義:

1
2
3
def terminal_velocity(m0, mf, Isp, A, rho_c, Cd):
   g = 9.8  # 重力加速度 (m/s^2)
   return np.sqrt((2 * g * Isp * np.log(m0 / mf)) / (A * rho_c * Cd))
  • ロケットのターミナル速度を計算する関数です。
    関数は初期質量 m0、最終質量 mf、比推力 Isp、断面積 A、大気密度 rho_c、抗力係数 Cdを引数に取ります。

2. パラメータの設定:

1
2
3
4
5
6
m0 = 1000  # 初期質量 (kg)
mf = 500   # 最終質量 (kg)
Isp = 300  # 比推力 (s)
A = 5      # 断面積 (m^2)
rho_c = 1  # 大気密度 (kg/m^3)
Cd = 0.5   # 抗力係数
  • ロケットと環境に関するパラメータを設定します。

3. 初期質量と最終質量の範囲の設定:

1
2
m0_values = np.linspace(800, 1200, 50)
mf_values = np.linspace(400, 600, 50)
  • グラフに表示する初期質量と最終質量の範囲を設定します。

4. 3Dグリッドの生成:

1
m0_grid, mf_grid = np.meshgrid(m0_values, mf_values)
  • 3Dグラフの横軸と縦軸に対応する初期質量と最終質量の組み合わせを作成します。

5. ターミナル速度の計算:

1
vt_values = terminal_velocity(m0_grid, mf_grid, Isp, A, rho_c, Cd)
  • 3Dグリッド上の各点でのターミナル速度を計算します。

6. 3Dプロット:

1
2
3
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(m0_grid, mf_grid, vt_values, cmap='viridis')
  • matplotlib を使用して、初期質量、最終質量、および対応するターミナル速度の3Dプロットを作成します。

7. グラフの装飾:

1
2
3
4
ax.set_xlabel('Initial Mass (kg)')
ax.set_ylabel('Final Mass (kg)')
ax.set_zlabel('Terminal Velocity (m/s)')
ax.set_title('Terminal Velocity vs. Initial and Final Mass')
  • グラフに軸ラベルやタイトルを追加します。

8. グラフの表示:

1
plt.show()
  • 最終的に作成された3Dグラフを表示します。

結果解説

[実行結果]

この3Dグラフは、ロケットの運動方程式から計算されるターミナル速度を初期質量と最終質量の両方に依存して表示しています。

以下にグラフの主な要素を説明します。

X軸（Initial Mass）:

グラフの底面にあたり、ロケットの初期質量を表します。
この軸を変化させると、グラフの左右方向が変わります。

Y軸（Final Mass）:

グラフの底面にあたり、ロケットの最終質量を表します。
この軸を変化させると、グラフの前後方向が変わります。

Z軸（Terminal Velocity）:

グラフの縦方向に伸びる軸で、計算されたターミナル速度を表します。
この軸を変化させると、グラフの上下方向が変わります。

グラフの表面は、各点でのターミナル速度を色で表しており、色が濃いほど速度が大きいことを示しています。

このような3Dグラフを通じて、初期質量と最終質量がターミナル速度に与える影響を視覚的に理解できます。