時間項の離散化導入-陰解法 超でかい行列作った

非圧縮性流体

空間の離散化はやったが、時間項の離散化はまだだったので、解説していきます。

簡単の為に、1次元熱拡散方程式で考えていきます。1次元は下のようになるのでした。

$$
\begin{align}
&\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \tag{1} \label{eq:scale-factor-1}
\end{align}
$$

空間微分の方は離散化すればいいです。

$$
\begin{align}
&\alpha \frac{\partial^2 T}{\partial x^2} =\alpha \frac{T_{i+1} – 2T_i + T_{i+1}} {\Delta x^2}\tag{2} \label{eq:scale-factor-2}
\end{align}
$$

時間項はどうしたらいいのでしょうか?とりあえず\eqref{eq:scale-factor-1}の右辺をFと置いて積分を試してみます。

$$
\begin{align}
&\int_t^{t+\Delta t}\frac{\partial T}{\partial t} dt = \int_t^{t+\Delta t} F dt\tag{3} \label{eq:scale-factor-3}
\end{align}
$$

左辺はいいとして、右辺は無理ゲーです。Fは時間的に変化するに決まってるので、時間の関数にもなっていることは確実ですが、その形が全然わからないので、積分しようがありません。そこで、時間に対しても離散化を試みます。

方法としては、単純に考えると上の3通りの方法が考えられます。
空間ではi点での値を\(T_i\)と表したように、n時点での値を\(T^n\)と表します。n乘ではないので注意してください。

一番左の近似を考えると、\eqref{eq:scale-factor-3}式は以下のように近似できることになります。

$$
\begin{align}
&T(t+\Delta t) -T(t) = \int_t^{t+\Delta t} F dt\\\\
&T^{n+1} – T^n = \Delta t \cdot F^n \\\\
&T^{n+1} = T^n + \Delta t \cdot F^n
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$

このようなn時点での値を使ってn+1時点の値を求めるような、近似の仕方を「陽解法」と言います。陽解法では、前の値を代入してやればそのまま解けるシンプルな方法ですが、どんどん誤差が大きくなってしまったり、適切に時間幅等を決めないと発散してしまったりします。
次に真ん中の考え方でやると以下のようにも近似することができます。

$$
\begin{align}
&T^{n+1} – T^n = \Delta t \cdot F^{n+1} \\\\
&T^{n+1} – \Delta t \cdot F^{n+1} = T^n
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$

この方法を「陰解法」と言います。この式を満たすように、連立方程式を解いて値を決めてやると言うイメージです。発散の心配がなく、こっちの方法を採用することがほとんどです。
もう一つ、右の考え方でやると以下のようになります。

$$
\begin{align}
&T^{n+1} – T^n = \Delta t \frac{F^n + F^{n+1}}{2} \\\\
&T^{n+1} -\Delta t \frac{F^{n+1}}{2}= T^n + \Delta t \frac{F^n}{2}
\tag{6} \label{eq:scale-factor-6}
\end{align}
$$

この近似方法は陰解法に分類されるのですが、この場合を「クランク・ニコルソン法」と言います。


陰解法の行列の形

陰解法の解き方の前にまず、どのようなものを解くのかを把握していきます。ひとまず、\eqref{eq:scale-factor-1}の1次元熱拡散で考えます。
格子番号、格子数をそれぞれi、Iとして陰解法で離散化していくと、次のようになります。

 

$$
\begin{align}
&\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \\\\
&\frac{ T_i^{n+1} – T_i^n }{\Delta t} = \alpha \frac{ T_{ i-1 }^{n+1} – 2T_{i}^{n+1} + T_{i+1}^{n+1} }{\Delta x^2} \\\\
&-r T_{i-1}^{n+1} +(1+2 r) T_{i}^{n+1} -r T_{i+1}^{n+1}=T_i^n \quad (r=\alpha \Delta t/\Delta x^2 ) \tag{7}\label{eq:scale-factor-7} \\\\
&(i=1,2,…,I)
\end{align}
$$

 

この式にi=1,2…,Iまで入れていくと、I個の連立方程式ができるのが分かるでしょうか?それを行列にしてやります。i=1,Iは境界条件にあたり表現がめんどくさいので省いて行列でかくと、下のようになります。本来は、I=1,Iも入れてやって、そこに当たる行列の成分には境界条件から導くものが入ります。

$$\begin{align}
& \left[ \begin{array}{cccccc} 1+2r & -r &  &  &  & 0 \\ -r & 1+2r & -r &  &  & \\  & \ddots & \ddots & \ddots &  &  \\  &  & -r & 1+2r & -r & \\ 0 &  &  &  & -r & 1+2r \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} T_{2}^{n+1} \\ T_{3}^{n+1} \\ \vdots \\ T_{I-2}^{n+1} \\ T_{I-1}^{n+1} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} T_{2}^{n} \\ T_{3}^{n} \\ \vdots \\ T_{I-2}^{n} \\ T_{I-1}^{n} \end{array} \right]
\end{align}$$

