航空機の横風突風応答

航空機の横風突風応答について説明する。

はじめに

この記事では、航空機が横方向の突風を受けたときに、横滑り角、ヨーイング角速度、ロール角速度、バンク角がどのように立ち上がるかについて説明する。

前記事では、ラダーのみ定常旋回において、ある速度と旋回率を保ちながら、エルロン舵角をゼロにしたときに、迎角、横滑り角、バンク角、ピッチ角、エレベータ舵角、ラダー舵角、推力がどのように決まるかを整理した。

航空機のラダーのみ定常旋回

ラダーのみを用いた定常旋回について説明する。

mtkbirdman.com

2026.07.11

この記事では横風突風応答性能として、横風突風を受けた直後から有限時間内にバンク角がどれだけ生じるかを考える。

横風突風応答では、以下の 2 つの経路でロール角速度が生じる。

\begin{align}
\beta_g &\to& C_{l\beta} &\to& \hat{p} \\
\beta_g &\to& C_{n\beta} &\to& \hat{r} &\to& C_{l\hat{r}} &\to& \hat{p} \end{align}

したがって、横風突風応答性能を評価するには、横風突風を受けてから有限時間内に \(\beta\)、\(\hat{p}\)、\(\hat{r}\)、\(\phi\) がどのように発達するかを見る必要がある。

この記事では、次の順序で横風突風応答の簡易指標を導く。

1. 横風突風を 1-cosine gust として定式化し、有効横滑り角への入り方を整理する。
2. 横滑り角、無次元ロール角速度、無次元ヨーイング角速度、バンク角を状態量とする線形横・方向モデルを作る。
3. 1-cosine gust に対する状態方程式の解を、行列指数関数で表す。
4. 行列指数関数と突風入力を Taylor 展開し、横風突風がロール角速度を生させる作用を取り出す。
5. 長スパン・大上反角航空機で副次的な項を省略し、主要なロール生成経路を整理する。
6. 安定微係数と慣性モーメントで表される簡易横風突風ロール指数を定義する。
7. 簡易横風突風ロール指数を、6 自由度有限時間応答と比較する。

それではいってみよう。

横風突風の定式化

定常飛行速度を \(V\)、機体の横方向速度擾乱を \(v\) とすると微小擾乱では、横滑り角は次で表される。

\begin{align} \beta \simeq \frac{v}{V} \end{align}

ここでは、横風突風速度 \(U_g\) を、機体軸 \(+y_b\) 方向、すなわち右から吹いてくる方向を正とする大気速度成分として定義する。

このとき、機体が感じる有効な横方向相対速度は \(v+U_g\) である。

横風突風による外乱横滑り角は、

\begin{align}
\beta_g &=\tan^{-1}\left(\frac{U_g}{V}\right)
\simeq \frac{U_g}{V}
\end{align}

である。

したがって、微小擾乱の範囲では、有効横滑り角は

\begin{align} \beta_a \simeq \beta +\beta_g \end{align}

である。

このように、横風突風は空力係数の中で \(\beta_a\) を通じて入る。

1-cosine gust の定式化

今回は、離散突風の代表形として、14 CFR 25.341 に示される 1-cosine gust の形を基準にしてみる。

横方向の 1-cosine gust は、突風場内の進行距離 \(s\) を用いて

\begin{align} U_g(s)=\frac{U_{ds}}{2}\left[1-\cos\left(\frac{\pi s}{H}\right)\right],\quad 0 \le s \le 2H \end{align}

とする。

ここで、本記事の運動方程式に用いる \(U_{ds}\) は、解析高度における有次元の最大突風速度である。正負の符号により、反対方向の横風突風を表すことができる。

なお、14 CFR 25.341 に規定される \(U_{ds}\) は等価対気速度であるため、その値を運動方程式へ入力する場合は、解析高度における真突風速度へ変換する必要がある。

\begin{align} U_{ds,\mathrm{TAS}} =U_{ds,\mathrm{EAS}} \sqrt{ \frac{\rho_0}{\rho} } \end{align}

以降の \(U_{ds}\) は、特に断らない限り、運動方程式に用いる真突風速度 \(U_{ds,\mathrm{TAS}}\) を表す。

\(H\) は、突風がゼロから最大値に達するまでに航空機が飛行経路方向へ進む距離である。

したがって、1-cosine gust 全体の空間長さは \(2H\) である。

機体が速度 \(V\) で突風場を通過するとすれば、

