空気抵抗がある場合の最遠投射角

こんにちは。Math。です。

パソコンのフォルダを整理していたら，数学などに関するメモやら PDF やらがいくつか出てきました。ほとんどはネタにもできないような内容であったり，途中で書き終えていたりと使えないんですが，中には使えそうなものもありました。

というわけで，今回はその中の一つ，「空気抵抗がはたらく場合における最遠投射角の式」について書いていきます。

運動方程式を立てる

図 1 のように，鉛直上向きに $y$ 軸をとり，水平方向に $x$ 軸をとります。このとき， $x$ 軸がちょうど地表を表しているとします。

そして，質量 $m$ の物体を原点 $\mathrm{O}$ から $x$ 軸正の向きに角度 $\theta$ ，速度 $v_0$ で投げたとします。ただし， $\displaystyle 0\lt\theta\lt\frac{\pi}{2}$ ， $v _ 0\gt 0$ です。

f:id:MathMaru:20210413154223p:plain:w400 — 図 1．斜方投射

さらに，この物体には重力加速度 $g$ と，速度に比例する空気抵抗（比例定数： $k\gt 0$ ）がはたらくものとします。

f:id:MathMaru:20210413155646p:plain:w250 — 図 2．運動方向と外力の向き

時刻 $t$ における物体の $x$ 座標， $y$ 座標をそれぞれ $x(t)$ ， $y(t)$ と表すとき，物体の運動方程式は次のようになります。なお， $\dot{x}$ は時刻 $t$ による $x$ の 1 階微分， $\ddot{x}$ は 2 階微分を表します¹。

運動方程式