と、このようになります。この連立方程式を解いてやればいいです。まあこれはこれで難しいのですが。とりあえず、もし逆行列が求まるなりすれば、これは解けます。1次元はこれでできました。
1次元ならできることがわかりますが、3次元になったらどのようにすればいいのでしょうか?以外と1次元で分かった気になって、3次元の場合が分からなくなったりします。3次元だから、3つ行列を解くのかな?みたいな勘違いをしたり。3次元の場合も同じだが、ちゃんと1つの連立方程式にします。でかいけど、行列を求めてみます。
まず、3次元の場合の離散化式はx,y,z方向の格子番号をそれぞれi,j,kとして次のようになります。また、座標が\(x_1, x_2, x_3\)だと書きにくいので、x,y,zを使うことにします。

 

$$
\begin{align}
&\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x_j^2} \\\\
&\frac{\partial T}{\partial t} = \alpha \left[ \frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right] \\\\
&\frac{ T_{i,j,k}^{n+1} – T_{i,j,k}^{n} }{\Delta t} = \alpha \left[ \frac{ T_{ i-1,j,k }^{n+1} – 2T_{i,j,k}^{n+1} + T_{i+1,j,k}^{n+1} }{\Delta x^2} + \frac{ T_{ i,j-1,k }^{n+1} – 2T_{i,j,k}^{n+1} + T_{i,j+1,k}^{n+1} }{\Delta y^2} + \frac{ T_{ i,j,k-1 }^{n+1} – 2T_{i,j,k}^{n+1} + T_{i,j,k+1}^{n+1} }{\Delta z^2} \right]\\\\
& -X T_{i-1,j,k}^{n+1} -Y T_{i,j-1,k}^{n+1}-Z T_{i,j,k-1}^{n+1} +(1+2X+2Y+2Z) T_{i,j,k}^{n+1} -X T_{i+1,j,k}^{n+1} -Y T_{i,j+1,k}^{n+1} -Z T_{i,j,k+1}^{n+1}=T_{i,j,k}^{n}\tag{8}\label{eq:scale-factor-8} \\\\
&(X = \alpha \Delta t/\Delta x^2, Y = \alpha \Delta t/\Delta y^2, Z = \alpha \Delta t/\Delta z^2  ) (i=1,2,…,I, j=1,2,…,J, k=1,2,…,K)
\end{align}
$$

 

メッシュに振られた番号としては、このようになっているはずです。

ここで、連立方程式にするときにどう考えるかというと、x,y,zを順番に番号をふって、\(I \times J \times K\)の行列を作ってやるしかないです。

つまり、上の図のような順番で行列に入れてやります。最初はi,j,k=1,1,1から初めて、次はi,j,k=2,1,1、i,j,k=I,1,1まで行ったら、i,j,k=1,2,1、と繰り返していき、i,j,k = I,J,1まで行ったら、i,j,k=1,1,2としてと繰り返して、超巨大な行列にするのです。

でかすぎるが、行列にするとこんな感じです。横スクロールできるはずなので見てください。かなりイメージがつかめると思います。ただし、1+2X+2Y+2Z=Wと置いてあります。