\begin{align} s \simeq Vt \end{align}

である。

したがって、時間領域では

\begin{align} U_g(t)=\frac{U_{ds}}{2}\left[1-\cos\left(\frac{\pi Vt}{H}\right)\right],\quad 0 \le t \le \frac{2H}{V} \end{align}

となる。

外乱横滑り角は

\begin{align} \beta_g(t)\simeq\frac{U_{ds}}{2V}\left[1-\cos\left(\frac{\pi Vt}{H}\right)\right],\quad 0 \le t \le \frac{2H}{V} \end{align}

である。

無次元時間を

\begin{align} \tau=\frac{2V}{b}t \end{align}

とする。

このとき、

\begin{align} s&=Vt=\frac{b}{2}\tau \end{align}

なので、1-cosine gust は

\begin{align} \beta_g(\tau)\simeq\frac{U_{ds}}{2V}\left[1-\cos\left(\frac{\pi b}{2H}\tau\right)\right],\quad 0 \le \tau \le \frac{4H}{b} \end{align}

である。

以降、

\begin{align}
\beta_{g0}&=\frac{U_{ds}}{2V}, \qquad \omega=\frac{\pi b}{2H}
\end{align}

と置く。

ここで、\(\beta_{g0}\) は 1-cosine gust の半振幅係数であり、外乱横滑り角の最大値は

\begin{align} \beta_{g,\max}\simeq2\beta_{g0}=\frac{U_{ds}}{V} \end{align}

である。

小角近似の範囲では、横風突風入力は

\begin{align} \beta_g(\tau)=\beta_{g0}\left[1-\cos(\omega\tau)\right] \end{align}

である。

突風外では

\begin{align} \beta_g(\tau)=0 \end{align}

として扱う。

1-cosine gust の通過終了時刻は、無次元時間で

\begin{align} \tau_e&=\frac{2V}{b}\frac{2H}{V}=\frac{4H}{b} \end{align}

である。

また、

\begin{align} \omega \tau_e&=\frac{\pi b}{2H}\frac{4H}{b}=2\pi \end{align}

である。

線形横・方向モデルの定式化

機体軸は、

\begin{align} +x_b &: \text{前方}, \\ +y_b &: \text{右方}, \\ +z_b &: \text{下方} \end{align}

とする。

また、慣性テンソルは、

\begin{align} \mathbf{I}_b =\begin{bmatrix} I_x & 0 & -I_{xz} \\ 0 & I_y & 0 \\ -I_{xz} & 0 & I_z \end{bmatrix} \end{align}

とする。

状態量

横・方向の線形状態量を次で定義する。

ここで、\(\beta\)、\(\hat{p}\)、\(\hat{r}\)、\(\phi\)、\(\beta_g\) は、特に断らない限り基準飛行状態からの摂動量として扱う。

\begin{align} \mathbf{x} = \begin{bmatrix} \beta \\ \hat{p} \\ \hat{r} \\ \phi \end{bmatrix} \end{align}

ここで、無次元角速度は以下の式で表される。

\begin{align} \hat{p}&=\frac{pb}{2V}, \\ \hat{r}&=\frac{rb}{2V} \end{align}

無次元時間 \(\tau=\dfrac{2V}{b}t\) の微分は以下のようになる。

\begin{align} \frac{d}{d\tau} =\frac{b}{2V}\frac{d}{dt} \end{align}

以降、\(t\) に関する微分をドット、\(\tau\) に関する微分をプライムで表す。

横滑り角の式

横方向の並進運動は、機体軸速度成分を用いると次で表される。

\begin{align} \dot v=pw-ru+\frac{Y}{m}+g\sin\phi\cos\theta \end{align}

小迎角・小横滑り・小バンク近似を用い、\(u\simeq V\)、\(v\simeq V\beta\)、\(w\simeq V\alpha_0\) とすれば、

\begin{align} \dot{\beta}\simeq\frac{Y}{mV}-r+\alpha_0p+\frac{g}{V}\phi \end{align}

となる。

横力を

\begin{align} Y=q_\infty S C_Y \end{align}

とし、横力係数を

\begin{align} C_Y=C_{Y\beta}(\beta+\beta_g)+C_{Y\hat{p}}\hat{p}+C_{Y\hat{r}}\hat{r} \end{align}

と線形化する。

