非定常流に対する有限体積法

流体力学において、非定常流とは、その性質が時間依存する流体の輸送現象を指します。この種の現象は連続体力学の枠組みに属し、支配方程式には微小差分、すなわち連続的な空間変数と時間変数に関する関数の微分が含まれます。

有限体積法は、シミュレーション対象の形状を離散化することでシステムを数値的に解くことができる計算流体力学の手法であり、これにより、解の段階的な計算において流体特性の時間微分も考慮に入れることができる。^{[ 1 ]}

非定常流におけるスカラーの輸送に関する保存方程式は、一般的に次の形式をとる^{[ 2 ]。}

${\displaystyle {\frac {\partial \rho \phi }{\partial t}}+\operatorname {div} \left(\rho \phi \upsilon \right)=\operatorname {div} \left(\Gamma \operatorname {grad} \phi \right)+S_{\phi },}$

どこ${\displaystyle \rho }$密度であり、${\displaystyle \phi }$これは、あらゆる流体の流れの保存的な形態である。${\displaystyle \Gamma }$は拡散係数であり、${\displaystyle S}$はソース項です。${\displaystyle \phi }$流体要素（対流）から外れると${\displaystyle \operatorname {div} \left(\rho \phi \upsilon \right),}$そして${\displaystyle \operatorname {div} \left(\Gamma \operatorname {grad} \phi \right)}$増加率は${\displaystyle \phi }$拡散によるもの。$${\displaystyle S_{\phi }}$$増加率は${\displaystyle \phi }$情報源と${\displaystyle {\frac {\partial \rho \phi }{\partial t}}}$増加率は${\displaystyle \phi }$流体要素（過渡状態）の。

方程式の第1項は流れの非定常性を反映しており、定常流の場合には存在しない。支配方程式の有限体積積分は、制御体積内および有限時間ステップΔtにわたって実行される。

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left({\frac {\partial \rho \phi }{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n\cdot {\rho \phi u}\,\mathrm {d} A\right)\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n\cdot \left(\Gamma \operatorname {grad} \phi \right)\,\mathrm {d} A\right)\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S_{\phi }\,\mathrm {d} V\,\mathrm {d} t}$

方程式の定常部分の制御体積積分は、定常状態支配方程式の積分と同様です。ここでは、方程式の非定常成分の積分に焦点を当てる必要があります。積分手法を理解するために、一次元非定常熱伝導方程式を参照します。^{[ 3 ]}

${\displaystyle \rho c{\frac {\partial T}{\partial t}}={\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}+S}$

${\displaystyle \int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} V\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}{\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}\,\mathrm {d} V\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S\,\mathrm {d} V\,\mathrm {d} t}$