$$\begin{align}
& \left[ \begin{array}{ccccccccccccccccccccccccccccccc} W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  & \dots &  &  &  &  &  &  &   &  &  &  &  & 0 \\ -X & W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  & \dots &  &  &  &  &  &  &  &  &  &  &  \\ & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  &  &  & \dots &  &  &  &  &  &  &  &  &  & \\ &  & -X & W &  & & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  & \dots &  &  &  &  &  &  & \\ -Y &  & \dots &  & W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  & \dots &  &  &  &   &  &  &  &  \\ & -Y &  & \dots & -X & W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  & \dots &  &  &  &  &  &  & \\ &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  &  &  & \dots &  &  &  &  &  &  \\ &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  &  &  &  &  &  &  & \\  &  &  &  & -Y &  & \dots &  & W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  &  &  &  &  &  &  \\ -Z &  &  & \dots &  & -Y &  & \dots & -X & W & -X &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  &  &  &  &  &  \\ & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  &  &  &  &  \\ &  & -Z &  &  & \dots &  & -Y &  & \dots & -X & W &  &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  &  &  &  &  &  &  \\ &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  &  & \\ &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  & \\ &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  &  &  \\ &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  & \\  &  &  &  &  &  &  & -Z &  &  & \dots &  & -Y &  & \dots & -X & W &  &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  &  &  & \\  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots &  &  \\  &  &  &  &  &  &  &  &  & -Z &  &  & \dots &  & -Y &  & \dots & -X & W &  &  & \dots &  & -Y &  &  & \dots &  &  & -Z &  \\ &  &  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \dots &  &  & \ddots \\  &  &  &  &  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  &  & \\  &  &  &  &  &  &  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots &  & \ddots &  \\ &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots  &  & \dots & \\ &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  & -Z &  &  & \dots &  & -Y &  & \dots & W & -X & \dots & \\ &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  & \ddots &  &  &  &  & \ddots  &  &  & \ddots & \ddots & \ddots \\ &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  & -Z &  &  & \dots &  & -Y & \dots & -X & W & -X \\ 0 &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  &  & -Z &  &  & \dots &  & -Y & \dots & -X & W \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} T_{2,2,2}^{n+1} \\ \vdots \\ T_{I-1,2,2}^{n+1} \\ T_{2,3,2}^{n+1} \\ \vdots \\ T_{I-1,3,2}^{n+1} \\ \vdots \\ \vdots \\ T_{I-1,j-1,2}^{n+1}  \\ T_{2,2,3}^{n+1} \\ \vdots \\ \vdots \\ T_{I-1,2,3}^{n+1} \\ T_{2,3,3}^{n+1} \\ \vdots \\ T_{I-1,3,3}^{n+1} \\ \vdots \\ \vdots \\ T_{I-1,j-1,3}^{n+1}\\ \vdots \\ \vdots \\ \vdots  \\ \vdots \\ \vdots \\ T_{I-1,J-1,K-1}^{n+1} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c}T_{2,2,2}^{n} \\ \vdots \\ T_{I-1,2,2}^{n} \\ T_{2,3,2}^{n} \\ \vdots \\ T_{I-1,3,2}^{n} \\ \vdots \\ \vdots \\ T_{I-1,j-1,2}^{n}  \\ T_{2,2,3}^{n} \\ \vdots \\ \vdots \\ T_{I-1,2,3}^{n} \\ T_{2,3,3}^{n} \\ \vdots \\ T_{I-1,3,3}^{n} \\ \vdots \\ \vdots \\ T_{I-1,j-1,3}^{n}\\ \vdots \\ \vdots \\ \vdots  \\ \vdots \\ \vdots \\ T_{I-1,J-1,K-1}^{n} \end{array} \right]
\end{align}$$
 

この行列は、(I-1)(J-1)(K-1)行、(I-1)(J-1)(K-1)列の行列です。境界条件もぶっこむ場合は、I*J*K行I*J*K列の行列を作ることになります。また、WとYは、I-1分離れていて、WとZは、(I-1)(J-1)分離れています。このイメージを頭に入れておくといいと思います。


陰解法の解き方

先ほど、行列を出しましたが、陰解法で時間ステップを進めるとは、先ほどの逆行列を解くという行為に他なりません。しかし、コンピュータを使ったとしても、あんな行列の逆行列を求められるでしょうか?まあ一応、できます。直接法と呼ばれ、完全LU分解とかが代表例ですかね。しかし、非常に計算時間がかかります。

そこで、反復法というものを使用します。反復法の思想自体は非常に簡単です。\eqref{eq:scale-factor-5}式で考えてみます。簡単な考え方としては、こうです。まず、\(T_{n+1}\)に\(T_{n}\)の値を入れてとりあえず計算してみます。そうして出てきた\(T_{n+1}\)を\({}^{2}T^{n+1}\)として、再び代入します。左上の添字は反復回数です。

 

$$
\begin{align}
T^{n+1} &= T^n + \Delta t \cdot F^{n+1} \\\\
T^{n+1} &\longleftarrow {}^{1}T^{n+1} = T^{n} \\
&\searrow \\
T^{n+1} &\longleftarrow {}^{2}T^{n+1}  \\\\
\tag{9} \label{eq:scale-factor-9}
\end{align}
$$

 

これを繰り返すうちに、\({}^{l+1}T_{n+1} \simeq {}^{l}T_{n+1} \)となったら、この方程式を解けたと言えないでしょうか?

反復法の原理自体はこんな感じです。もちろん、これと同じ作業を行列でやります。ただ、こんな単純なプログラムでは多分収束するがすごく遅いです。例えば、\(A \mathbb{ x } = \mathbb{ b } \)を解くときに、\(\mathbb{ r } = \mathbb{ b } – A \mathbb{ x }\)として、これを掛けて収束を早くしたりします。このアイデアが使われているアルゴリズムは、共役勾配法と呼ばれます。もうちょっとシンプルな方法は、ヤコビの反復法や、ガウス-ザイデル法とかがあります。

もっとも、収束を早くするには不完全LU分解等の前処理が必須で、前処理付き共役勾配法をPCGと言って流体計算だけ出なくいろんな行列計算に広く使われるものとなっています。プログラムを実装する時に、ちゃんとやりますが、東京大学さんの資料が割とわかりやすいので読んでみるといいと思います。

http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdf
http://nkl.cc.u-tokyo.ac.jp/13n/SolverDirect.pdf

 



[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]