VSPAERO の .stab を用いる場合、本文の \(C_Y\) とその微係数には、\(C_L\)、\(C_D\)、\(C_S\) から機体軸へ変換した横力係数を用いる。

\begin{align} C_Y =-C_D\sin\beta +C_S\cos\beta \end{align}

無次元時間で書くと、横滑り角の式は次のようになる。

\begin{align} \beta' =\mu_y\left[C_{Y\beta}(\beta+\beta_g)+C_{Y\hat{p}}\hat{p}+C_{Y\hat{r}}\hat{r}\right]+\alpha_0\hat{p}-\hat{r}+\lambda_g\phi \end{align}

ここで、

\begin{align} \mu_y&=\frac{q_\infty S b}{2mV^2}=\frac{\rho S b}{4m} \\ \lambda_g&=\frac{gb}{2V^2} \end{align}

である。

ロール角速度とヨーイング角速度の式

慣性乗積 \(I_{xz}\) を残し、イナーシャカップリング項を省略すると、ロール・ヨーの線形回転運動は次で表される。

\begin{align} I_x\dot{p}-I_{xz}\dot{r}&=L, \\ -I_{xz}\dot{p}+I_z\dot{r}&=N \end{align}

である。

ここで、

\begin{align} L&=q_\infty S b C_l, \\ N&=q_\infty S b C_n \end{align}

である。

空力モーメント係数を

\begin{align} C_l&=C_{l\beta}(\beta+\beta_g)+C_{l\hat{p}}\hat{p}+C_{l\hat{r}}\hat{r}, \\ C_n&=C_{n\beta}(\beta+\beta_g)+C_{n\hat{p}}\hat{p}+C_{n\hat{r}}\hat{r} \end{align}

とする。

ここで、慣性係数を

\begin{align} \mu_I =\frac{q_\infty S b^3}{4V^2\left(I_xI_z-I_{xz}^2\right)} \end{align}

と定義する。

このとき、無次元ロール角速度と無次元ヨーイング角速度の式は

\begin{align} \hat{p}' &=\mu_I\left(I_zC_l+I_{xz}C_n\right) \\ \hat{r}' &=\mu_I\left(I_{xz}C_l+I_xC_n\right) \end{align}

である。

\(I_{xz}=0\) とすれば、

\begin{align} \mu_x&=\frac{q_\infty S b^3}{4I_xV^2}, \\ \mu_z&=\frac{q_\infty S b^3}{4I_zV^2} \end{align}

を用いて、

\begin{align} \hat{p}'&=\mu_x C_l, \\ \hat{r}'&=\mu_z C_n \end{align}

である。

バンク角の式

\(\phi \simeq 0\) まわりの小角近似では、バンク角の運動学は

\begin{align} \dot{\phi}=p+\tan\theta_0\,r \end{align}

である。

無次元時間では、

\begin{align} \phi'=\hat{p}+\tan\theta_0\,\hat{r} \end{align}

となる。

水平飛行に近ければ \(\theta_0 \simeq 0\) なので、

\begin{align} \phi' \simeq \hat{p} \end{align}

である。

状態行列と入力ベクトル

以上をまとめると、横・方向線形モデルは次の状態方程式で表せる。

\begin{align} \mathbf{x}'=\mathbf{A}\mathbf{x}+\mathbf{B}_g\beta_g(\tau) \end{align}

状態行列 \(\mathbf{A}\) は、次の式で表される。

\begin{align} \mathbf{A} =\begin{bmatrix} \mu_y C_{Y\beta} & \mu_y C_{Y\hat{p}}+\alpha_0 & \mu_y C_{Y\hat{r}}-1 & \lambda_g \\ \mu_I\left(I_zC_{l\beta}+I_{xz}C_{n\beta}\right) & \mu_I\left(I_zC_{l\hat{p}}+I_{xz}C_{n\hat{p}}\right) & \mu_I\left(I_zC_{l\hat{r}}+I_{xz}C_{n\hat{r}}\right) & 0 \\ \mu_I\left(I_{xz}C_{l\beta}+I_xC_{n\beta}\right) & \mu_I\left(I_{xz}C_{l\hat{p}}+I_xC_{n\hat{p}}\right) & \mu_I\left(I_{xz}C_{l\hat{r}}+I_xC_{n\hat{r}}\right) & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix} \end{align}

横風突風は有効横滑り角 \(\beta+\beta_g\) を通じて入るため、突風入力ベクトルは