$\left\{\begin{align} m\ddot{x} &= -k \dot{x}\\ m\ddot{y} &= -mg - k \dot{y} \end{align}\right.$

ここで，物体を時刻 $t=0$ に投げたとすれば，この運動方程式の初期条件は

$x(0) = y(0) = 0, \quad \dot{x}(0) = v_0 \cos\theta, \quad \dot{y}(0) = v_0 \sin\theta$

となります。以下では，この初期条件のもとで運動方程式を解いていきます。

運動方程式を解く

どちらも 2 階線形微分方程式ですが，速度を未知関数とすることで 1 階線形微分方程式に帰着できます。

$x$ 軸方向について

$v_x(t):=\dot{x}(t)$ とおくと， $x$ 軸方向の運動方程式は次のように書き換えられます。

$m\dot{v}_x = -kv_x$

これはすぐに解くことができて，初期条件 $v _ x(0)=\dot{x}(0)=v _ 0\cos\theta$ に注意すれば

$\displaystyle \dot{x}(t) = v_x(t) = v_0 e^{-\frac{k}{m}t} \cos\theta$

と求まります。さらに，初期条件 $x(0)=0$ のもとでこの両辺を $t$ で積分すれば

$\displaystyle x(t) = \frac{mv_0}{k} (1 - e^{-\frac{k}{m}t}) \cos\theta$

と求まります。

これより， $x$ 軸方向の変位 $x(t)$ と速度 $\dot{x}(t)$ のグラフの概形は図 3 のようになります。

f:id:MathMaru:20210414122907p:plain — 図 3． $x$ 軸方向の運動（左：変位，右：速度）

グラフを見てもらうと分かるように，時間が経つにつれて速度は 0 に近づいていきます。また，いくら時間が経とうとも距離 $\displaystyle\frac{mv_0}{k}\cos\theta$ までは飛ばないことも分かります。

$y$ 軸方向について

$x$ 軸方向に比べてやや複雑ですが，手順は全く同じです。先ほどと同様に $v_y(t):=\dot{y}(t)$ とおいて， $y$ 軸方向の運動方程式を次のように書き換えます。

$m\dot{v}_y = -mg - kv_y$

これは比較的簡単に解くことができて，初期条件 $v _ y(0)=\dot{y}(0)=v _ 0\sin\theta$ に注意すれば

$\displaystyle \dot{y}(t) = v_y(t) = \biggl(\frac{mg}{k} + v_0\sin\theta\biggr) e^{-\frac{k}{m}t} - \frac{mg}{k}$

と求まります。さらに，初期条件 $y(0)=0$ のもとで，この両辺を $t$ で積分すれば

$\displaystyle y(t) = \frac{m}{k}\biggl(\frac{mg}{k} + v_0\sin\theta\biggr) (1 - e^{-\frac{k}{m}t}) - \frac{mg}{k}t$

と求まります。

これより， $y$ 軸方向の変位 $y(t)$ と速度 $\dot{y}(t)$ のグラフの概形は図 4 のようになります。

f:id:MathMaru:20210414122913p:plain — 図 4． $y$ 軸方向の運動（左：変位，右：速度）

これもグラフから，時間の経過とともに落下距離が増えていくことや，落下速度が $\displaystyle\frac{mg}{k}$ を超えないことなどが分かります。

以上より，この物体の運動（の軌跡）は次のように表せます。

物体の運動

$\left\{\begin{align} x(t) &= \frac{mv_0}{k} (1 - e^{-\frac{k}{m}t}) \cos\theta\\ y(t) &= \frac{m}{k}\biggl(\frac{mg}{k} + v_0\sin\theta\biggr) (1 - e^{-\frac{k}{m}t}) - \frac{mg}{k}t \end{align}\right.\tag{1}$

最遠投射角

得られた式 $(1)$ を使って，物体の飛距離が最大となるような投射角（最遠投射角） $\theta_m$ を求めていきます。ただ，式 $(1)$ のままでは文字が多くてややこしいので，無次元量

$\displaystyle \hat{x} := \frac{kx}{mv_0}, \quad \hat{y} := \frac{k^2 y}{m^2 g}, \quad r := \frac{kv_0}{mg}, \quad \tau := \frac{k}{m}t$

を導入して，式 $(1)$ を次のように無次元化してあげます。

物体の運動（無次元化 ver.）

$\left\{\begin{align} \hat{x}(\tau) &= (1 - e^{-\tau}) \cos\theta\\ \hat{y}(\tau) &= (1 + r\sin\theta) (1 - e^{-\tau}) - \tau \end{align}\right.$

さて，最遠投射角を求めるためには，飛距離と投射角の関係を知る必要があります。そこで，まずは無次元化飛距離を求めていきます。

…と言いたいところなんですが，そう上手くはいきません。

飛距離を求めるためには， $\hat{y}(\tau)$ が再び 0 となる無次元化時刻 $\tau\gt 0$ を求める必要があります。しかし， $\hat{y}(\tau)$ の式には $\tau$ の 1 次関数と指数関数が混在していて，初等的に解くことができません。なので，少し工夫します。

どうするかというと， $\hat{x}(\tau)$ と $\hat{y}(\tau)$ の式を連立させて $\tau$ を消去します。 $\hat{x}(\tau)$ の式から

$\displaystyle 1 - e^{-\tau} = \frac{\hat{x}}{\cos\theta}, \quad \tau = - \log \biggl(1 - \frac{\hat{x}}{\cos\theta}\biggr)$

と変形できるので，それぞれを $\hat{y}(\tau)$ の式に代入することで

$\displaystyle \hat{y} = (1 + r\sin\theta)\frac{\hat{x}}{\cos\theta} + \log \biggl(1 - \frac{\hat{x}}{\cos\theta}\biggr)$

を得ます。よって，この式で $\hat{y}=0$ とすることで，無次元化飛距離 $\hat{x}$ の満たす方程式が得られます。

無次元化飛距離が満たす式

$\displaystyle 0 = (1 + r\sin\theta)\frac{\hat{x}}{\cos\theta} + \log \biggl(1 - \frac{\hat{x}}{\cos\theta}\biggr) \tag{2}$

ただし， $\hat{x}=0$ は自明な解なので，以下では $\hat{x}\gt 0$ とします。

このとき，投射角 $\theta$ を決めてあげると，飛距離 $\hat{x}$ は唯一つに定まります²。すなわち， $\hat{x}$ は $\theta$ の関数と考えることができます。特に，最遠投射角 $\theta _ m$ においては $\displaystyle\frac{d\hat{x}}{d\theta}(\theta _ m)=0$ を満たします。

このことに注意して，式 $(2)$ の両辺を $\theta$ で微分して $\theta=\theta _ m$ を代入すると

$\displaystyle 0 = r\hat{x} + \frac{\hat{x}\sin\theta_m}{\cos^2\theta_m}\Biggl[1 + r\sin\theta_m - \biggl(1 - \frac{\hat{x}}{\cos\theta_m}\biggr)^{-1}\Biggr]$

を得ます。 $\hat{x}\gt 0$ なので，この式はさらに

$\displaystyle 0 = r + \frac{\sin\theta_m}{\cos^2\theta_m}\Biggl[1 + r\sin\theta_m - \biggl(1 - \frac{\hat{x}}{\cos\theta_m}\biggr)^{-1}\Biggr]$

と書けます。よって，この式を $\hat{x}$ について解き，改めて式 $(2)$ に代入することで， $\theta_m$ の満たす式が得られます。

最遠投射角が満たす式

$\displaystyle \biggl(1 + \frac{r}{\sin\theta_m}\biggr)\biggl[\log\biggl(1 + \frac{r}{\sin\theta_m}\biggr) - 1\biggr] = r^2 - 1 \tag{3}$

$r=1$ の場合

式 $(3)$ の右辺が 0 となるので，この式を満たすのは

$\displaystyle 1 + \frac{1}{\sin\theta_m} = 0 \quad \mathrm{or} \quad \log\biggl(1 + \frac{1}{\sin\theta_m}\biggr) = 1$

のときです。ただし，左の式を満たす $\theta _ m$ は存在しないので，もし存在するならば右の式を満たします。そこで，右の式を $\sin\theta _ m$ について解くと

$\displaystyle \sin\theta_m = \frac{1}{e-1}$

となり，これを満たす $\theta _ m$ は存在します。実際， $\theta _ m\approx 35^\circ$ です。

$r \ne 1$ の場合

$\displaystyle\rho:=1+\frac{r}{\sin\theta _ m}$ とおいて式 $(3)$ を式変形していくと

$\displaystyle \frac{r^2-1}{\rho} e^{\frac{r^2-1}{\rho}} = \frac{r^2-1}{e}$

とできます³。

あいにく，ここから $\theta _ m$ について初等的に解くことはできませんが，Lambert の $W$ 関数⁴とよばれる特殊関数を用いても良いなら，次のようにして解くことができます。

（i） $r ^ 2-1\gt 0$ の場合

$\displaystyle \frac{r^2-1}{\rho} = W\biggl(\frac{r^2-1}{e}\biggr)$

と表せるので，さらに計算していくことで

$\displaystyle \sin\theta_m = r\biggl(\frac{r^2-1}{W((r^2-1)e^{-1})} - 1\biggr)^{-1}$

となります。 $r ^ 2-1\gt 0$ だったので， $W$ を主枝 $W_0$ に代えても同じです。

（ii） $r ^ 2-1\lt 0$ の場合

$\displaystyle -1 \lt -\frac{1}{\rho} \lt \frac{r^2-1}{\rho} \lt 0$

であることに注意すれば，主枝 $W _ 0$ を用いて

$\displaystyle \frac{r^2-1}{\rho} = W_0\biggl(\frac{r^2-1}{e}\biggr)$

と表せます。よって，全く同様に

$\displaystyle \sin\theta_m = r\biggl(\frac{r^2-1}{W_0((r^2-1)e^{-1})} - 1\biggr)^{-1}$

を得ます。

（i）と（ii）より， $r\ne 1$ の場合における最遠投射角 $\theta _ m$ は次のように書けます。

最遠投射角の式

$\displaystyle \theta_m = \sin^{-1} \Biggl[ r\biggl(\frac{r^2-1}{W_0((r^2-1)e^{-1})} - 1\biggr)^{-1} \Biggr]$

おわりに

かなり大変でしたが，空気抵抗がはたらく場合における最遠投射角の式を求めることができました。記事として書き直すために改めて自分でも計算していたんですが，過去の自分のモチベーションとか技巧に感心してしまいました。今の自分にできるとは到底思えないです。

ちなみに，いろんな数値で計算してみると分かりますが， $\theta _ m$ が 45 以外になることは普通にあります。すなわち，空気抵抗を考慮した場合，一番よく飛ぶ角度は 45 とは限らないということが分かります。

参考文献

ランベルトのW関数 - Wikipedia（閲覧日：2021/4/14）

$\dot{x}$ は $x$ 軸方向の速度， $\ddot{x}$ は $x$ 軸方向の加速度です。↩
$\hat{y}(\tau)$ のグラフを描くと， $\hat{y}=0$ となる $\tau$ が 2 つ（一つは自明な解 $\tau=0$ ）あることが分かります。そして， $\hat{x}$ と $\tau$ が一対一に対応していることから， $\hat{y}=0$ となる $\hat{x}$ も 2 つあります（一つは $\hat{x}=0$ ）。↩
式変形の際， $r ^ 2-1$ の符号に注意する必要があります。↩
大雑把に言うと， $y=xe ^ x$ の逆関数 $x=W(y)$ のことです。↩