【流束ヤコビ行列の導出】1次元圧縮性性オイラー方程式

こんにちは（@t_kun_kamakiri）
本記事では、流体力学の基礎式を簡単にしたオイラー方程式に関する基礎理論のうちの流束ヤコビ行列の導出を解説します。
こちらは計算力学技術者の熱流体1級（単相流）で問われる知識です。
導出過程も含め物理的意味を押さえておくと良いです。

まず、本記事で示したい内容は以下です。

ここにタイトルを入力

システム方程式から流束ヤコビ行列を導出する

システム方程式

\begin{align*}
\frac{\partial \boldsymbol{Q}}{\partial t} + \frac{\partial \boldsymbol{F}}{\partial x} = 0
\end{align*}

流束ヤコビ行列

\begin{align*}
\boldsymbol{A}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}=
\begin{pmatrix}
0 & 1 & 0 \\
-\frac{3-\gamma}{2}u^2 & (3-\gamma)u & \gamma -1 \\
\big(\frac{\gamma-1}{2}-h\big)u & h-\big(\gamma-1\big)u^2 & \gamma u
\end{pmatrix}
\end{align*}

$\boldsymbol{A}$を導出するのが本記事の内容になります。

システム方程式とは
流束ヤコビ行列の導入
理想気体の状態方程式
流束ヤコビ行列の成分の導出
まとめ
参考書
計算力学技術者資格のための問題アプリ

システム方程式とは

圧縮性流れを扱う場合、レイノルズ数$Re=\frac{UL}{\nu}$が大きいことが多く、粘性項目を無視したオイラー方程式での数値計算法の解法を基本として考えられています。

以下に示すのが1次元圧縮性の非粘性に対するオイラー方程式です。

質量保存則

\begin{align*}
\frac{\partial \rho}{\partial t}+\frac{\partial \rho u}{\partial x}=0\tag{1}
\end{align*}

運動量保存則

\begin{align*}
\frac{\partial \rho u}{\partial t}+\frac{\partial }{\partial x}(p+\rho u^2)=0\tag{2}
\end{align*}

エネルギー保存則

\begin{align*}
\frac{\partial e}{\partial t}+\frac{\partial }{\partial x}(e+p)u=0\tag{3}
\end{align*}

$\rho$：密度
$u$：流速
$p$：圧力
$e$：単位体積当たりの気体の全エネルギー

これらをまとめると以下になります。

\begin{align*}
\frac{\partial \boldsymbol{Q}}{\partial t} + \frac{\partial \boldsymbol{F}}{\partial x} = 0\tag{4}
\end{align*}

ここで、

\begin{align*}
\boldsymbol{Q} =
\begin{bmatrix}
\rho \\
\rho u \\
e \\
\end{bmatrix}
=
\begin{bmatrix}
q_1 \\
q_2 \\
q_3 \\
\end{bmatrix}, \quad
\boldsymbol{F} =
\begin{bmatrix}
\rho u \\
p+\rho u^2 \\
(e+p)u \\
\end{bmatrix}=
\begin{bmatrix}
f_1 \\
f_2 \\
f_3 \\
\end{bmatrix}
\end{align*}

スカラー輸送方程式（移流方程式など）に対して、(4)を総称してシステム方程式と呼ぶことがあります。流速がマッハ数に対して大きくなる場合、圧縮性流れとして現象をモデル化する必要が出てきますが、レイノルズ数も大きくなるため、粘性項を無視したオイラー方程式の解法を基本として計算方法が展開されているということが大事な点です。

流束ヤコビ行列の導入

(4)の特徴を理解するために以下のように式を変換します。
以下$\boldsymbol{A}$のような行列を導入し、行列の対角化を行うと現象を理解しやすくなるというメリットがありますが、本記事ではそこまでは踏み込みません。

\begin{align*}
\frac{\partial \boldsymbol{Q}}{\partial t} + \frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}\frac{\partial \boldsymbol{Q}}{\partial x} = 0\tag{5}
\end{align*}

ここで、$\boldsymbol{A}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}$と置きます。
$\boldsymbol{A}$を流束ヤコビ行列と呼びます。