\begin{align}
\mathbf{B}_g =\begin{bmatrix} \mu_y C_{Y\beta} \\ \mu_I\left(I_zC_{l\beta}+I_{xz}C_{n\beta}\right) \\ \mu_I\left(I_{xz}C_{l\beta}+I_xC_{n\beta}\right) \\ 0 \end{bmatrix} \end{align}

状態行列 \(\mathbf{A}\) は、横滑り角、ロール角速度、ヨーイング角速度、バンク角の相互作用を、入力ベクトル \(\mathbf{B}_g\) は、外乱横滑り角が各状態へ入る経路を表す。

以降は、状態行列と入力ベクトルの成分を次のように表す。

\begin{align} \mathbf{A} =\begin{bmatrix} a_{\beta\beta} & a_{\beta p} & a_{\beta r} & a_{\beta\phi} \\ a_{p\beta} & a_{pp} & a_{pr} & 0 \\ a_{r\beta} & a_{rp} & a_{rr} & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}, \qquad \mathbf{B}_g =\begin{bmatrix} B_{g,\beta} \\ B_{g,p} \\ B_{g,r} \\ 0 \end{bmatrix} \end{align}

ここで、

\begin{align} a_{\beta\beta}&=\mu_y C_{Y\beta}, \\ a_{\beta p}&=\mu_y C_{Y\hat{p}}+\alpha_0, \\ a_{\beta r}&=\mu_y C_{Y\hat{r}}-1, \\ a_{\beta\phi}&=\lambda_g, \\ a_{p\beta}&=\mu_I\left(I_zC_{l\beta}+I_{xz}C_{n\beta}\right), \\ a_{pp}&=\mu_I\left(I_zC_{l\hat{p}}+I_{xz}C_{n\hat{p}}\right), \\ a_{pr}&=\mu_I\left(I_zC_{l\hat{r}}+I_{xz}C_{n\hat{r}}\right), \\ a_{r\beta}&=\mu_I\left(I_{xz}C_{l\beta}+I_xC_{n\beta}\right), \\ a_{rp}&=\mu_I\left(I_{xz}C_{l\hat{p}}+I_xC_{n\hat{p}}\right), \\ a_{rr}&=\mu_I\left(I_{xz}C_{l\hat{r}}+I_xC_{n\hat{r}}\right) \end{align}

また、

\begin{align}
B_{g,\beta}&=a_{\beta\beta}, \\
B_{g,p}&=a_{p\beta},
\\ B_{g,r}&=a_{r\beta}
\end{align}

である。

1-cosine gust に対する線形応答

線形状態方程式

\begin{align} \mathbf{x}'=\mathbf{A}\mathbf{x}+\mathbf{B}_g\beta_g(\tau) \end{align}

において、時刻 \(\tau=0\) で 1-cosine gust を入力したときの時間応答の解を求める。

初期状態は以下の式で表される。

\begin{align} \mathbf{x}(0)=\mathbf{0} \end{align}

行列指数関数を用いると、状態量は

\begin{align} \mathbf{x}(\tau) =\int_0^\tau \exp\left[\mathbf{A}(\tau-\sigma)\right]\mathbf{B}_g\beta_g(\sigma)\,d\sigma \end{align}

となる。

外乱横滑り角の最大値を

\begin{align} \beta_{g,\max} \simeq \frac{U_{ds}}{V} \end{align}

とすれば、1-cosine gust は

\begin{align} \beta_g(\sigma) =\frac{\beta_{g,\max}}{2} \left[ 1- \cos\left( \frac{\pi b}{2H}\sigma \right) \right] \end{align}

と書ける。

したがって、

\begin{align} \mathbf{x}(\tau) =\frac{\beta_{g,\max}}{2}\int_0^\tau\exp\left[\mathbf{A}(\tau-\sigma)\right]\mathbf{B}_g\left[1-\cos\left(\frac{\pi b}{2H}\sigma\right)\right]d\sigma \end{align}

である。

さらに、\(s=\tau-\sigma\) と置くと、\(ds=-d\sigma\) であり、積分範囲は \(\sigma=0\) で \(s=\tau\)、\(\sigma=\tau\) で \(s=0\) となる。

よって、

