物理計算の単位

ライブスクリプトを開く

この例では、物理計算において単位を扱う方法を説明します。落下するパラシュート兵の終端速度を、SI 単位と帝国単位の両方で計算します。重力と抗力を考慮して、パラシュート兵の動きを解きます。

はじめに

航空機から飛び出すパラシュート兵を想像してください。パラシュート兵に働く力は、重力とパラシュートからの逆向きの抗力の 2 つしかないと仮定します。抗力は、パラシュート兵の速度の 2 乗に比例します。

パラシュート兵に働く合力は次のように表現できます。

$mass \cdot acceleration = drag force - gravitational force$ ,

$m \frac{\partial}{\partial t} v (t) = c_{d} {v (t)}^{2} - m g$ ,

ここで、

$m$ はパラシュート兵の質量
g は重力加速度
$v (t)$ はパラシュート兵の速度
$c_{d}$ は空気抵抗係数

運動方程式の定義と求解

運動方程式を表現する微分方程式を定義します。

syms g m c_d
syms v(t)
eq = m*diff(v(t),t) + m*g == c_d*v(t)^2

eq = 
 $m \frac{\partial}{\partial t} v (t) + g m = c_{d} {v (t)}^{2}$

パラシュートは $t = 0$ で即座に開くため、方程式 eq は $t \geq 0$ のすべての値に対して有効です。dsolve を初期条件 $v (0) = 0$ で使用して微分方程式を解析的に解きます。解では、パラシュート兵の速度を時間の関数として表現しています。

velocity = simplify(dsolve(eq, v(0) == 0))

velocity = 
 $- \frac{\sqrt{g} \sqrt{m} \tanh (\frac{\sqrt{c_{d}} \sqrt{g} t}{\sqrt{m}})}{\sqrt{c_{d}}}$

空気抵抗係数の単位系を求める

SI 単位系における空気抵抗係数 $c_{d}$ の単位を求めます。

力の SI 単位はニュートン $(N)$ です。ニュートンを基本単位で表すと、 $(\frac{kg \cdot m}{s^{2}})$ です。これらは等しいため、単位換算率は 1 です。

u = symunit;
unitConversionFactor(u.N, u.kg*u.m/u.s^2)

ans =  $1$

抗力 $c_{d} {v (t)}^{2}$ はニュートン $(N)$ で重力 $m g$ と同じ単位でなければなりません。次元解析を使用して $c_{d}$ の単位について解きます。

syms drag_units_SI
drag_units_SI = simplify(solve(drag_units_SI * (u.m / u.s)^2 == u.N))

drag_units_SI = 
 $\frac{kg "kilogram - a physical unit of mass."}{m "meter - a physical unit of length."}$

終端速度の推定

以下の値を定義してパラシュート兵の動きを記述します。

パラシュート兵の質量 $m = 70 kg$
重力加速度 $g = 9.81 m / s^{2}$
空気抵抗係数 $c_{d} = 40 kg / m$

これらの値を速度方程式に代入して、結果を単純化します。

vel_SI = subs(velocity,[g,m,c_d],[9.81*u.m/u.s^2, 70*u.kg, 40*drag_units_SI])

vel_SI = 
 $- \frac{\tanh (\frac{t \sqrt{40 \frac{kg "kilogram - a physical unit of mass."}{m "meter - a physical unit of length."}} \sqrt{\frac{981}{100} \frac{m "meter - a physical unit of length."}{{s "second - a physical unit of time."}^{2}}}}{\sqrt{70 kg "kilogram - a physical unit of mass."}}) \sqrt{70 kg "kilogram - a physical unit of mass."} \sqrt{\frac{981}{100} \frac{m "meter - a physical unit of length."}{{s "second - a physical unit of time."}^{2}}}}{\sqrt{40 \frac{kg "kilogram - a physical unit of mass."}{m "meter - a physical unit of length."}}}$

vel_SI = simplify(vel_SI)

vel_SI = 
 $- \frac{3 \sqrt{763} \tanh (\frac{3 \sqrt{763} t}{35} \frac{1}{s "second - a physical unit of time."})}{20} \frac{m "meter - a physical unit of length."}{s "second - a physical unit of time."}$

有効桁数 3 桁で速度の数値近似を計算します。

digits(3)
vel_SI = vpa(vel_SI)

vel_SI = 
 $- 4.14 \tanh (2.37 t \frac{1}{s "second - a physical unit of time."}) \frac{m "meter - a physical unit of length."}{s "second - a physical unit of time."}$

重力と抗力のバランスが取れるとパラシュート兵の速度は一定に近づきます。これは終端速度と呼ばれ、パラシュートからの抵抗力が重力を相殺する時点で生じます (それ以降は加速なし)。 $t ⟶ \infty$ の極限を得て終端速度を求めます。

vel_term_SI = limit(vel_SI, t, Inf)

vel_term_SI = 
 $- 4.14 \frac{m "meter - a physical unit of length."}{s "second - a physical unit of time."}$

速度の帝国単位への変換

最後に、速度関数を SI 単位から帝国単位に変換します。

vel_Imperial = rewrite(vel_SI,u.ft)

vel_Imperial = 
 $- 13.6 \tanh (2.37 t \frac{1}{s "second - a physical unit of time."}) \frac{ft "foot - a physical unit of length."}{s "second - a physical unit of time."}$

終端速度を変換します。

vel_term_Imperial = rewrite(vel_term_SI,u.ft)

vel_term_Imperial = 
 $- 13.6 \frac{ft "foot - a physical unit of length."}{s "second - a physical unit of time."}$

時間ごとの速度をプロット

速度を時間の関数としてプロットするには、時間 t を秒単位で表現し、t を T で置換します。T は無次元のシンボリック変数です。

syms T
vel_SI = subs(vel_SI, t, T*u.s)

vel_SI = 
 $- 4.14 \tanh (2.37 T) \frac{m "meter - a physical unit of length."}{s "second - a physical unit of time."}$

vel_Imperial = rewrite(vel_SI, u.ft)

vel_Imperial = 
 $- 13.6 \tanh (2.37 T) \frac{ft "foot - a physical unit of length."}{s "second - a physical unit of time."}$

separateUnits を使用して式を単位から切り離します。fplot を使用して式をプロットします。プロットのラベルとして用いるため、symunit2str を使用して単位を string に変換します。

[data_SI, units_SI] = separateUnits(vel_SI);
[data_Imperial, units_Imperial] = separateUnits(vel_Imperial);

パラシュート兵の速度は $t > 1$ のとき一定に近づきます。範囲 $0 \leq T \leq 2$ の速度変化をプロットすることで、速度がどのように終端速度に近づくかを示します。

subplot(1,2,1)
fplot(data_SI,[0 2])
title('Velocity in SI Units')
xlabel('Time in s')
ylabel(['Velocity in ' symunit2str(units_SI)])
subplot(1,2,2)
fplot(data_Imperial,[0 2])
title('Velocity in Imperial Units')
xlabel('Time in s')
ylabel(['Velocity in ' symunit2str(units_Imperial)])

Figure contains 2 axes objects. Axes object 1 with title Velocity in SI Units contains an object of type functionline. Axes object 2 with title Velocity in Imperial Units contains an object of type functionline.