$\boldsymbol{A}$の成分は後で具体的に計算するとして、結果だけ先に示しておきましょう。

どのようにして(6)を導くのかを今から見ていきます。
$\frac{\partial \boldsymbol{Q}}{\partial t} + \frac{\partial \boldsymbol{F}}{\partial t}=0$のままでは少しわかりにくいので、$\boldsymbol{Q}, \boldsymbol{F}$は3つの成分をもつ縦ベクトルなので、

\begin{align*}
\frac{\partial Q_i}{\partial t} + \frac{\partial F_i}{\partial x}= 0\tag{7}
\end{align*}

のように書くことにします。
そうすると、(7)は$\frac{\partial Q_i}{\partial t} + \sum_{j}\frac{\partial F_i}{\partial Q_j}\frac{\partial Q_j}{\partial x}= 0$となり、$A_{ij}=\frac{\partial F_i}{\partial Q_j}$と置いていることになります。
$j$が2回続いているので$\sum_{j}$は省略して$\frac{\partial Q_i}{\partial t} +\frac{\partial F_i}{\partial Q_j}\frac{\partial Q_j}{\partial x}= 0$と書くことが多いです。

わかりやすく$i$成分だけを書くと$j=1,2,3$だけを足すので、$\frac{\partial Q_i}{\partial t} +\frac{\partial F_i}{\partial Q_1}\frac{\partial Q_1}{\partial x}+\frac{\partial F_i}{\partial Q_2}\frac{\partial Q_2}{\partial x}+\frac{\partial F_i}{\partial Q_3}\frac{\partial Q_3}{\partial x}= 0$となります。

行列で書くと以下のようになっていることからも理解できることでしょう。

\begin{align*}
\frac{\partial}{\partial t}
\begin{bmatrix}
\rho u \\
p+\rho u^2 \\
(e+p)u \\
\end{bmatrix}+
\begin{pmatrix}
A_{11} & A_{12} & A_{13} \\
A_{21} & A_{22} & A_{23} \\
A_{31} & A_{32} & A_{33} \\
\end{pmatrix}
\frac{\partial}{\partial x}
\begin{bmatrix}
\rho \\
\rho u \\
e \\
\end{bmatrix}=0
\end{align*}

つまり、$\boldsymbol{A}$の成分は以下のようになっています。

\begin{align*}
\boldsymbol{A}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}=
\begin{pmatrix}
\frac{\partial F_1}{\partial Q_1}& \frac{\partial F_1}{\partial Q_2} & \frac{\partial F_1}{\partial Q_3} \\
\frac{\partial F_2}{\partial Q_1}& \frac{\partial F_2}{\partial Q_2} & \frac{\partial F_2}{\partial Q_3} \\
\frac{\partial F_3}{\partial Q_1}& \frac{\partial F_3}{\partial Q_2} & \frac{\partial F_3}{\partial Q_3} \\
\end{pmatrix}\\
=
\begin{pmatrix}
A_{11} & A_{12} & A_{13} \\
A_{21} & A_{22} & A_{23} \\
A_{31} & A_{32} & A_{33} \\
\end{pmatrix}
\end{align*}

例えば、$A_{23}=\frac{\partial F_1}{\partial Q_2}=\frac{\partial (e+p)u}{\partial (\rho u)}$と調子で計算していくことになります。

理想気体の状態方程式

理想気体の状態方程式$pV=nRT$を少し変形しておきましょう。
$pV=nRT$を$p=\frac{m}{V}\frac{R}{M}T$と置くと、$p=\rho \tilde{R}T$となります。
ここで、$M$は気体のモル質量、$\rho=\frac{m}{V}$は密度、$\tilde{R}=\frac{R}{M}$は比気体定数です。

さらに、$\tilde{R}=\tilde{C_p}-\tilde{C_v}=\tilde{C_v}(\gamma-1)$と置けば、理想気体の状態方程式は$p=\rho \tilde{C_v}(\gamma-1)T$となります。

$E$は単位体積あたりの内部エネルギー$E=_rho \tilde{C_v}T$とも書くことができ、単位体積あたりの全エネルギー$e=E+\frac{1}{2}\rho u^2$であるので、