\begin{align} \mathbf{x}(\tau) &=\frac{\beta_{g,\max}}{2} \int_\tau^0 \exp(\mathbf{A}s) \mathbf{B}_g \left[1-\cos\left(\frac{\pi b}{2H}(\tau-s)\right)\right](-ds) \\ &=\frac{\beta_{g,\max}}{2}\int_0^\tau\exp(\mathbf{A}s)\mathbf{B}_g\left[1-\cos\left(\frac{\pi b}{2H}(\tau-s)\right)\right]ds \end{align}

となる。

ロール角速度応答とバンク角応答

ロール角速度成分とバンク角成分を取り出すベクトルを次で定義する。

\begin{align} \mathbf{e}_p =\begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \end{bmatrix}, \qquad \mathbf{e}_\phi =\begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \end{bmatrix} \end{align}

ここで、\(\mathbf{e}_p\) は無次元ロール角速度 \(\hat{p}\) を取り出し、\(\mathbf{e}_\phi\) はバンク角 \(\phi\) を取り出す。

横風突風に対するロール角速度応答とバンク角応答は以下の式で表される。

\begin{align} \hat{p}_{\mathrm{lin}}(\tau) &=\frac{\beta_{g,\max}}{2} \mathbf{e}_p^{\mathrm T} \int_0^\tau \exp\left[ \mathbf{A}(\tau-\sigma) \right] \mathbf{B}_g \left[1-\cos\left(\frac{\pi b}{2H}\sigma\right)\right]d\sigma, \\ \Delta\phi_{\mathrm{lin}}(\tau) &=\frac{\beta_{g,\max}}{2} \mathbf{e}_\phi^{\mathrm T} \int_0^\tau \exp\left[\mathbf{A}(\tau-\sigma)\right] \mathbf{B}_g\left[1-\cos\left(\frac{\pi b}{2H}\sigma\right)\right]d\sigma \end{align}

したがって、ピーク外乱横滑り角あたりのロール角速度応答とバンク角応答は、以下のようになる。

\begin{align} \frac{\hat{p}_{\mathrm{lin}}(\tau)}{\beta_{g,\max}} &=\frac{1}{2} \mathbf{e}_p^{\mathrm T} \int_0^\tau \exp\left[ \mathbf{A}(\tau-\sigma) \right] \mathbf{B}_g \left[ 1- \cos\left( \frac{\pi b}{2H}\sigma \right) \right] d\sigma, \\ \frac{\Delta\phi_{\mathrm{lin}}(\tau)}{\beta_{g,\max}} &=\frac{1}{2} \mathbf{e}_\phi^{\mathrm T} \int_0^\tau \exp\left[ \mathbf{A}(\tau-\sigma) \right] \mathbf{B}_g \left[ 1- \cos\left( \frac{\pi b}{2H}\sigma \right) \right] d\sigma \end{align}

突風通過終了時刻は、

\begin{align} \tau_e =\frac{4H}{b} \end{align}

であり、

\begin{align} \cos\left[ \frac{\pi b}{2H} (\tau_e-s) \right] =\cos\left( 2\pi-\frac{\pi b}{2H}s \right) =\cos\left( \frac{\pi b}{2H}s \right) \end{align}

である。

したがって、突風通過終了時のピーク外乱横滑り角あたりの応答は、

\begin{align} \frac{\hat{p}_{\mathrm{lin}}(\tau_e)}{\beta_{g,\max}} &=\frac{1}{2} \mathbf{e}_p^{\mathrm T} \int_0^{4H/b} \exp(\mathbf{A}s) \mathbf{B}_g \left[ 1- \cos\left( \frac{\pi b}{2H}s \right) \right] ds, \\ \frac{\Delta\phi_{\mathrm{lin}}(\tau_e)}{\beta_{g,\max}} &=\frac{1}{2} \mathbf{e}_\phi^{\mathrm T} \int_0^{4H/b} \exp(\mathbf{A}s) \mathbf{B}_g \left[ 1- \cos\left( \frac{\pi b}{2H}s \right) \right] ds \end{align}

となる。

Taylor 展開による代数近似

線形有限時間応答は行列指数を含み、このままではよくわからないので、Taylor 展開を使って多項式の形に書き下してみる。

行列指数関数を Taylor 展開すると、次のようになる。

\begin{align} \exp\left[ \mathbf{A}(\tau-\sigma) \right] =\sum_{k=0}^{\infty} \frac{ \mathbf{A}^k(\tau-\sigma)^k }{ k! } \end{align}

である。

また、ピーク外乱横滑り角で正規化した 1-cosine gust は、