${\displaystyle \int _{e}^{w}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\int _{t}^{t+\Delta t}\left[\left(kA{\frac {\partial T}{\partial x}}\right)_{e}-\left(kA{\frac {\partial T}{\partial x}}\right)_{w}\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

ここで、節点の温度が制御体積全体に共通しているという仮定を置くと、方程式の左辺は次のように書ける^{[ 4 ]。}

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\rho c\left(T_{P}-{T_{P}}^{O}\right)\Delta V}$

1次後退差分法を用いると、等式の右辺は次のように表せる。

${\displaystyle \rho c\left(T_{P}-{T_{P}}^{0}\right)\Delta V=\int _{t}^{t+\Delta t}\left[\left(K_{e}A{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}A{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

方程式の右辺を評価するために、重み付けパラメータを使用します。${\displaystyle \theta }$0から1の間であり、積分を次のように書きます。${\displaystyle T_{P}}$

${\displaystyle I_{T}=\int _{t}^{t+\Delta t}T_{P}\,\mathrm {d} t=\left[\theta T_{P}+\left(1-\theta \right){T_{P}}^{0}\right]\Delta t}$

さて、最終的な離散化方程式の正確な形式は、の値に依存します。${\displaystyle \Theta }$分散が${\displaystyle \Theta }$は${\displaystyle 0<\Theta <1}$計算に使用されるスキーム${\displaystyle T_{P}}$値によって${\displaystyle \Theta .}$したがって

${\displaystyle \rho c{\frac {\left(T_{P}-{T_{P}}^{0}\right)}{\Delta t}}\Delta x=\theta \left[\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]+(1-\theta )\left[\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]+{\bar {S}}\Delta x}$

• 明示的スキーム：明示的スキームでは、ソース項は次のように線形化されます。${\displaystyle b=S_{u}+{S_{P}}{T_{P}}^{0}}$. 代替品${\displaystyle \theta =0}$明示的な離散化を得るには、次のようになります。^{[ 5 ]}

${\displaystyle a_{P}T_{P}=a_{w}{T_{w}}^{0}+a_{e}{T_{e}}^{0}+\left[{a_{P}}^{0}-\left(a_{w}+a_{e}-S_{P}\right)\right]{T_{P}}^{0}+S_{u}}$

どこ${\displaystyle a_{P}={a_{P}}^{0}}$注目すべき点の1つは、右辺には古いタイムステップの値が含まれているため、左辺は時間方向の順方向マッチングによって計算できるということです。このスキームは後方差分に基づいており、テイラー級数の打ち切り誤差は時間に関して1次です。すべての係数は正である必要があります。一定のkと均一なグリッド間隔の場合、${\displaystyle \delta x_{PE}=\delta x_{WP}=\Delta x}$この条件は次のように記述できます。

${\displaystyle \rho c{\frac {\Delta x}{\Delta t}}>{\frac {2K}{\Delta x}}}$

この不等式は、使用できる最大時間ステップに厳しい条件を課し、このスキームに対する重大な制約となります。空間精度を向上させるには、最大時間ステップを2乗に比例して小さくする必要があるため、非常にコストがかかります。${\displaystyle \Delta x}$^{[ 6 ]}

• クランク・ニコルソン方式 ：クランク・ニコルソン法は、以下の設定によって得られる。${\displaystyle \theta ={\frac {1}{2}}}$離散化された非定常熱伝導方程式は次のようになる。

${\displaystyle a_{P}T_{P}=a_{E}\left[{\frac {T_{E}+{T_{E}}^{0}}{2}}\right]+a_{W}\left[{\frac {T_{W}+{T_{W}}^{0}}{2}}\right]+\left[{a_{P}}^{0}-{\frac {a_{E}}{2}}-{\frac {a_{W}}{2}}\right]{T_{P}}^{0}+b,}$

どこ${\displaystyle a_{P}={\frac {a_{W}+a_{E}}{2}}+{a_{P}}^{0}-{\frac {S_{P}}{2}}.}$

新しい時間レベルでのTの未知値が複数あるため、この方法は陰解法であり、各時間ステップで全てのノード点の連立方程式を解く必要があります。${\displaystyle {\frac {1}{2}}<\theta <1}$Crank-Nicolson スキームを含め、すべての時間ステップの値に対して無条件に安定であるため、物理的に現実的で有界な結果を得るためには、すべての係数が正であることを保証することがより重要です。これは、係数が${\displaystyle {T_{P}}^{0}}$以下の条件を満たす

${\displaystyle {a_{P}}^{0}=\left[{\frac {a_{E}+a_{W}}{2}}\right],}$

これは

${\displaystyle \Delta t<\rho c{\frac {\Delta x^{2}}{K}},}$

したがって、クランク・ニコルソン法は中心差分法に基づいており、時間方向の精度は2次である。計算全体の精度は空間差分法にも依存するため、クランク・ニコルソン法は通常、空間中心差分法と組み合わせて使用​​される。

3.完全陰解法：Ѳの値を1に設定すると、完全陰解法が得られます。離散化された方程式は次のとおりです。 ^{[ 7 ]}

${\displaystyle a_{P}T_{P}=a_{W}T_{W}+a_{E}T_{E}+{a_{P}}^{0}{T_{P}}^{0}+S_{u}}$

${\displaystyle a_{P}={a_{P}}^{0}+a_{W}+a_{E}-S_{P}}$

方程式の両辺には新しい時間ステップにおける温度が含まれており、各時間レベルで代数方程式系を解く必要があります。時間進行手順は、与えられた初期温度場から始まります。 ${\displaystyle T^{0}}$時間ステップを選択した後、連立方程式を解く。${\displaystyle \Delta t}$次に解決策${\displaystyle T}$割り当てられています${\displaystyle T^{0}}$そして、この手順を繰り返して解をさらに時間ステップ進めます。すべての係数が正であることがわかります。これにより、陰解法は任意の時間ステップサイズに対して無条件に安定します。このスキームの精度は時間に関して1次であるため、結果の精度を確保するには小さな時間ステップが必要です。陰解法は、その堅牢性と無条件の安定性から、汎用的な過渡計算に推奨されます。