\begin{align*}
p=(\gamma-1)(e-\frac{1}{2}\rho u^2)\tag{8}
\end{align*}

となります。

単位体積あたりの全エネルギー$e= E+\frac{1}{2}\rho u^2$
エンタルピー$H=\rho h$
単位体積当たりの全エンタルピー$h=e+p$

わざわざ気体定数$R$に対して、比気体定数$\tilde{R}=\frac{R}{M}$や$\tilde{C_p}, \tilde{C_v}$などチルダ（~）を頭につけたかというと、きっちり区別しておいた方がわかりやすいからです。
まず、気体定数$R$や定圧比熱$C_p$、定積比熱$C_v$の単位は[J/K mol]ですが、$\tilde{R}=\frac{R}{M}$や$\tilde{C_p}, \tilde{C_v}$の単位は[J/K kg]です。どちらの単位が好まれて使われているかは、分野によりますが、単位が違うと値が異なります。

流束ヤコビ行列の成分の導出

さて、ようやく$\boldsymbol{A}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}=
\begin{pmatrix}
A_{11} & A_{12} & A_{13} \\
A_{21} & A_{22} & A_{23} \\
A_{31} & A_{32} & A_{33} \\
\end{pmatrix}$の計算に移ります。

\begin{align*}
\frac{\partial f_1}{\partial q_1} = \frac{\partial \rho u}{\partial q_1} = \frac{\partial q_2}{\partial q_1} = 0
\end{align*}

\begin{align*}
\frac{\partial f_2}{\partial q_2} &= \frac{\partial (\rho u^2 + p)}{\partial m_2}\\
&= \frac{\partial \left( \frac{(q_2)^2}{\rho} + (\gamma – 1) \left[ e – \frac{(q_2)^2}{2 \rho} \right] \right)}{\partial q_2}\\
&= \frac{\partial}{\partial q_2} \left( \frac{(q_2)^2}{\rho} – (\gamma – 1) \frac{(q_2)^2}{2 \rho} \right) \\
&= \frac{2 q_2}{\rho} – (\gamma – 1) \frac{q_2}{\rho}\\
&= (3 – \gamma) u
\end{align*}

\begin{align*}
\frac{\partial f_3}{\partial m_3} &= \frac{\partial ((e + p)u)}{\partial m_3}\\
&= \frac{\partial \left( \left( e + (\gamma – 1) \left[ e – \frac{(m_2)^2}{2 \rho} \right] \right) u \right)}{\partial m_3}\\
&= \frac{\partial}{\partial m_3} \left( \left( e + (\gamma – 1) \left[ e – \frac{(m_2)^2}{2 \rho} \right] \right) u \right)\\
&= \gamma u
\end{align*}

※計算途中は読者に任せるとして結果はこちらです。

流束ヤコビ行列

まとめ

本記事では、オイラー方程式に基づき、流束ヤコビ行列の導出を解説しました。以下に主要なポイントをまとめます。

1. システム方程式の基本構造

圧縮性流体の流れをモデル化する際に用いられる基礎方程式であり、質量保存則、運動量保存則、エネルギー保存則の統合形として表現されます。
システム方程式：$\frac{\partial \boldsymbol{Q}}{\partial t} + \frac{\partial \boldsymbol{F}}{\partial x} = 0$
状態変数ベクトル$\boldsymbol{Q}$と流束ベクトル$\boldsymbol{F}$を具体化することで、圧縮性流れの各保存則に対応。

2. 流束ヤコビ行列の定義

流束$\boldsymbol{F}$を状態変数 $\boldsymbol{Q}$で偏微分した行列：$\boldsymbol{A}=\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{Q}}=\begin{pmatrix}0 & 1 & 0 \\
-\frac{3-\gamma}{2}u^2 & (3-\gamma)u & \gamma -1 \\
\big(\frac{\gamma-1}{2}-h\big)u & h-\big(\gamma-1\big)u^2 & \gamma u
\end{pmatrix}$
流束ヤコビ行列$\boldsymbol{A}$は、オイラー方程式をより簡単に扱える形へ変換し、流れの特性を表すために用いられます