\begin{align} \frac{\beta_g(\sigma)}{\beta_{g,\max}} &=\frac{1}{2} \left[ 1- \cos\left( \frac{\pi b}{2H}\sigma \right) \right] \\\\
&=\frac{\pi^2b^2}{16H^2}\sigma^2 -\frac{\pi^4b^4}{768H^4}\sigma^4 +\frac{\pi^6b^6}{92160H^6}\sigma^6 -\cdots \end{align}

である。

これをロール角速度応答とバンク角応答に代入すると

\begin{align} \frac{\mathbf{x}(\tau)}{\beta_{g,\max}} ={}& \frac{\pi^2b^2}{48H^2} \mathbf{B}_g\tau^3 +\frac{\pi^2b^2}{192H^2} \mathbf{A}\mathbf{B}_g\tau^4 +\left[ \frac{\pi^2b^2}{960H^2} \mathbf{A}^2\mathbf{B}_g -\frac{\pi^4b^4}{3840H^4} \mathbf{B}_g \right]\tau^5 +\left[ \frac{\pi^2b^2}{5760H^2} \mathbf{A}^3\mathbf{B}_g -\frac{\pi^4b^4}{23040H^4} \mathbf{A}\mathbf{B}_g \right]\tau^6 +\cdots \end{align}

である。

1-cosine gust は \(\tau^2\) から立ち上がるため、横滑り角、無次元ロール角速度、無次元ヨーイング角速度は \(\tau^3\) から立ち上がる。

バンク角は無次元ロール角速度を積分した状態量なので、さらに 1 次遅れて \(\tau^4\) から立ち上がる。

Taylor 係数とロール生成経路

ここで、

\begin{align} \mathbf{e}_\phi^{\mathrm T}\mathbf{A} &=\begin{bmatrix} 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a_{\beta\beta} & a_{\beta p} & a_{\beta r} & a_{\beta\phi} \\ a_{p\beta} & a_{pp} & a_{pr} & 0 \\ a_{r\beta} & a_{rp} & a_{rr} & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}\\\\ &=\begin{bmatrix} 0 & 1 & 0 & 0 \end{bmatrix} =\mathbf{e}_p^{\mathrm T} \end{align}

が成り立つ。

そこで、ロール角速度応答とバンク角応答に共通する係数列を、

\begin{align} K_n &=\mathbf{e}_\phi^{\mathrm T}\mathbf{A}^n\mathbf{B}_g =\mathbf{e}_p^{\mathrm T}\mathbf{A}^{n-1}\mathbf{B}_g, \qquad n=1,2,3,\cdots, \\\\ K_0 &=\mathbf{e}_\phi^{\mathrm T}\mathbf{B}_g =0 \end{align}

と定義すると、ロール角速度応答とバンク角応答の式は以下のように整理できる。

\begin{align} \frac{\hat{p}_{\mathrm{lin}}(\tau)}{\beta_{g,\max}} ={}& \frac{\pi^2b^2}{48H^2}K_1\tau^3 +\frac{\pi^2b^2}{192H^2}K_2\tau^4 +\left(\frac{\pi^2b^2}{960H^2}K_3-\frac{\pi^4b^4}{3840H^4}K_1\right)\tau^5 +\cdots, \\ \frac{\Delta\phi_{\mathrm{lin}}(\tau)}{\beta_{g,\max}} ={}& \frac{\pi^2b^2}{192H^2}K_1\tau^4 +\frac{\pi^2b^2}{960H^2}K_2\tau^5 +\left(\frac{\pi^2b^2}{5760H^2}K_3-\frac{\pi^4b^4}{23040H^4}K_1\right)\tau^6 +\cdots \end{align}

となる。

各係数は、

\begin{align} K_1 &=B_{g,p}, \\
K_2 &=a_{p\beta}B_{g,\beta}+a_{pp}B_{g,p}+a_{pr}B_{g,r}, \\
\end{align}

である。

ここで、\(K_1\) は、

\begin{align} K_1 =\mathbf{e}_p^{\mathrm T}\mathbf{B}_g \end{align}

であり、突風入力ベクトルが最初に無次元ロール角速度へ入る直接的な作用を表す。

すなわち、

\begin{align} \beta_g \rightarrow C_{l\beta} \rightarrow \hat{p} \end{align}

という経路である。

\(K_2\) は、

\begin{align} K_2 =\mathbf{e}_p^{\mathrm T}\mathbf{A}\mathbf{B}_g \end{align}

