Implicit time integration

陰的時間スキーム

Last Update:2009年05月28日

概要

質量行列と剛性行列からなる次のような系に対する陰的な時間積分スキームを $\b{C}$ が線形，非線形の場合について紹介する.

$\b{M}\b{a}+\b{C}\b{v}=\b{F}$

陰的な時間増分法を行う場合、時刻、時間ステップの速度 $\b{v}_n$ と加速度 $\b{a}_n$ から時刻 $t+\Delta t$ ,時間ステップ n+1 での速度 $\b{v}_{n+1}$ と加速度 $\b{a}_{n+1}$ を求める際には以下の時刻 n+1 における平衡方程式と時間積分関係式の２つの式を連立させて解く

時刻 $t+\Delta t$ ,時間ステップにおける平衡方程式

$\b{M}\b{a}_{n+1}+\b{C}(\b{v}_{n+1})\b{v}_{n+1}=\b{F}$

時間積分関係式

$\b{v}_{n+1} = \b{v}_n + \int^{t+\Delta t}_t \b{a} dt$

時間積分関係式を数値的には次のように時間方向に離散化する

時間積分関係式
$\b{v}_{n+1}=\b{v}_{n}+\{(1-\gamma)\b{a}_{n} + \gamma\b{a}_{n+1} \}\Delta t$

ここで $\gamma$ は 0.5 からの値をとるパラメータで、時間積分の精度、数値粘性などを制御することができる。

$\gamma$ を0.5とした場合はCrank-Nicolson法と等価となり二次精度を持つ。 $\gamma$ を1.0とした場合は完全陰解法、または後退オイラー法(Backward-Euler Method)と呼ばれ一次精度であるが安定性が高い

線形問題

熱拡散問題，非定常移流拡散問題を解く時は線形問題として取り扱うことができる．この時平衡方程式における係数行列 $\b{C}$ が定数なので，平衡方程式は次のようになる．

$\b{M}\b{a}_{n+1}+\b{C}\b{v}_{n+1} = \b{F}$

これに時間積分関係式を代入して変形することで，前ステップの解から次ステップの解を求めることができる．線形性から一度の行列計算で解を求めることができる．前ステップの解からの差を求める増分解析と，直接解を求める方法について方法を示す

直接解を求める方法

時間積分関係式を離散化された平衡方程式に代入すると次のようになる。

$(\b{M}+\gamma\Delta t\b{C})\b{a}_{n+1}=\b{F}-\b{C}\{\b{v}_n-\Delta t(1-\gamma)\b{a}_n\}$

これを解いて $\b{a}_{n+1}$ を求める。また、時間積分関係式から $\b{v}_{n+1}$ が求まる

増分解析

上の増分をとると、次のようになる。

$(\b{M}+\gamma\Delta t\b{C})\Delta\b{a}=\b{F}-\b{C}(\b{v}_n+\Delta t\b{a}_n)-\b{M}\b{a}_n$

$\b{a}_{n+1}=\b{a}_n+\Delta\b{a}$

$\Delta\b{v}=\Delta t\b{a}_n+\gamma\Delta t\Delta \b{a}$

$\b{v}_{n+1}=\b{v}_n+\Delta\b{v}$

非線形問題

特に速度 $\b{v}$ に対して係数行列 $\b{C}$ が非線形であることを考える．流体や，非線形熱伝導問題がこのような形を取る．

$\b{M}\b{a}_{n+1}+\b{C}(\b{v}_{n+1})\b{v}_{n+1} = \b{F}$

このような方程式は反復的な計算をさせて解へ収束させる方法を取る．代表的な方法であるPMA法について説明する．

ＰＭＡ法

ＰＭＡ法(Predictor Multicorrector Algorithm)は,一度次のステップの解について時間積分関係式を満たす範囲内で予測をして、それが時刻 $t+\Delta t$ での平衡方程式を満足するようにするために、Newton-Raphson法による増分解析を用いて、時間積分関係式を満たすように解を修正していく手法である。

時間積分関係式を満たす時刻 $t+\Delta t$ での流速、加速度の予測 $\b{v}_{n+1}^{pred}$ , $\b{a}_{n+1}^{pred}$ は無数に考えられるが代表的な予測は

加速度０の予測

$\{ \b{a}_{n+1}^{pred}=0\\ \b{v}_{n+1}^{pred}=\b{v}_{n}+(1-\gamma)\b{a}_{n}\Delta t$

加速度変化なしの予測

$\{ \b{a}_{n+1}^{pred}=\b{a}_{n}\\ \b{v}_{n+1}^{pred}=\b{v}_{n}+\b{a}_{n}\Delta t$

代入すると簡単にわかるが、これらの式は時間積分関係式を満たしている。

これらの予測をもちいたＰＭＡの手続きは以下のとおり

PMAのアルゴリズム

速度、加速度の初期化

$\{\b{a}_{n+1}^1=\b{a}_{n+1}^{pred}\\\b{v}_{n+1}^1=\b{v}_{n+1}^{pred}$

iを１にセット

増分量 $\Delta\b{a}^i$ の計算

$\b{M}^*\Delta\b{a}^i=\b{R}^i$

$\{\b{M}^*=\b{M}+\gamma\Delta t\b{C}(\b{v}_{n+1}^i)\\ \b{R}^i = \b{F}-\b{M}\b{a}^i_{n+1}-\b{C}(\b{v}^i_{n+1})\b{v}^i_{n+1}$

解の更新

$\{\b{v}_{n+1}^{i+1}=\b{v}_{n+1}^i+\gamma\Delta t\Delta\b{a}^i\\ \b{a}_{n+1}^{i+1}=\b{a}_{n+1}^i+\Delta\b{a}^i$

さらに反復が必要ならiを１進めて３番に進む

$\Delta t$

修正ＰＭＡ法

ＰＭＡ法では系の非線形性が強い場合不安定になり、収束しない場合がある。これは１回目の増分量の計算において非線形項を予測値を元に計算しているためである。そこで１回目の増分計算に限り、非線形項を前の時間ステップの速度 $\b{v}_n$ を元に計算する。つまり、一回目の反復で $\b{C}(\b{v}_{n+1}^i) \right \b{C}(\b{v}_n)$ と置き換える。