工学系大学院生のブログ

2人の大学院生による雑記ブログ

第7-2回 陰解法と有限体積法[python, julia, fortran]



有限体積法における陰解法を用いた離散化

第7-2回では有限体積法における陰解法を用いた離散化を行います。



例として,一次元の拡散方程式は次のように書けます。

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)$$

ここで\(T=T(t,x), \lambda=\lambda(x) \)とします。


セルjにおける離散化を行います。


左辺時間項の離散化は一次精度で通常通り行います。

$$\frac{\partial }{\partial t}(\rho C T)=\rho C \frac{\Delta T^n _i}{\Delta t}$$

ここで,\(\Delta T^n = T^{n+1}-T^n \)です。




次に右辺ですが,今までと違い現在の時間\(n\)ではなく次の時間\(n+1\)を用います。

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)= \frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) ^{n+1}$$



ここから詳しく変形を行います。まず下記の様に時間方向に対してテイラー展開を行います。

$$T^{n+1}=T^n+\Delta t \frac{\partial T}{\partial t} +O(\Delta t^2)$$

$$\simeq T^n+\Delta t \frac{\partial T}{\partial t}$$

$$\simeq T^n +\Delta t \frac{T^{n+1}-T^n}{\Delta t}$$

$$= T^n +\Delta T^n $$

先程と同様に,\(\Delta T^n = T^{n+1}-T^n \)です。



以上のことから

$$ \frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) ^{n+1}= \frac{\partial}{\partial x} \left( \lambda \frac{\partial}{\partial x}( T^n +\Delta T^n)\right) $$


少し変な形ですが,分かりやすいように段階を踏みながら書いていきます。

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right) =\frac{1}{\Delta x}(\lambda_{i+1/2}\frac{\partial T}{\partial x}_{x=i+1/2}-\lambda_{i-1/2}\frac{\partial T}{\partial x}_{x=i-1/2})$$

$$=\frac{1}{\Delta x} ( \lambda_{i+1/2} \frac{T_{i+1}-T_i}{\Delta x}-\lambda_{i-1/2} \frac{T_{i}-T_{i-1}}{\Delta x} )$$


よって

$$ \frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right) ^{n+1}=\frac{1}{\Delta x^2} (\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2} (T_{i}^n -T_{i-1}^n))$$ $$+ \frac{1}{\Delta x^2} ( \lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2} (\Delta T_{i}^n -\Delta T_{i-1}^n))$$




以上の結果をまとめます。

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)$$

$$⇔ \rho C \frac{\Delta T^n_i}{\Delta t} = \frac{1}{\Delta x^2}(\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2}(T_{i}^n -T_{i-1}^n))$$ $$+ \frac{1}{\Delta x^2}(\lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2}(\Delta T_{i}^n -\Delta T_{i-1}^n))$$


これで一次元の場合の離散化ができました。



しかし,まだ終わりではありません。


\(\Delta T^n= T^{n+1}- T^{n} \)は未来の値のため,これが未知の変数となり上記の式では3つあります。


未知の変数を左辺に持っていき,右辺を既知の値にまとめます。

$$⇔ \rho C \frac{\Delta T^n_i}{\Delta t}- \frac{1}{\Delta x^2}(\lambda_{i+1/2} (\Delta T_{i+1}^n-\Delta T_i^n )-\lambda_{i-1/2}(\Delta T_{i}^n -\Delta T_{i-1}^n))$$ $$=\frac{1}{\Delta x^2}(\lambda_{i+1/2} (T_{i+1}^n-T_i^n )-\lambda_{i-1/2}(T_{i}^n -T_{i-1}^n)) $$

$$⇔\frac{\lambda_{i-1/2}}{\Delta x^2}\Delta T_{i-1}^n+( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2})\Delta T_{i}^n+\frac{\lambda_{i+1/2}}{\Delta x^2}\Delta T_{i+1}^n=R_i^n$$

$$\left[ \begin{array}{ccc} \frac{\lambda_{i-1/2}}{\Delta x^2} & ( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2}) & \frac{\lambda_{i+1/2}}{\Delta x^2} \end{array}\right] \left[ \begin{array}{c} \Delta T_{i-1}^n \\ \Delta T_{i}^n \\ \Delta T_{i+1}^n  \end{array}\right] =R_i^n$$


上記の式はセルiにおける式です。つまりセル分だけ上記の式を作ることができます。


例えば,空間セル10個のときは未知の変数は10個です。(\(\Delta T ^n_i,i=1~10\))


そして作れる式はセル分なので10個です。


未知数10個,式10個なので,これらを連立すれば答えを求めることができます。


連立方程式を解くための行列の扱い



連立方程式を解くには行列を用います。


行列を用いて空間全体の式をまとめると次のようになります。

$$\left[\begin{array}{ccccc} \ddots& \ddots & \ddots && \\ & \frac{\lambda_{i-1/2}}{\Delta x^2} & ( \frac{\rho C}{\Delta t}-\frac{\lambda_{i-1/2}}{\Delta x^2}-\frac{\lambda_{i+1/2}}{\Delta x^2}) & \frac{\lambda_{i+1/2}}{\Delta x^2} & \\ && \ddots & \ddots & \ddots \\ \end{array}\right] \left[\begin{array}{c} \vdots \\ \Delta T_{i-1}^n \\ \Delta T_{i}^n \\ \Delta T_{i+1}^n  \\ \vdots \end{array}\right] =\left[\begin{array}{c} \vdots \\ \vdots \\ R_i^n \\ \vdots \\ \vdots \end{array}\right] $$

$$⇔\bf A \bf \Delta \bf T = \bf R$$

$$⇔\bf \Delta \bf T = \bf A^{-1} \bf R$$


この行列\(\bf A\)の一番上と一番下は境界条件のため少し変わりますが,このような形になります。


逆行列を求めることができれば,\(\Delta T^n\)がわかり,下記の計算で最終的な答えが出せます。

$$T^{n+1}=T^{n}+\Delta T^n$$


次回はこの逆行列の求め方を行いたいと思います。


補足. 二次元の場合



二次元の場合は座標も増え大変ですが同様にセル\((i,j)\)において,陰解法で離散化を行うと次のようになります。

$$\frac{\partial }{\partial t}(\rho C T)=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x}\right)+ \frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y}\right)$$

$$⇔ \rho C \frac{\Delta T^n_i }{\Delta t}=\frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (T_{i+1,j}^n-T_{i,j}^n )-\lambda_{i-1/2,j}(T_{i,j}^n -T_{i-1,j}^n))$$ $$+ \frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (\Delta T_{i+1,j}^n-\Delta T_{i,j}^n )-\lambda_{i-1/2,j}(\Delta T_{i,j}^n -\Delta T_{i-1,j}^n))$$ $$+\frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (T_{i,j+1}^n-T_{i,j}^n )-\lambda_{i,j-1/2}(T_{i,j}^n -T_{i,j-1}^n))$$ $$+ \frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (\Delta T_{i,j+1}^n-\Delta T_{i,j}^n )-\lambda_{i,j-1/2}(\Delta T_{i,j}^n -\Delta T_{i,j-1}^n))$$


今回は一般的に使える様に,熱伝導率\(\lambda\)を定数にしなかったため複雑になりました。



次回は反復法を用いた逆行列の求め方を行います。


*コードは第7-5回にあります。


by hide

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

トラックバックURL: https://teru-hide.com/cfd-implicit-discretization/trackback/
にほんブログ村 英語ブログへ
にほんブログ村