であり、突風入力によって最初の状態変化が生じ、その状態変化が無次元ロール角速度へ入る作用を表す。

代表的には、

\begin{align}
\beta_g &\rightarrow C_{n\beta} \rightarrow \hat{r} \rightarrow C_{l\hat{r}} \rightarrow \hat{p}
\end{align}

という経路である。

簡易横風突風ロール指数の導出

主要なロール生成経路を読み取るため、次の近似を置く。

\begin{align} I_{xz} &\simeq 0, \\ \alpha_0 &\simeq 0, \\ \theta_0 &\simeq 0, \\ C_{Y\hat{p}} &\simeq 0, \\ C_{Y\hat{r}} &\simeq 0, \\ C_{n\hat{p}} &\simeq 0 \end{align}

また、

\begin{align} q_\infty =\frac{1}{2}\rho V^2 \end{align}

を用いる。

入力ベクトルが直接作る作用

\(I_{xz}\simeq0\) とすれば、突風入力ベクトルの各成分は、

\begin{align}
B_{g,\beta} &\simeq \frac{\rho S b}{4m}C_{Y\beta}, \\
B_{g,p} &\simeq \frac{\rho S b^3}{8I_x}C_{l\beta}, \\
B_{g,r} &\simeq \frac{\rho S b^3}{8I_z}C_{n\beta} \end{align}

となる。

したがって、入力ベクトルが直接作るロール角速度応答の係数は、

\begin{align}
K_{1,\mathrm{simple}} =\frac{\rho S b^3}{8I_x}C_{l\beta}
\end{align}

である。

これは、

\begin{align} \beta_g \rightarrow C_{l\beta} \rightarrow \hat{p} \end{align}

というように、突風による有効横滑り角の変化が \(C_{l\beta}\) を通じて直接ロールへ入る経路を表す。

状態行列を 1 回介する作用

同じ近似のもとで、ロール角速度方程式に関係する状態行列の成分は、

\begin{align} a_{p\beta} &\simeq \frac{\rho S b^3}{8I_x}C_{l\beta}, \\ a_{pp} &\simeq \frac{\rho S b^3}{8I_x}C_{l\hat{p}}, \\ a_{pr} &\simeq \frac{\rho S b^3}{8I_x}C_{l\hat{r}} \end{align}

である。

したがって、

\begin{align}
K_{2,\mathrm{simple}} ={}& a_{p\beta}B_{g,\beta} +a_{pp}B_{g,p} +a_{pr}B_{g,r} \\\\
={}&
\frac{\rho^2S^2b^4}{32mI_x} C_{l\beta}C_{Y\beta}
+\frac{\rho^2S^2b^6}{64I_x^2} C_{l\hat{p}}C_{l\beta}
+\frac{\rho^2S^2b^6}{64I_xI_z} C_{l\hat{r}}C_{n\beta} \end{align}

となる。

第 1 項は

• 横力を介して機体側の横滑りが生じ、
• 上反角効果でロールする経路

第 2 項は

• 突風によって生じたロール角速度に
• ロール減衰 \(C_{l\hat{p}}\) が作用する経路

第 3 項は

• 方向静安定で生じたヨーイング角速度が
• ヨーイング角速度によるロールを通じて無次元ロール角速度へ入る経路

である。

簡易横風突風ロール指数

以上より、初期のバンクの生成のみにかかわる主要な経路を、 Taylor 展開の次数をいったん脇に置いてまとめると、以下のようになる。

\begin{align}
K_{\mathrm{gust,path}} \simeq
+\frac{\rho S b^3}{8I_x} C_{l\beta}
+\frac{\rho^2S^2b^6}{64I_xI_z} C_{l\hat{r}}C_{n\beta}
\end{align}

ここで、第 1 項は、突風による有効横滑り角が \(C_{l\beta}\) を通じて直接ロールへ入る経路、第 2 項は、方向静安定で生じたヨーイング角速度が、\(C_{l\hat{r}}\) を通じてロール角速度へ入る経路である。

\(K_{\mathrm{gust,path}}<0\) のときは、右から吹いてくる横風突風 \(\beta_{g}>0\) に対して左のローリング角速度 \(\hat{p}<0\) が生じるということで、横風に煽られやすいということになる。

上式は次のようにもまとめられる。

