工学系大学院生のブログ

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

第7-4回 二次元拡散方程式(実装までの式変形)[python, julia, fortran]


第7-4回では二次元拡散方程式を実装できる形に変形していきます。

$$\rho C \frac{\partial T}{\partial 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)$$


時間一次精度、陰解法で離散化します。第7-2回で行った方法を用います。


左辺は今までと同じです。

$$\rho C \frac{\partial T}{\partial t} \simeq \rho C \frac{\Delta T^n_{i,j}}{\Delta t}$$

$$\Delta T^n_{i,j}= T^{n+1}_{i,j}- T^{n}_{i,j}$$


右辺第一項は次のように書きます。詳しくは第7-2回をご覧ください。

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

$$ \simeq \frac{ \lambda_{i+1/2,j} }{\Delta x^2} (\Delta T_{i+1,j}^n +T_{i+1,j}^n-\Delta T_{i,j}^n -T_{i,j}^n )$$

$$-\frac{ \lambda_{i-1/2,j} }{\Delta x^2} (\Delta T_{i,j}^n +T_{i,j}^n -\Delta T_{i-1,j}^n -T_{i-1,j}^n) $$


右辺第二項も同様に行います。

$$\frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y} \right) $$

$$\simeq \frac{ \lambda_{i,j+1/2} }{\Delta y^2} (\Delta T_{i,j+1}^n +T_{i,j+1}^n-\Delta T_{i,j}^n -T_{i,j}^n )$$

$$- \frac{ \lambda_{i,j-1/2} }{\Delta y^2} (\Delta T_{i,j}^n +T_{i,j}^n -\Delta T_{i,j-1}^n -T_{i,j-1}^n) $$


まとめると次のようになります。

$$\rho C \frac{\Delta T^n_{i,j}}{\Delta t}=\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} (\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))+R_{i,j}^n$$



$$ R_{i,j}^n =\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 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))$$




さらに、反復数kを用いてガウスザイデル法の形に変形します。詳しくは第7-3回を参照してください。

$$ \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2,j} }{\Delta x^2} + \frac{ \lambda_{i-1/2,j} }{\Delta x^2}+ \frac{ \lambda_{i,j+1/2} }{\Delta y^2} + \frac{ \lambda_{i,j-1/2} }{\Delta x^2} \right ) \Delta T^{k+1}_{i,j} =\frac{1}{\Delta x^2} \left[\lambda_{i+1/2,j} \Delta T_{i+1,j}^k -\lambda_{i-1/2,j} \Delta T_{i-1,j}^k \right]$$ $$+ \frac{1}{\Delta y^2}\left[\lambda_{i,j+1/2} \Delta T_{i,j+1}^k -\lambda_{i,j-1/2} \Delta T_{i,j-1}^k \right]+R_{i,j}^n$$

$$ \Delta T^{k+1}_{i,j} =\left ( \frac{1}{\Delta x^2} \left[\lambda_{i+1/2,j} \Delta T_{i+1,j}^k -\lambda_{i-1/2,j} \Delta T_{i-1,j}^k \right]+ \frac{1}{\Delta y^2} \left [ \lambda_{i,j+1/2} \Delta T_{i,j+1}^k -\lambda_{i,j-1/2} \Delta T_{i,j-1}^k \right]+R_{i,j}^n \right ) $$ $$ / \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2,j} }{\Delta x^2} + \frac{ \lambda_{i-1/2,j} }{\Delta x^2}+ \frac{ \lambda_{i,j+1/2} }{\Delta y^2} + \frac{ \lambda_{i,j-1/2} }{\Delta x^2} \right ) $$




これで式としては実装できる形となりました。(ガウスザイデル法の場合,格子の並び順によっては一部kではなくk+1と書きますが,理解しやすい形のためこのまま行きます。)




しかし、反復法の収束性を見極める基準をまだ設けていません。



今回は2-ノルムが\(10^{-4}\)以下となったら計算を終了させたいと思います。


定常反復法は一定割合で収束していきます。そこで、真の解との差がほとんどなくなったら収束したとみなします。その基準を2-ノルムが\(10^{-4}\)以下としました。



下記のように行列があったとします。\(\bf A\)は\(i × i\)行列, \(\bf x\)と\(\bf b\)は\(i\)の列行列です。

$$\bf A \bf x =\bf b$$


反復数\(k\)のときを考えると誤差があるため、以下の式の右辺は0となりません。

$$\bf A \bf x^{k} -\bf b=\bf e$$


そこで、2-ノルムをとります。

$$2-norm=\sqrt{\sum_{k=0}^i (a_{k,j}x_{k}-b_{k})^2}$$


いつでも使えるように\(\bf\)を用いて規格化します。

$$2-norm=\sqrt{\frac{\sum_{k=0}^i (a_{k,j}x_{k}-b_{k})^2}{\sum_{k=0}^i (b_{k})^2}}$$


以上の値を用いて

$$2-norm \le 10^{-4}$$


となったとき計算終了とします。



実装までの式変形がすべて終わったので、次回は実際にコードに落とし込みたいと思います。


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


By hide

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

コメントを残す

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

トラックバックURL: https://teru-hide.com/cfd-trans-2d-diffusion-equation/trackback/
にほんブログ村 英語ブログへ
にほんブログ村