\begin{align} K_{\mathrm{gust,path}} =\frac{\rho^2S^2b^6}{64I_xI_z} \left( \frac{8I_z}{\rho S b^3}C_{l\beta} +C_{l\hat{r}}C_{n\beta} \right) \end{align}

したがって、主指標は次の形になる。

\begin{align}
K_{\mathrm{gust,roll}} \simeq
\frac{8I_z}{\rho S b^3}C_{l\beta}+C_{l\hat{r}}C_{n\beta}
\end{align}

ここで、通常の航空機では \(C_{l\beta}<0\)、\(C_{l\hat{r}}C_{n\beta}>0\) なので、横風突風に対して、\(|C_{l\beta}|\)が大きければ大きいほど煽られやすく、\(C_{n\beta}\)が大きければ大きいほど煽られにくいということがわかる。

VvGammaChart.py の postprocess\_vv\_gamma\_cases() では、この値を crosswind\_gust\_roll\_index として出力する。

6 自由度有限時間応答との比較

最後に、ここまでに得られた簡易横風突風ロール指数が、有限時間の 6 自由度応答とどの程度対応するかを確認する。

比較には、VSPAERO の .stab から得た線形空力モデルを用いて、非線形剛体 6 自由度運動方程式を数値積分した結果を用いる。

1-cosine gust を入力し、初期バンクを \(\phi_0\) から評価時間内に生じたバンク角変化の最大絶対値を用い、6 自由度有限時間横風突風応答指標を次で定義する。

\begin{align} |\Delta\phi|_{\max} =\max_{t_0\le t\le t_f} \left| \phi(t)-\phi_0 \right| \end{align}

この指標は、突風通過中と突風通過後を含む評価時間内で、機体がどこまでバンクしたかを表す量である。

6 自由度剛体運動シミュレーションの詳細は、以下の記事に譲る。

航空機の横風突風応答の 6 自由度剛体運動をシミュレーションする Python スクリプト

今回の比較では、垂直尾翼容積比 \(V_v\) と等価上反角 \(\Gamma_{\mathrm{eff}}\) をパラメトリックに変化させ、次の条件で横風突風応答を計算した。

\begin{align} V&=10\ \mathrm{m/s}, \\ U_{ds}&=3\ \mathrm{m/s}, \\ H&=9.144\ \mathrm{m}, \\ I_z&=1000\ \mathrm{kg\,m^2}, \\ t_0&=0\ \mathrm{s}, \\ t_f&=3\ \mathrm{s} \end{align}

突風通過終了時刻は、

\begin{align}
t_e &=\frac{2H}{V} =1.8288\ \mathrm{s}
\end{align}

である。

したがって、今回の評価時間には、突風通過中の応答に加えて、突風通過後の過渡応答も含まれる。

全 231 ケースのうち、安定微係数の後処理に成功したのは 228 ケースであり、そのすべてで横風突風に対する 6 自由度応答が得られた。

最大バンク角変化は、

\begin{align} 0.081\ \mathrm{deg} \le |\Delta\phi|_{\max} \le 3.058\ \mathrm{deg} \end{align}

の範囲にある。

\begin{align} |\Delta\phi|_{\max,\mathrm{mean}} \simeq 1.264\ \mathrm{deg} \end{align}

である。

最大値は、

\begin{align} V_v&=0.001, \\ \Gamma_{\mathrm{eff}} &=9.693\ \mathrm{deg}, \\ |\Delta\phi|_{\max} &=3.058\ \mathrm{deg}, \\ t\left(|\Delta\phi|_{\max}\right) &=1.616\ \mathrm{s} \end{align}

で得られた。

簡易横風突風指数と 6 自由度有限時間横風突風応答指標の Pearson の相関係数は、

\begin{align} r\left( K_{\mathrm{gust,roll}}, |\Delta\phi|_{\max} \right) \simeq 0.939 \end{align}

である。

したがって、

\begin{align}
K_{\mathrm{gust,roll}} =
\frac{8I_z}{\rho S b^3}C_{l\beta} +C_{l\hat{r}}C_{n\beta}
\end{align}

は、今回の解析範囲では、垂直尾翼容積比と上反角の組合せによる横風突風時のロール応答の主傾向を把握する簡易指標として機能している。

設計初期には crosswind\_gust\_roll\_index で応答の大きい領域を確認し、その後で 6 自由度時刻歴、最大応答時刻、突風通過後の残留応答を確認する流れが使いやすい。

おわりに

航空機の横風突風応答について説明した。