流体シミュレーション基礎

Transcription

1 中部 CAE 懇話会 : 流体伝熱基礎講座 第 4 回 : 流体計算の時間進行法と計算アルゴリズム 名古屋工業大学 大学院工学研究科 機能工学専攻教授森西洋平 講義内容 :. 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域. 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラ格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 HSMAC 法 (SOLA 法 フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

2 第 章時間進行法 目的 : 階常微分方程式に対する時間進行法を理解する. 本章では変数 ( に関する次の 階常微分方程式の初期値問題の時間進行法 ( 数値時間積分法 を考える. d d ( ( ( ( > (. ( (. つまり 時間刻み幅 毎の時間離散点上の値 ( を初期条件として式 (. を満足させながら求めたい. ( ( L を 式 (. 時間進行法を構築するため まず ( を ( と次式を得る. ( ( d d まわりにテイラー展開して式 (. を使用する ( d ( d ( d ( d! d L L! d ( ( ( d ( ( d ( ( d ( ( ( d L ( L! d! d (. 時間に関する離散点上の値を表現するため ( ( ( の表記を使用する. 以降も同様. 最も単純な時間進行法はテイラー展開の右辺第 項以降を打ち切って構成される陽的オイラー法である. あるいは ( 陽的オイラー法では における勾配 を計算する. と から ( ( Eler exlc 図. 陽的オイラー法による時間進行

3 式 (. と陽的オイラー法をテイラー展開を通して比べると 打ち切誤差の初項が ( 陽的オイラー法は 次精度しか持たないことがわかる. O となり 時間進行法の精度を向上させる方法として以下がしばしは採用される. 線形多段解法 : を計算する際に 時刻 および での情報のみならず さらに過去 の情報 ( 時刻 も利用して高次精度化を行う. 計算の出発点近傍で別の時間進行法を使用する必要があること および時間刻み幅の変更が容易でないことに注意が必要. 代表は アダムス バシュフォース (AB 公式 アダムス モルトン (AM 公式 後退差分 (BD 公式. 段解法 : を計算する際に 基本的に時刻 での情報 ( およびそれから計算される値 のみを利用して高次精度化を行うもの. 代表は ルンゲ クッタ (RK 法. 以下本章では 線形多段階法としてアダムス バシュフォース (AB 公式 アダムス モルトン (AM 公式 後退差分 (BD 公式 また 段階法としてルンゲ クッタ (RK 法を取り上げ紹介する. また これら時間進行法の数値安定性についても示す. 補足 : 階常微分方程式を扱う理由流体運動の支配方程式は一般に偏微分方程式である輸送方程式として記述されるが なぜ 階常微分方程式を扱うのかを説明する. 例として外力を伴う 次元移流拡散方程式を考える. U ν f x x 変数 ( x が独立変数 と x の双方に依存するのでこの方程式は偏微分方程式である. またより一般的に移流速度 U U ( x 拡散係数 ν ν ( x および外力 f f ( x も時間と空間の双方に依存するものとする. ここで 空間に等間隔格子 ( x x N を使用し 離散値を ( x ( のように表現し 空間 階微分を 次精度中心差分で近似する. d U ν f ( N d x x 空間離散化後の変数は ( なので偏微分が常微分に変わる. さらに変数を [ ] T L N とベクトル表記し 端点 N に対する境界条件も含めて上式を書き直すと次の離散システム方程式を得る. d d ( A( ( ( ( ( A( は差分の係数を要素とする ( N N の行列 ( は外力と境界条件を含むベクトルである. このように流体運動の支配方程式を空間離散化した後に変数をベクトル表記すると 空間離散化後の時間発展方程式は基本的に式 (. の形となる. 空間多次元の場合も同様である. なお 離散システム方程式の変数はベクトル 式 (. の変数はスカラー であるが 時間進行法の構成に違いは生じない.

4 線形多段階法 ( N 段の線形多段階法 (Ler ml-se mehods は一般に次式で記述される. ただし. N N β 特に β ( とする次の時間進行法はアダムス バシュフォース公式 (Adms -Bshforh forml と呼ばれる陽解法 ( β となる. ( ( β β L β アダムス バシュフォース公式 ( Adms Bshforh ( また 図. アダムス バシュフォース公式による時間進行 β ( 公式 (Adms -Molo forml と呼ばれる陰解法となる. ( ( β β β L とする次の時間進行法はアダムス モルトン β アダムス モルトン公式 ( Adms Molo ( 図. アダムス モルトン公式による時間進行 4

5 さらに β β ( dfferece forml と呼ばれる. 後退差分公式も陰解法である. とする次の時間進行法は後退差分公式 (Bckwrd L ( 後退差分公式 ( Bckwrd dfferece ( 図.4 後退差分公式による時間進行 以下でそれぞれに対する 4 次精度までの公式を示す. 5

6 6 アダムス バシュフォース公式時間 次精度から 4 次精度までのアダムス バシュフォース (AB 公式は以下のとおりである. ( 陽的オイラー法 (AB ( ( AB ( ( 5 6 AB ( ( AB4 次精度の公式 (AB は陽的オイラー法である. 特に精度が明示されていない場合 流体計算では 次精度の公式 (AB をアダムス バシュフォース法と称する場合が多い. 補足 :AB 公式の導出 AB 公式を ( ( L β β β とおき 時刻 まわりのテイラー展開 ( L 4 6 d d d d d d ( ( ( L 6 q d d q d d q d d q L q から AB 公式が指定した時間精度 ( 次精度ならば ( O で近似されるように係数 β β β を定める. 例として 次精度アダムス バシュフォース法 (AB を取り上げる. ( β β AB の係数 β β を定めるためテイラー表を利用する. ( 6 / β β β β β β d d d d ( O と ( O の列それぞれの係数の合計を とする つの条件から 係数 β β に対する連立方程式が得られる.

7 β β β これを解き β β を得る. これは AB 公式を与える. ( ( なお テイラー表の 列目 (O( の合計は AB 公式の打ち切り誤差の初項 を与える. d 5 d ε β O( 6 d d 補足 : 陽解法による計算陽解法は既知量の単純代入計算のみで時間進行が行えるので計算コードの作成が容易であり また陰解法に比べて時間ステップあたりの計算量が非常に少なくて済む. 例えば AB は実際には上式を変形した次の形で使用される. ( これは既知量 ( 時刻 以前の量 の単純代入計算のみで未知量 が得られることを意味し 時間ステップの進行に繰り返し計算を必要としない. このような時間進行法は陽解法と呼ばれる. 7

8 8 アダムス モルトン公式時間 次精度から 4 次精度までのアダムス モルトン (AM 公式は以下のとおりである. ( 陰的オイラー法 (AM ( ( クランク ニコルソン法 (AM ( ( 8 5 AM ( ( AM4 次精度の公式 (AM は陰的オイラー法であり 次精度の公式 (AM は通常クランク ニコルソン (Crk-NcolsoCN 法と呼ばれる.AM 公式は の値を得るために ( を必要とするが このような時間進行法は陰解法と呼ばれる. 補足 :AM 公式の導出 AM 公式を ( ( L β β β とおき 時刻 まわりのテイラー展開 ( L 4 6 d d d d d d L 6 d d d d d d ( ( ( L 6 q d d q d d q d d q L q から AM 公式が指定した時間精度 ( 次精度ならば ( O で近似されるように係数 β β β を定める. 例としてクランク ニコルソン法 (AM を取り上げる. ( β β AM の係数 β β を定めるためテイラー表を利用する. ( 6 / β β β β β β d d d d

9 O( と O( の列それぞれの係数の合計をとするつの条件から 係数 β β に対する連立方程式が得られる. β β β これを解き β β を得る. これはクランク ニコルソン法 (AM を与える. ( ( なお テイラー表の 列目 (O( の合計はクランク ニコルソン法の打ち切り誤差の初項 を与える. d d ε β O( 6 d d 補足 : 陰解法による計算例えばクランク ニコルソン法 (AM は実際には上式を変形した次の形で使用される. ( ( 右辺は既知量であるが左辺は未知量 を含む関数であり 上式を満足する を特定するために通常繰り返し計算が必要となる. 9

10 後退差分公式時間 次精度から 4 次精度までの後退差分 (BD 公式は以下のとおりである. ( 陰的オイラー法 (BD ( 4 BD ( BD ( BD4 次精度の公式 (BD は陰的オイラー法である. これら後退差分公式はギアー (Ger 法と呼ばれることもある.BD 公式も陰解法である. 補足 :BD 公式の導出 BD 公式を ( L とおき 時刻 まわりのテイラー展開 ( ( [ ] ( [ ] L 6 q d d q d d q d d q ( ( [ ] ( [ ] L 6 q d d q d d q L q から BD 公式が指定した時間精度 ( 次精度ならば ( O で近似されるように係数 を定める. 例として 次精度後退差分公式 (BD を取り上げる. ( BD の係数 を定めるためテイラー表を利用する. 6 8 / 6 / / d d d d

11 O ( O( および O( の列それぞれの係数の合計をとするつの条件から 係数 に対する連立方程式が得られる. これを解き を得る. 4 これは BD 公式を与える. ( 4 なお テイラー表の4 列目 (O( の合計は BD 公式の打ち切り誤差の初項 を与える. ( 8 d d ε O( 6 d d

12 ルンゲ クッタ法式 (. を時間進行するための s 段ルンゲ クッタ (Rge-K 法 (RK 法 の一般公式は次式で与えられる. s ( c s c ただし 係数に対し s c ( s の条件を与える. ( ( ( s ( ( (.. (.. 係数 c および は要求精度を満たすように式 (. のテイラー展開との比較から決定される.RK 法の係数は次のブッチャー配列 (Bcher le で表現すると便利である. c c c M s M s s M L L O L L s s M ss s (.. なお に対して とするものを ( 陽的 RK 法と呼ぶ.s 段 ( 陽的 RK 法の到達可 は次表となることが知られている. 能精度 ( s s s ( s ( s s 8 9 精度と計算量の兼ね合いからは 4 段 4 次精度公式が最も効率が良いことがわかる. また > に対して とするものを半陰的 RK 法と呼ぶ.s 段半陰的 RK 法の到達可能精 s s となることが知られている. さらに ある > について のものを陰的 s である. 度は ( ( RK 法と呼ぶが s 段陰的 RK 法の到達可能精度は ( s 補足 : 係数に対する条件式 (.. に対しての必然性はないが このように設定すると扱い易い形となるので通常この条件を加える. また 常微分方程式の数値解法の教科書における RK 法の一般公式は式 (.. よりも次式によるものが多い. s K c ( K s s K しかしここでは 計算流体力学への応用における明快性から式 (.. を一般公式として採用する.

13 陽的ルンゲ クッタ法陽的 RK 法は に対して と設定するものであり s 段陽的 RK 法に対するブッチャー配列は以下となる. s s ss s s s c c c L L M O O M M M L L L L L L この場合 ( c として ( ( s c ( ( ( ( ( s と記述すると s 段陽的 RK 法は次のように表現できる ( 更新方式. ( ( ( ( ( ( ( s s c なおこの公式から c は 段目の計算時点での から の間の到達割合いを意味することがわかる. 以下 4 次精度までの陽的 RK 法の公式の具体例を挙げる. まず 陽的オイラー法は 段 次精度陽的 RK 法 (RK と見なすことができ 対応するブッチャー配列と計算式は 陽的オイラー法 (RK ( ( ( となる. これが陽的オイラー法 と一致することは容易に理解できる.

14 段 次精度陽的 RK 法 (RK 段陽的 RK 法の一般形は次式で与えられる. ( ( ( ( ( ( ( ( c ここで を自由パラメータに選ぶと 段 次精度陽的 RK 法 (RK に対する係数のブッチャー配列は以下で与えられる. c /( /( RK の一般係数 ( は自由パラメータ 特に の場合を修正オイラー法 / の場合を He 法 ( 改良オイラー法 と呼ぶ. また / 4 の場合は Rlso の最適公式である. 修正オイラー法 ( He 法 ( 改良オイラー法 / Rlso の最適公式 ( / 4 / / / / / / / 4 / 4 ( ( RK ( Modfed Eler ( ( ( ( / / ( ( 図.5 RK( 修正オイラー法 による時間進行 4

15 補足 :RKの係数まず ( ( の時間微分は d d D d d d d d ( ( D D ( 4 D D D d D ( D と書ける. ただし ( ( / および D D D とする. これより ( のテイラー展開は次式で記述される. [ ] ( ( D D ( D! 4 D ( D ( D ( D D 4! [ ] L 次に RKの一般公式における ( ( を すると次式を得る. ( ( c ( ( ( 4 および ( c c L ( ( これをRKの時間進行の式 ( ( ( ( ( c L まわりにテイラー展開 に代入すると RK は の計算を行っていることに対応する. ただし / である. これを ( のテイラー展開 と比較して O( まで係数を一致させ さらに式 (.. の条件も加えると RK の係数に対する条件式を得る. c c / これは自由度 の連立方程式 ( 式 つ 変数 4つ であるので 係数の解を表現すると c /( を得る. を自由パラメータとして 5

16 段 次精度 RK 法 (RK 段陽的 RK 法の一般形は次式で与えられる. ( ( ( ( ( ( ( ( ( ( ( ( ( ( c c 段陽的 RK 法に対する係数のブッチャー配列は以下となる. c c 代表的な 段 次精度 RK 法 (RK に対する係数のブッチャー配列を挙げておく. 古典的 RK( / c He のRK( / 4 c / / / / 4 / 4 / / 9 8 / 9 / 6 / / 6 / 4 / 4 Rlso のRK( / c / 4 Wry のRK( 8/ 5 c / / / 8/5 8/5 / 4 / 4 / / 4 5/ / 9 / 4 / 9 / 4 / 4 となるRK( / 6 c / 4 Wllmso のRK( / c / 4 / 6 / 6 / / / 4 7 / 4 / 4 /6 5/6 / 7 4 / 7 / 6 / 8/5 これらのうち となるRKは更新方式における低容量型 RKである. つまり によって の計算に ( が必要なくなるので 更新計算において他のRK 公式よりも配列を つ減らせる.Wry のRKは積み上げ方式における低容量型 RKである.Wllmso のRK は Wllmso の計算法に従えば低容量 (N 容量 型となるRKである. 6

17 補足 :RK の係数 RK と同様にして RK に対する係数の条件式も得られる. c c c c / c c / c / 6 これは自由度 の連立方程式 ( 式 6つ 変数 8つ である. c と c を自由パラメータとして係数の解を表現すると c c c c c / に対して c c c ( c c c c ( c c 6cc c c c c c c 6c ( c c c 6c ( c c を得る. また を自由パラメータとして c c / の場合に対する係数の解 / / / / /(4 / 4 /(4 / 4 および c / c の場合に対する係数の解 / / /(4 / 4 /(4 / 4 が知られている. 7

18 4 段 4 次精度 RK 法 (RK4 4 段陽的 RK 法の一般形は次式で与えられる. ( ( ( ( ( ( ( ( 4 ( ( ( ( ( ( ( ( ( ( 4 ( ( ( 4 4 c c c 4 4 段陽的 RK 法に対する係数のブッチャー配列は以下となる. c c c 次精度を満足する 4 段陽的 RK 法の係数の一般解はいくつか知られているが ここでは代表的な 4 段 4 次精度 RK 法 (RK4 に対する係数のブッチャー配列を挙げておく. 古典的 RK4 K のRK4(K の /8 公式 / / / / / / / / / 6 / / / 6 / 8 / 8 / 8 / 8 Gll のRK4 / / / ( / ( / / ( / / 6 ( / 6 ( / 6 / 6 補足 :RK4 の係数 RK あるいは RK と同様にして RK4 に対する係数の条件式も得られる. c c c c c 4c4 / 8

19 c c 4c4 / c c 4c4 / 4 c 4 ( 4c 4c / 6 c 4 4c 4c / cc 4 4c 4c c4 / 4 4c / 4 ( ( 8 これは自由度 の連立方程式 ( 式 変数 である. まず c と 係数の解を表現すると c を自由パラメータとし c ( c ( c ( c c c( c c 4c ( c c c [ c c c ] c の場合についてRK4の c 4 c( c c c c ( c ( ( ( c ( c ( c 4 4 c ( c c [ 6cc 4( c c ] c( c c [ 6cc 4( c c ] ( c c c c cc c ( c c ( c c( c c ( c また ( ( を得る. この解は c / c / で K の RK4を与える. c c / に対しては を自由パラメータとする次のパラメータ解が知られている. / / / / /(6 / 6 /(6 / / 6 この解は / で古典的 RK4 ( / 6 で Gll のRK4を与える. 他にも4つのパラメータ解が知られているが そのうちのつを示しておく ( ただし. / / 8 /(4 / 6 / 8 /( / 6 /( / / / /( / 6 / 6 /( / / 6 / 6 9

20 陰的ルンゲ クッタ法陰的 RK 法はある > について なので 基本的に式 (.. をそのまま適用する必要があ る. 例えば 段陰的 RK 法に対するブッチャー配列と計算式は c c および ( ( ( { ( c ( c } ( ( ( ( c ( c } ( ( { ( c ( c } であり ( と ( の連立方程式を解いた後に を計算する. 一般に s 段陰的 RK 法では ( ( ( s L の連立方程式をニュートン法等の繰り返し法で時間ステップ毎に解く必要があり計算量が膨大となる. 半陰的 RK 法は > に対して と設定するものであり 例えば 段半陰的 RK 法に対するブッチャー配列と計算式は c c および ( ( ( c ( ( ( { ( c ( c } ( ( { ( c ( c } ( となる. この場合は および ( に対する単独の陰的方程式を順番に解いた後に を計算すれば良いので 段あたりの計算負荷はAM 法と同程度となる. 陰的 RK 法は s 段で最大 s 次精度まで 半陰的 RK 法では s 段で最大 (s 次精度まで到達可能である.

21 陰的 RK 法 (IRK 法 および半陰的 RK 法 (SIRK 法 の例を挙げておく. 陰的オイラー法 ( 段 次 IRK 法 陰的中点法 ( 段 次 IRK 法 CN 法 ( 段 次 SIRK 法 / / / / / / 段 4 次 IRK 法 段 次 SIRK 法 ( ノルセット法 ( / 6 / 4 ( / ( / 6 ( / 6 ( / 6 ( / / 4 ( / 6 / ( / / / / / 6 段 6 次 IRK 法 ( 5 5 / ( 5 5 / / 5 / 6 5 / 6 5 / 6 5 /8 5 / 4 5 / / 9 / 9 / 9 4 / 9 5 /5 5 /5 5 / 6 5 / 6 5 / 6 5 /8 5 / 5 / 4 ここで 陰的オイラー法 陰的中点法 クランク ニコルソン法 (CN 法 について補足しておく. 陰的オイラー法は 段 次精度陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ ( 陰的オイラー法 および ( ( ( ( ( ( と記述できる. この場合のつの計算式は等価で であり これを用いると AM 公式での陰的オイラー法 が得られる. 陰的中点法は 段 次精度陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ / / ( 陰的中点法 および ( ( ( ( と記述できる.つの計算式を比べると ( / となっており これが陰的中点法

22 であることがわかる. クランク ニコルソン法 (CN 法 は 段 次精度半陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ および / / ( / / ( クランク ニコルソン法 ( ( ( ( ( ( ( ( ( ( と記述できる. この場合 ( であり これが AM 公式でのクランク ニコル ソン法 ( であることがわかる.

23 時間進行法の安定性解析 線形テスト問題 d λ λ λ R λi d ( に対する時間進行法の数値安定性を考える. まずこの問題の解析解は λ λ ( R λ e e e I で与えられ λr ならば > で ( は有界となる. λ I Sle Usle λ R ω ( k s ( ω( k x x 図.6 線形テスト問題の安定領域 ( 複素平面 補足 : なぜ線形テスト問題を考えるのか 流体運動の支配方程式は輸送方程式であるので その本質を表現するモデル方程式 ( 次元の 移流 拡散方程式の空間離散式 を考える. d δ δ U ν d δx δx ここで U : 移流速度 ν : 拡散係数 δ : δx x に対する空間離散オペレータ 解析的な扱いを容易にするため x 方向に周期的な場を仮定し 離散フーリエ変換を導入する. N k ( π ˆ ( x ex[ ω ( k x ] k N / ~N / ω( k k Lx Lx N x N / ( x ( ex[ ω( k x ] ~N ˆ k k N / 離散フーリエ変換によりモデル方程式を対角化する. dˆ k ( [ νω ( k Uω ( k ] ˆ k ( k k N / ~N / d ここで ω ( k と ω ( k は空間 階および 階微分に対する離散オペレータの修正波数で 空間離 散化手法に応じて定まる. 例えば フーリエ スペクトル法 : 次精度中心差分 : ( k ω( k ω ( k ω( k ω ω ( k s ( ω( k x / ( x / 以上より モデル方程式は線形テスト問題 ( 固有値問題 に帰着され 固有値の実部は拡散項に 固有値の虚部は移流項に関係することが理解される. dˆ k ( R I R I λ kˆ k ( λk λ k λ k λ k ν ω ( k λ k U ω ( k d なお 純粋移流問題の安定性は ω ( k の最大値 純粋拡散問題の安定性は ω ( k の最大値に対応して定まる.

24 線形テスト問題の解の有界性が数値計算においてどのような条件で満足されるかを検討するため まず陽的オイラー法を線形テスト問題に適用してみる. ( λ ( λ この操作は 初期値 数値解の複素増幅率をξ として ξ と記述するとき 陽的オイラー法の複素増幅率が ( z ξ ただし z λ となることを意味する. ここで z λ λr λi としている. これから明らかなとおり 陽的オイラー法による数値解が有界であるための条件は ξ z で与えられる. これを満足する z は複素平面において ( を中心とする半径 の円の内側に存在する. この領域を陽的オイラー法の絶対安定領域と呼ぶ. これは 陽的オイラー法を用いて安定な数値解を得るためには 時間刻み幅 の設定に制限があることを示している. 具体的には λ を負の実数とすると λ の範囲に時間刻み幅 を設定する必要がある. Eler exlc Usle λ I Sle λ R 図.7 陽的オイラー法の絶対安定領域 4

25 次に 陰的オイラー法を線形テスト問題に適用し絶対安定領域を調べる. ( λ これより 陰的オイラー法の複素増幅率 ξ は次式で与えられる. ξ z λ ( z これに対応する絶対安定領域は次式で与えられる. ξ z これを満足する z は複素平面において ( を中心とする半径 の円の外側の全領域である. この領域は λr となる複素平面の左半分を全て含むので 陰的オイラー法は λr に対し時間刻み幅 によらず無条件安定な数値解を与える. λ I Sle Eler mlc Usle Sle Usle Im(z Re(z O λ R 図.8 陰的オイラー法の絶対安定領域 図.9 A( 安定 陰的オイラー法のように絶対安定領域が複素平面の左半分 ( λ R 全てを含む時間進行法は A 安定 ( 絶対安定 という. なお 絶対安定の緩和表現として 絶対安定領域が複素 z 平面上で z ; < π rg z < < < π / を含む場合 A( 安定と呼ぶ. { } ( 無限扇形領域 ( ところで 陽的オイラー法および陰的オイラー法は 次精度の時間進行法であるため打ち切り誤差の影響が懸念される. そこで 次精度の時間進行法である中点蛙跳び法 (le-frog 法 を線形テスト問題に適用してみる. ( λ これに ξ および z λ を代入して整理すると 中点蛙跳び法の複素増幅率 ξ は 5

26 ξ z ξ の解として次式で得られる. ここでは複素増幅率 ξ の解を求めるのではなく 絶対安定領域の境 θ 界 ξ を与えるzの範囲を調べてみる. そのため上式を変形しξ e を代入すると ξ z ξ θ θ θ ( e e e θ e θ cos を得る. これは z が虚軸上の ( の範囲のみで安定 つまり中点蛙跳び法は λ R 以外では数値解が無条件に不安定となることを示している. 以上のように 時間進行法は必ずしも精度が高いものが適切なわけではなく 数値安定性も含めて検討する必要がある. 線形多段階法の絶対安定領域 線形多段階法の一般形 ( ただし β N N を線形テスト問題に適用し ξ および z λ を代入すると次式を得る. ξ z β ξ N N この式による複素増幅率 ξ の根が全て ξ を満足する場合に線形多段階法は絶対安定である. 線形多段階法の安定限界曲線 ( 絶対安定領域の境界 つまり ξ を与える z は LM 法の安定限界関数 z ϕ に対して ( ξ N N θ ξ β ξ LM 法の安定限界関数 ξ e を代入し からπの範囲のθ に対する z を描けば良い. 以下に AM 公式 z ϕ ξ を示しておく. AM 公式および BD 公式に対する安定限界関数 ( AB 公式 z ξ 陽的オイラー法 (AB ( ξ ξ z AB ξ 6

27 AM 公式 ( ξ ξ z AB ξ 6ξ 5 ( ξ 4ξ z AB4 55ξ 59ξ 7ξ 9 ξ ξ z 陰的オイラー法 (AM ( ξ z クランク ニコルソン法 (AM ξ ( ξ ξ z AM 5ξ 8ξ ( ξ 4ξ z AM4 9ξ 9ξ 5ξ BD 公式 ξ z 陰的オイラー法 (BD ξ ξ 4ξ z BD ξ ξ 8ξ 9ξ z BD 6ξ 4 5ξ 48ξ 6ξ 6ξ z BD4 4 ξ 7

28 AB 公式 AM 公式および BD 公式の安定限界関数に対して π の θ の値を代入して描いた安定限界曲線を以下に示す. θ e を代入し さらにから 線形多段解法 (AB 法 AM 法 BD 法 では 精度が上がると絶対安定領域が狭くなることがわかる. また これらの中で A 安定な時間進行法は 陰的オイラー法 (AM クランク ニコルソン法 (AM および 次精度後退差分法 (BD のみである. なお BDは A(86 o 安定 BD4は A(7 o 安定である. ξ λ I AB Eler exlc s AB s AB4 s s λ R 図. AB 法の安定限界曲線 λ I 4 s AM AM4 s AM s λ R s AM s Eler exlc 図. AM 法の安定限界曲線 8 BD s BD4 6 λ I 4 Eler exlc s 7 o 86 o s s s λ I 4 BD BD BD4 s BD s λ R 図. BD 法の安定限界曲線 λ R 図. BD と BD4 の原点付近の安定限界曲線 8

29 RK 法の絶対安定領域線形テスト問題に s 段 RK 法を適用し行列表示すると次式を得る. ただし z λ である. ( z z L L z ( s z z z L z ( ( s M O O O M M M M z s s L zs zs s z ( ( ss z s s L L zss z ( ( ss z L L L z s この連立方程式を Crmer の公式で解き ( 解の s 成分 を求めると次式を得る. T ( I za z de R( z R( z e de( I za ここで s s 行列 A と I の成分を ( A および ( I δ また s 元ベクトル とe の成分 T T を ( L s および e ( L としている. R(z は RK 法の安定性因子と呼ば れ RK 法の絶対安定領域は R( z で与えられる. 補足 :RK 法の安定性因子の導出まず RK 法に対するブッチャー配列を c A T と表現し s s 行列 A と I の成分を ( A および ( I δ また s 元ベクトルの成分を T T ( L s c ( c L c s T T T ( ( s e ( L ( L ( L とすると 線形テスト問題に適用された s 段 RK 法は次のように記述される. ( I za e T z 上記連立方程式を Crmer の公式で解いた の解 ( 解の s 成分 はまず次式となる. ( I za e de T z ( I za de T z ここで ブロック行列の行列式に関する数学公式 A C de de( D de( A CD B B D を適用すると 分母および分子に現れている行列式はそれぞれ ( I za ( I za e T de T ( I za z de de T ( I za ze z de となるので T T de( I za ze de( I za ze R( z R( z de( I za de( I za を得る. 9

30 陽的 RK 法に対しては で なので de ( I za であり 安定性因子 R(z は T R( z de( I za ze 陽的 RK 法の安定性因子で与えられる. これは一般にzの s 次多項式 s R ( z qz qz L q s z となり 陽的 RK 法では A 安定なものが存在しないこと ( z で R(z となるため を示している. また 線形テスト問題に対する の解析解のテイラー展開 m λ z z e z L L m!! と比べると s 段 次精度陽的 RK 法の安定性因子は z z s R ( z z L q z L qs z!! となることが解る. 特に s 4 に対しては s 段 s 次精度陽的 RK 法が構成できるが この場合の安定性因子はそれぞれ以下となる. R ( z z 陽的オイラー法 (RK z R ( z z! RK z z R ( z z!! RK z z z R ( z z RK4!! 4! 4 陽的 RK 法に対する安定限界曲線を図.4 に示す. 線形多段階法と異なり 陽的 RK 法は精度が上がると絶対安定領域が広がる. RK4 λ I RK s s s Eler exlc RK s 4 λ R 図.4 陽的 RK 法の安定限界曲線

31 陰的および半陰的 RK 法の安定性因子 (z R は z の有理関数となるので A 安定となる可能性がある. 節で紹介した陰的および半陰的 RK 法は全て A 安定であり それぞれに対する安定性因子は以下となる. z z R ( 陰的オイラー法 ( 段 次 IRK 法 ( z z z R 陰的中点法 ( 段 次 IRK 法 および CN 法 ( 段 次 SIRK 法 6 6 ( z z z z R 段 次 SIRK 法 ( ノルセット法 ( z z z z z R 段 4 次 IRK 法 ( z z z z z z z R 段 6 次 IRK 法輸送法方程式の時間進行を考える場合 移流項の数値安定性は固有値の虚部に また拡散項の数値安定性は固有値の実部に支配されるので 時間進行法の安定性のまとめとして 実数および純虚数の固有値に対する安定領域を表. に示しておく.

32 表. 実数および純虚数の固有値に対する安定領域 (ABAMBDRKIRK AB(Eler exlc AB AB AB4 AM(Eler mlc AM(Crk - Ncolso AM AM4 BD BD BD4 RK RK RK4 Rel ( λ Imgry - λ - λ λ λ -. λ λ Sle Sle Sle Sle -6 λ - λ Sle Sle Sle R R R R R R R - λ R Usle( λ I Usle( λ Usle( λ I Usle( λ Sle I λ.94 I λ 4.7 I Usle( λ -.5 λ λ R R -.79 λ λ Remrks A - sle A - sle A - sle A(86 - sle A(7 - sle - sge IRK Sle Sle A - sle - sge IRK4 Sle Sle A - sle - sge IRK6 Sle Sle A - sle I I I I I I o o なお 時間進行法の安定性に関して次の定理が知られている. ダルキスト (Dhlqs の定理 : 陽的ルンゲ クッタ法は A 安定にはならない. 陽的な線形多段解法は A 安定にはならない. A 安定な線形多段階法の精度は を超えない. 4 次精度で A 安定な線形多段解法のうち 局所誤差の絶対値が最も小さいものはクランク ニコルソン法 (AM である. 第 章の参考文献 ( 三井斌友 数値解析入門 常微分方程式を中心に - 朝倉書店 (985. ( 三井斌友 常微分方程式の数値解法 岩波書店 (. ( U.M. アッシャー L.R. ペツォオルド共著 中森眞理尾雄監修 嘉村友作 三ツ間均共訳 常微分方程式と微分代数方程式の数値解法 培風館 (6. (4 Ger C.W. Nmercl Il Vle Prolems Ordry Dfferel Eqos Prece-Hll (97. (5 l S.O. Nmercl Mehods for Il Vle Prolems Ordry Dfferel Eqos Acdemc Press (988. (6 Bcher J.C. Nmercl Mehods for Ordry Dfferel Eqos Joh Wley d Sos (. (7 Wllmso J.H. Low-Sorge Rge-K Schemes J. Com. Phys. Vol.5 (

33 . 非圧縮性流れの計算アルゴリズム 目的 : 非圧縮性流れの計算アルゴリズム ( 圧力ベース解法 を理解する. 非圧縮性流れの運動はナビエ ストークス (NS 方程式および連続の式を用いて記述される. ( ν s (. (. 上式中で は時間 は速度ベクトル は圧力 / 密度 ( 本章では以降圧力と呼ぶ ν は動粘度 は勾配演算子 s は外力等によるソース項のベクトルである. 動粘度は定数とする. また 式 (. の左辺第 項は移流項 右辺第 項は圧力勾配項 右辺第 項は粘性項と呼ばれる. 本章では 式 (. および式 (. を解いて速度 および圧力 の解を求める非圧縮性流れの一般的な解法を示す. 空間および時間離散化手法を定めれば式 (. および式 (. は離散化できるが 非圧縮性流れでは連続の式 (. がNS 方程式 (. に対する拘束条件となるので これらの式が同時に満足されるように速度 および圧力 の数値解を求める解法が必要となる. 速度 については式 (. を基に時間発展を行えば良いが 非圧縮性流れでは圧力 に関する時間発展方程式が存在しない. このため 通常は式 (. に発散 ( を作用させて得られる圧力のポアソン方程式 ( [ ν ( s] (. あるいはこれに関連した方程式と連続の式 (. とを基に圧力を求めながら式 (. の時間発展を行い速度 が求められる.NS 方程式 (. の時間進行法に応じてそれぞれ経済的な解法がこれまでに開発されている. 表. 非圧縮性流れの計算アルゴリズム NS 式の各項の時間進行法計算アルゴリズム移流項粘性項圧力項 MAC 法系統陽解法陽解法フラクショナル ステップ法系統陽解法陰解法 ( 陰解法 陰解法陽解法 SIMPLE 法系統陰解法陰解法 それぞれの計算アルゴリズムについて 節 節 節で説明する. 解法の説明に具体的な空間離散式が必要となる場合には 次元の流れ場を考える. その場合のNS 方程式 (. および連続の式 (. は次式で与えられる. ν s x x y x (.4 x y ν s y x y y x y (.4 x y (.5

34 上式中で は速度ベクトル の x y 方向成分 s x 方向成分である. 離散値の表現方法として ( x y 点上に配置された変数の離散値を例えば s y はソース項ベクトルs の x y と表記する. さらに 時間離散点 における を と表記する. は時間刻み幅である. 節で示されるとおり非圧縮性流体の差分解析ではスタガード格子の使用が推奨されるので 本章での各解法の説明は基本的にスタガード格子を用いて行われる. なお 有限体積法に基づいて空間離散化が行なわれる場合には 通常式 (. の代りに次式 T ( { ν [( ( ]} s (.6 が基礎式として用いられる. ( T は速度勾配テンソル ( の転置テンソルである. 粘性係数 ν が定数の場合には連続の式 (. から式 (. と式 (.6 とは等しくなる. 渦動粘度 ν を用いる乱流モデルの運動方程式は式 (.6 のν をν で置き換えたものになり またν が空間的な変 ν 数である場合も扱えるので 非圧縮性流れに対する運動方程式としては式 (. よりも式 (.6 の方が一般的である. しかし式 (.6 を用いると解法の統一的な説明が煩雑になるので 節 節 節の解法の説明では式 (. を運動方程式として用いる. 乱流モデルを用いる場合および有限体積法により空間的離散化が行われる場合に対しては 粘性項の空間離散式が式 (.6 中の形を基に構成されているものと解釈して頂きたい. また 輸送量 φ の輸送方程式の一般形は次式となる. φ ( φ ( Γφ φ sφ (.7 ここで Γ φ はφ の拡散係数 s φ はソース項である. 式 (.7 の形の輸送方程式を流れ場と同時に解く場合 通常時間進行の各ステップにおいて速度 および圧力 が求められた後に式 (.7 を単純に時間進行してφ が求められるので 式 (.7 の時間進行法については特に説明しない. 差分格子系とMAC 法系統の解法 NS 方程式 (. の移流項と粘性項を陽的に時間離散化し 連続の式 (. を新たな時刻 において満足させると次式を得る. ( ν s (.8 (.9 式 (.8 では時間に関して 次精度の陽解法が用いられているが MAC 法系統の解法で時間精度の向上を図りたい場合には移流項および粘性項に対してより高次精度の陽的時間進行法 ( 例えば 次精度アダムス バシュフォース法 が用いられる. 式 (.8 に発散 ( を作用させて式 (.9 を用いると圧力 のポアソン方程式が得られる. [ ( s] ν (. MAC 法系統の解法では 基本的には圧力のポアソン方程式 (. を解いて式 (.9 を満足する速度 場 を式 (.8 から求める. 4

35 レギュラー格子を用いる差分法の問題点速度成分 および圧力 を同じ空間離散点上に配置する差分法は基礎方程式の離散化が 最も容易である. この空間配置による差分格子をレギュラー格子と呼ぶ ( 図.. 式 (.8 および式 (.9 を 次元で記述し 空間微分をレギュラー格子上の差分で近似して整理すると次式を得る. y y y : y y y : : : y y y y y x x x x x y x x x x x x x 図. レギュラー格子 図. レギュラー格子の問題点 f x x (. g y y (. x x y y (. ここで f および g は f ν x y x y g ν x y x y s x s y であり f と g はそれぞれ ( x y 点上での f と g の離散近似式である. 圧力 は 式 (. により および が時間進行された後に式 (. を満足するものでなくてはならない. そ のような圧力 は次式を満足する. 5

36 6 ( R y y g g x x f f D y y y y y y x x x x x x (. ここで (. R y y x x D 式 (. は式 (. を式 (. に代入して得られる. しかしながら式 (. の左辺は隣り合う格子点の値を参照しない離散式であるため ( 図. この式から得られる圧力の数値解は空間的に つ飛びの振動を起こす可能性があり滑らかな物理的な圧力解が得られない場合がある. 式 (. は式 (. を空間的に離散化したものと見なせるので 式 (. の左辺を 階微分に対する次の差分式で置き換えて圧力を求めることも考えられる. ( R y y g g x x f f D y y y y y y x x x x x x (.4 上式による圧力 の解は滑らかで物理的なものとなるが 逆にこの圧力を用いて計算される速度 および は連続の式の離散式 (. を満足しなくなる. つまり レギュラ格子を用いた差分法による非圧縮性流体の解法では 圧力の解の物理性あるいは連続の式の満足度のいずれかが犠牲になる.

37 MAC 法とスタガード格子レギュラー格子の問題点を解消する差分格子として 非圧縮性流体の数値解析ではスタガード格子が用いられる ( 図.. x / x / x / x / x / y / y y / y / y y y / y y / y x x x x x : : : 図. スタガード格子 スタガード格子では 圧力の配置点に対しそれぞれの方向に半格子ずれた位置に速度成分が配置される. つまり 圧力 の配置点を ( x y とすると x 方向速度成分 は ( x / y 点 y 方向速度成分 は ( x y / 点上に配置される.NS 方程式 (.8 の各方向成分はそれぞれの速度成分の配置点上で また連続の式 (.9 は圧力の配置点で空間的に離散化される. / / f / x / (..5 / / g / y / (..5 / x / / y / (..6 式 (.5 を式 (.6 に代入すると次式が得られる. 7

38 x y ( D x S / f y / / f x x / / y g / / g y / (..7 ここで ( / / / D /. S x y 式 (.7 の左辺は式 (.4 の左辺と等価であり これを解いて求められる圧力の数値解 は滑 らかで物理的であり さらにこの圧力を用いて式 (.5 により時間発展される速度 および は式 (.6 を離散的に満足する. スタガード格子の使用が提案されたこの解法はMAC 法 (Mrker Ad Cell 法 と呼ばれる (. 原論文においてはマーカーを飛ばして自由表面の位置を決定しながら流れ場を解く解法としてMAC 法が提案されているが 現在では流れ場の解法の部分のみに対してMAC 法の名称が用いられている. MAC 法の計算手順は以下のとおりである. 初期条件として と を与え 式 (.7 を解いて を求めた後 式 (.5 から と を求め 4 と を と として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. 8

39 SMAC 法 MAC 法では圧力のポアソン方程式 (..7 の右辺の境界値および圧力の境界条件を適切に与えることが一般に煩雑である. これを解消する解法にSMAC 法 (Smlfed MAC 法 がある (.S MAC 法による解法を説明するため まず式 (.5 を次のように分解する. また ˆ / / f / x / ˆ / / g / y / (.8 (.8 / ˆ / / ˆ / δ δ δ x y δ / δ / (.9 (.9. (. û および ˆ は中間速度 δ は圧力補正量である. 式 (.9 を式 (.6 に代入するとδ のポアソン方程式が得られる. δ x x ( Dˆ S δ δ δ δ δ δ δ / x / y y / y / (. ここで ( D ˆ ˆ ˆ ˆ S / / / /. x y ˆ SMAC 法の計算手順は以下のとおりである. 初期条件として を与え 式 (.8 から û ˆ を計算し 式 (. を解いてδ を求めた後 4 式 (.9 から と および式 (. から を求め 5 を として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SMAC 法 ( およびMAC 法 ではポアソン方程式を解くために計算時間の大部分が費やされるので 楕円型方程式の離散式が経済的に解ける行列解法を使用する必要がある.MAC 法と比べた場合のSMAC 法の利点は境界条件の設定がわかり易い点にある. 一例として 境界 BC (x N/ y での速度の値 N / が与えられている場合を考える ( 図.4. この場合 まず中間速度 û に対する境界値として 9

40 BC ˆ N / N / (. BC N / N / を与える. 境界条件は であることを要求するので 式 (.9 から δ (. N δ N が与えられる. なお SMAC 法の原論文 ( ではポテンシャル関数 ϕ が導入されているが δ とϕ は ϕ δ で結び付けられる. また 境界条件も含 めて空間離散化がなされた離散式を基に解法が適切に構成されていれば SMAC 法は MAC 法を分離して解いているだけで両者は同じ解法となる ( 移流項および粘性項の双方が陽的に時間進行されているので分離誤差は生じない. y x x N / x N x N/ x N/ BC N/ y δ N δ N 図.4 SMAC 法での境界条件 補足 : 行列分解によるSMAC 法の計算アルゴリズムの導出速度 ( 圧力 および ( f g の離散値ベクトルをそれぞれ およびf と表記すると 式 (.5 および式 (.6 は次式のように表現される. Gδ G f D ただし 圧力は δ としている. またG は離散的勾配演算子 D は離散的発散演算子であり これら離散的空間演算子では境界条件も考慮されている. 境界条件から生じる項はN S 方程式および連続の式それぞれに対して および に含まれているものとする. これらをまとめて行列表示すれば次式となる. I G ( G f D δ ここで 左辺の行列は次のようにLU 分解できる. 4 ( G f I I G D DG I δ これはさらに次の 式へと分解できる. I D ( G f ˆ DG δˆ I G ˆ I δ δˆ これらの行列を展開して整理すると次の計算アルゴリズムを得る. ˆ G f ( ( δ DG D ˆ ˆ Gδ これは SMAC 法の計算アルゴリズムを与える. ここで 第 式中の DG は離散ラプラス演算子に対応する.

41 4 HSMAC 法 ( またはSOLA 法 MAC 法およびSMAC 法の計算時間の大部分はポアソン方程式の離散式を解くために費やされる.HSMAC 法 (Hghly Smlfed MAC 法 ( またはSOLA 法 ( は δ のポアソン方程式 (. を解かずに速度と圧力の同時補正計算により連続の式を満足させるSMAC 法の近似解法である. まず 式 (. の左辺の対角項 ( δ に関する項 のみを残してδ の式を作り ( 優 対角近似 緩和係数 β を乗じると次式を得る. β ( Dˆ S δ (.4 x x / x/ y y / y / 式 (.9 および式 (. から理解されるとおり 点 ( 上の圧力補正量 δ のみが与えられた場合 点 ( まわりの速度 / がその影響を受ける. これより 点 ( る. および ˆ / ˆ / δ x / / / 4 / および圧力 のみ まわりの速度および圧力をそれぞれ次式により補正す (.5 (.5 ˆ / ˆ / δ x/ (.5c ˆ / ˆ / δ y / (.5d ˆ / ˆ / δ y / δ (.6 HSMAC 法の計算手順は以下のとおりである. 初期条件として を与え 式 (.8 からû ˆ を計算し 式 (.4 式 (.5 式 (.6 による点 ( まわりの速度および圧力の補正計算を全格子点に対して実行し さらにそれを MAX ( D ˆ S < ε が満足されるまで繰り返す (ε は連続の式に対する許容誤差. MAX ( D ˆ S < ε が満足された時点の速度と圧力を とする. 4 を として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. なお式 (.4 から 連続の式 ( Dˆ が満足されれば速度と圧力は補正 されないことがわかる. S MAC 法およびSMAC 法と比べた場合のHSMAC 法の利点は圧力のポアソン方程式を解く必要が無いという点にある. しかし HSMAC 法はSMAC 法の近似解法であり解の収束が式 (..4 中の緩和係数 β の値に大きく依存することが欠点となる. なお β には < β の範囲の値が用いられ 経験的に β. 7 程度の値が推奨されている (.

42 . フラクショナル ステップ法先のMAC 法系統の解法ではNS 方程式の粘性項および移流項が陽的に時間進行されているので数値安定性条件から時間刻み幅に対する制限が生じる. 特に壁面近傍に差分格子を集中させて速度の滑り無し条件を課す場合 MAC 法系統の解法では粘性項に対する数値安定性条件が厳しくなり非常に小さな時間刻み幅を使用する必要が生じる. 数値安定性条件の問題を解決するには陰的な時間進行法を導入すれば良いので ここではNS 方程式の粘性項に対してクランク ニコルソン法を適用する解法を考える. この場合 式 (. と式 (. を時間離散化した方程式は次式となる. ν ( ( ν s (.7 (.8 NS 方程式の移流項には 次精度のアダムス バシュフォース法が適用されている. ところで 式 (.7 および式 (.8 を直接解いて速度 および圧力 の連立解を得るには大規模な行列計算が必要となるので経済的ではない. これはNS 方程式の粘性項あるいは移流項を陰解に時間進行する解法に共通の問題点である. フラクショナル ステップ法は 粘性項と圧力項との分離に基づき非圧縮性流れの陰的な非定常数値計算を経済的に実施可能とする解法である. フラクショナル ステップ法と分離誤差まず フラクショナル ステップ法の一例として Km & Mo (KM ( によるフラクショナル ステップ法を示す. この解法では 式 (.7 は次の様に分離される. ˆ ν ν ( ( ˆ s (.9 ˆ ϕ (. 式 (.7 と式 (.9 との差をとり式 (. と比べるとポテンシャル関数 ϕ と圧力 との関係が ν ϕ ϕ (. で与えられることがわかる. 式 (. と式 (.8 からポテンシャル関数 ϕ のポアソン方程式が得 られる. ϕ ˆ (. KMのフラクショナル ステップ法の計算手順は以下のとおりである. 初期条件として を与え 式 (.9 から ˆ を計算し 式 (. を解いてϕ を求めた後 4 式 (. から を求め 圧力が必要ならば式 (. から を求め 5 を として に戻る. 以上の手順を必要な時間まで繰り返す. ところで KM のフラクショナル ステップ法では式 (.7 を式 (.9 および式 (. へと分離して解法が構成されているが このような分離の際には分離誤差が生じてしまう. 加えて 上記のように空間的に離散化されていない方程式を基にフラクショナル ステップ法を構成すると中 4

43 4 間速度 ˆ の境界条件の設定が曖昧となってしまう. 分離誤差および中間速度の境界条件設定の問題を解決するため 時間および空間に関して完全に離散化された方程式を基に解法を議論する. ここでは式 (.7 と式 (.8 の空間微分演算子についても離散化した方程式を考える. s c c L L G ν ν (. D (.4 ただし 速度 ( および圧力 の離散値ベクトルをそれぞれ および 移流項の離散ベクトルを c と表記している. また L は離散的ラプラス演算子 G は離散的勾配演算子 D は離散的発散演算子である. これらの離散演算子は境界条件も考慮されたものであり 境界条件から生じる陽的な項は NS 方程式および連続の式それぞれに対して および で表現されている. 式 (. と式 (.4 を行列を用いて表示すれば次式となる. KM KM D G A r r (.5 ここで A は L I A ν I は単位行列であり また式 (.5 の右辺のベクトルは s c c r r KM KM L I ν で定義される. 式 (.5 の右辺は陽的に扱う項および境界条件により生じる項 左辺が陰的に扱う項である. フラクショナル ステップ法の導入に伴う分離誤差の評価を行うため 式 (.5 に対して次の近似式を導入する (. ( KM KM D AB G A r r (.6 上式中の行列 B は A に対する近似行列である. 式 (.6 の左辺の行列は次式のように LU 分解できる. KM KM I BG I DBG D A r r (.7 上式はさらに次の 式へと分離できる. KM KM DBG D A r r ˆ ˆ (.8 ˆ ˆ I BG I (.9 式 (.8 および式 (.9 の行列を展開して整理すると次式を得る. A KM r ˆ (.4 ( ˆ D DBG (.4 ˆ BG (.4 式 (.4~ 式 (.4 を解くことは式 (.6 を解くことと等価である. また 式 (.4~ 式 (.4 中に

44 現れる空間的な離散演算子は既に境界条件が考慮されたものであるので 適切に離散化されたフラクショナル ステップ法では中間速度 û に対する境界条件は必要ないことがわかる. 次に 行列 B の選択に応じたフラクショナル ステップ法を考える. まず B A の場合は式 (.5 を直接解くことに等しく分離誤差はゼロであるが A の計算に時間がかかるのであまり用いられない (Uzw 法.KMのフラクショナル ステップ法を含め従来のフラクショナル ν ステップ法の多くは B I を採用することに対応し その場合 ( AB I L となるので 時間 次精度の分離誤差を持つ.KMのフラクショナル ステップ法は一見して時間 次精度の ようにも受け取れるが 式 (.7 を分離して式 (.9~ 式 (.5 を構成する際に ( ϕ ( ϕ という解析的な互換性が用いられている. しかし 離散式について対応する互換性は一般的に成立しない ( GLϕ LGϕ.Pero ( は時間 次精度および 次精度のフラクショナル ステップ ν 法がそれぞれ B I L および ν ν B I L L で構成できることを示した. しかし これらの高次精度のフラクショナル ステップ法はあまり汎用的ではない. Dkowcz - Dsky (DD のフラクショナル ステップ法非定常問題に対してフラクショナル ステップ法を適用する場合には できれば分離誤差も含めて 次精度以上の時間精度を持つ解法を使用したい.Pero の時間 次精度フラクショナル ステップ法とは別に Dkowcz &Dsky (6 (DD は式 (.5 をデルタ形式 ( 解の補正量を変数とする形式 に書き直した次の行列方程式を基に より汎用的な時間 次精度のフラクショナル ステップ法が構成できることを示した. A G r DD (.4 D r DD ここで式 (.4 の右辺のベクトルは以下で定義される. r DD G νl c c s r DD D DDのフラクショナル ステップ法が時間 次精度となることは式 (.4 の圧力項がデルタ形式であることによる. 速度ベクトルもデルタ形式で記述されているが これは後に式 (.54 で紹介する近似因子化を導入するためであり DDのフラクショナル ステップ法の本質とは関係ない.D Dのフラクショナル ステップ法では 式 (.4 の左辺の行列に対してKMのフラクショナル ステップ法と同様な近似を導入する. A AG r DD (.44 D r DD 式 (.4 から式 (.44 への近似の誤差は ν AG( G( LG( O( (.45 となるので DDのフラクショナル ステップ法は時間 次精度の分離誤差を持つ. ここで ( O と評価されている. 近似式 (.44 は次のようにLU 分解が可能である. A D I DG G I r r DD DD (.46 44

45 これは次のように分離できる. A δ r D DG δ r DD DD (.47 I G δ I δ (.48 式 (.47 および式 (.48 を展開して整理すると次式を得る. ν I L δ r DD (.49 ˆ δ (.5 δ DG ( D (.5 ˆ Gδ (.5 δ (..5 式 (.49~ 式 (.5 を用いて速度および圧力の時間進行を行うのが DD のフラクショナル ステップ法である. DD のフラクショナル ステップ法の計算手順は以下のとおりである. 初期条件として および を与え 式 (.49 を解いてδ を求め 式 (.5 からû を計算し 式 (.5 を解いてδ を求めた後 4 式 (.5 から を求め 式 (.5 から圧力 を求め 5 および をそれぞれ および以上の手順を必要な時間まで繰り返す. として に戻る. DDのフラクショナル ステップ法の計算手順のうち 式 (.49 は速度補正量の離散式 式 (.5 は圧力補正量のポアソン方程式に対応する離散式であり いずれも行列を解く必要がある. 式 (.5 は連続の式の離散式を満足させるために解かれるものであり そのままの形で出来るだけ精度良く解かねばならない. 式 (.49 に対しては近似因子化 ν ν ν I Lx I Ly I Lz δ r DD (.54 を導入すれば経済的に解くことができる. ここで L x L y L z はそれぞれ x y z 各方向 の 階微分に対応する離散演算子であり L と同様に境界条件も考慮されているものとする. 式 (.49 から式 (.54 への演算子の近似誤差は ( O であるので この近似は式 (.7 が持つ時間精 度 ( 次精度 を低下させない. 式 (.54 の解 δ は 次元の行列を 回解こくとで得られる. ν I L x δ r DD (.55 ν I Ly δ δ (.55 ν I Lz δ δ (.55c 45

46 SIMPLE 法系統の解法 NS 方程式の時間進行に関して MAC 法系統の解法では粘性項および移流項の双方が陽的に時間進行され フラクショナル ステップ法では粘性項が陰的に移流項は陽的に時間進行される. SIMPLE 法系統の解法は粘性項および移流項の双方を陰的に時間進行してさらに大きな時間刻み幅の設定を可能とし収束性および経済性を高める解法である. しかしながら 速度に関して非線形である移流項を完全に陰的に扱うことは各時間ステップの運動方程式の時間積分に繰り返し計算が必要となり経済的でないので SIMPLE 法系統の解法では式 (. と式 (. を時間離散化した次式 ( ν s (.56 (.57 を基に解法が構成される. 次元流れでは式 (.56 および式 (.57 は次式となる. x x x y y y x y ν sx (.58 x y ν s y (.58 x y (.59 さらに式 (.58 および式 (.59 の空間微分を離散化して整理すれば次式を得る. ( / / P / (.6 x ( P / / / / / / y / (.6 / x / / y / 式 (.6 中の総和記号は周辺点 E W N S の総和をとることを意味し (.6 の左辺は次式と展開される. P ( P / / / / [ ] P P E E W W N N S S / (.6 に関する式 E W N S は離散式の係数で 粘性項および移流項の離散化および境界条件に応じて定められる. ここで 変数の添字 PEWNS は離散点の相対的な位置を表す. 図.5 に / に対するNS 方程式のスタガード格子上の離散点を示す. 速度 / の離散 式については E / W / N / S / と定義され P / 46

47 る.SIMPLE 法系統の解法の空間離散化は 通常図中の破線で示されるコントロール ボリュームに関する有限体積法で行われる. 式 (.6 中の右辺第 項は圧力勾配項の離散式である. 式 / (.6 中の は離散式の残りの項をまとめたものである. 速度 / の離散式も同様に構成される ( 図.6. x / x / x / x / y N x / N x / y / y / y y y W P E y / y W P E y / y / y y S y x x S x y / x x x x x 図.5 P / に対する離散点図.6 / に対する離散点 P SIMPLE 法式 (.6 は速度成分 に関する離散式であるが 右辺に現れる が与えられなければ速 度 は求められない.SIMPLE 法 (7(Sem-Imlc Mehod for Pressre-Lked * Eqos では まず次式を解いて速度の予測値 * を求める. ( * / / * P / (.6 x ( P / / / * / * / / y / (.6 しかし * と * は連続の式の離散式 (.6 を満足しない. そこで圧力および速度の補正を行って式 (.6 を満足させることを考える. δ (.6 * / / δ / (.64 / * / δ / (.64 式 (.6 と式 (.64 を式 (.6 に代入した後に式 (.6 との差を計算すると δ およびδ とδ 47

48 との関係式を得る. ( P ( P / / δ δ δ (.65 / δ / x / δ δ δ (.65 / δ / y / ここで 解法が簡単になるように式 (.65 の左辺第 項の総和の項を消去してδ およびδ とδ との関係式を簡略化する. δ δ δ / ( (.66 P x / / δ δ δ / ( (.66 P y / / 圧力補正量 δ の方程式は式 (.64(.66 を連続の式 (.6 に代入して得られる. ここで x ( P x / ( P y * ( D S ( D * S / δ ( P y / ( P / δ δ δ / 48 δ δ x/ δ δ y / * * * * / / / / x y. / (.67 式 (.67 を解いて得られる圧力補正量 δ を用いて式 (.64(.66 により速度補正を行うと 連 続の式を離散的に満足する速度 および が得られる. 圧力の補正には式 (.6 ではなく緩和係数 ( を導入した次式を用いる. < δ (.68 緩和係数 の導入は式 (.65 から式 (.66 への簡略化の際に生じる誤差のために必要となる. ま た 速度場が連続の式を満足するように式 (.67 が定められているので速度の補正に緩和係数は導入しない. SIMPLE 法の計算手順は以下のとおりである. 初期条件として を与え 式 (.6 を解いて * と * を求め 式 (.67 を解いてδ を求め 4 式 (.64 および式 (.66 による速度補正 および式 (.68 による圧力補正を通して を求め 5 を として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す.

49 なお 式 (.6 のような輸送方程式の離散式を繰り返し法で解く場合 SIMPLE 法系統の解法では減速緩和法を導入して数値計算の安定性を高めることが行われる. 式 (.6 に減速緩和法を導入すれば次式を得る. re ( P / / * / x / * / ( P re / / (.69 ここで / は繰り返し計算におけるつ前のステップでの * の値であり は減速緩和係数 ( < である. 速度 に対しても また一般的に輸送方程式の離散式を解く場合にも同様の減速緩和法が導入される. SIMPLE 法の問題点は つのパラメータ( と を流れ場に応じて定めなければならないこと および解の収束性が の値に大きく依存することである.Pkr (8 は.5.8 を推奨しているが その後の研究で - とすれば概ね後に述べるSIMP LEC 法と同等の収束性が得られることが示されている (9. 補足 : 行列分解によるSIMPLE 法の計算アルゴリズムの導出式 (.6 式 (.6 をまとめて行列表示すれば次式となる. A G D ただし G は離散的勾配演算子 D は離散的発散演算子であり これら離散的空間演算子では境界条件も考慮されているものとする. 境界条件から生じる項はNS 方程式および連続の式それぞれに対して および に含まれているものとする. ここで δ を代入し 左辺の行列を LU 分解すると以下を得る. A I A G G D DA G I δ この式は次の 式へと分解できる. A * D DA G * δ * I A G * I δ δ G さらに A - に対して次の近似を導入する. A [ dg( A ] dg( は行列の対角項を抽出する関数で は上記対角化近似における補正係数である. 最後に δ * δ とおいた後に分解された行列を展開して整理すると次の計算アルゴリズムを得る. 49

50 A D * G * [ dg( A ] Gδ D * [ dg( A ] Gδ δ これは SIMPLE 法の計算アルゴリズムである. また 以上の導出から SIMPLE 法の圧力補正係数 の導入は A - の近似に伴うものであることが理解される. SIMPLER 法 SIMPLE 法の収束性を改善する解法として SIMPLER 法 (SIMPLE Resed が提案されている (8.SIMPLE 法の問題点は圧力補正における減速緩和にあるので SIMP LER 法では圧力補正量は連続の式を満足させるためだけに用い 圧力は補正計算をせずに運動方程式から直接求める. この手順を説明するため まず式 (.6 を変形した次式を導入する. (.7 / / / / ( P ( P x / (.7 / / ( P ( P y / / / 式 (.7 の右辺第 項は速度のみで計算される. また解が収束した時点でを考慮し 次の中間速度 û および ˆ を定義する. となること ˆ / / / (.7 ( P / / ˆ (.7 ( P / この中間速度は既知量である. また 式 (.7 の右辺第 項を û および ˆ で置き換え この変形に伴う圧力を とすれば次式を得る. / ˆ / / ( P x / 5 (.7 / ˆ / / ( P y / (.7

51 式 (.7 を連続の式の離散式 (.6 に代入すれば圧力 ここで x ( P x / ( P y ( Dˆ S ( D ˆ / ( P y / ( P / ˆ ˆ の方程式が得られる. / ˆ S / / / /. x y ˆ x / / y / (.7 * あとはSIMPLE 法と同様に を用いて速度の予測値 および * の方程式を解いた後に圧力補正量 δ の方程式を解き速度補正を行う. ただしその際に圧力の補正は行わない. SIMPLER 法の計算手順は以下のとおりである. 初期条件として を与え 式 (.7 から中間速度 û および ˆ を計算し 式 (.7 を解いて圧力 を求め 4 式 (.6 を解いて * と * を求め 5 式 (.67 を解いてδ を求め 6 式 (.64 および式 (.66 による速度補正を通して速度成分 を求め 7 を として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SIMPLER 法では圧力方程式を 回 ( 式 (.67 と式 (.7 解く必要があるので時間進行 ステップあたりの計算時間は SIMPLE 法よりもかかるが 圧力の緩和係数は使用しないのでその最適値を探す必要がなく 定常計算における収束性は一般的に SIMPLE 法よりも優れている. なお 式 (.6 を繰り返し法で解く過程では式 (.69 と同様の減速緩和法が導入される. 5

52 SIMPLEC 法 SIMPLE 法において圧力の緩和係数 を導入しなければならないのは速度補正式の導出 の過程における式 (.65 から式 (.66 への大胆な近似に起因する. 速度補正式の導出の過程において SIMPLE 法よりは妥当な近似を用いて速度補正式を導出するのが SIMPLEC 法である (9.S IMPLEC 法の説明のため まず式 (.65 を次の様に書き直す. P P ( δ δ δ δ δ (.74 / P / / x / δ δ δ ( δ δ (.74 / P / / y / 速度補正量は空間に連続的に分布すると考えられるので式 (.74 の左辺第 項は他の項よりも値が小さいことが期待される. これより 式 (.74 の左辺第 項を消去してδ およびδ とδ との関係式が与えられる. δ δ δ / (.75 x / P / δ δ δ / (.75 y / P / 圧力補正量 δ の方程式は式 (.64 と式 (.75 を連続の式 (.6 に代入して得られる. x * ( D S P y P / δ δ x / / δ y δ / P P / δ δ x / / δ δ y / (.76 上式を解いて得られる圧力補正量 δ を用いて式 (.64 式 (.75 により速度補正を行うと連続の 式を離散的に満足する速度 および が得られる.SIMPLE 法とは異なり SIMP LEC 法では圧力の補正に対する緩和計算は必要なく式 (.6 を直接用いる. SIMPLEC 法の計算手順は以下のとおりである. 初期条件として を与え 式 (.6 を解いて * と * を求め 式 (.76 を解いてδ を求め 5

53 4 式 (.64 および式 (.75 による速度補正 および式 (.6 による圧力補正を通して を求め 5 を として に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SIMPLEC 法では SIMPLER 法と同様に圧力の緩和係数を導入する必要がないことが大きな利点である. また SIMPLEC 法で解かれる圧力補正量の式 (.76 および速度補正量の定義式 (.75 は SIMPLE 法で対応する式の係数が変更されているだけであるので SIMP LEC 法の時間進行 ステップあたりの計算時間は SIMPLE 法とあまり変わらず SIMPL ER 法よりは少なくなる. なお 式 (.6 を繰り返し法で解く過程では式 (.69 と同様の減速緩和法が導入される. 第 章の参考文献 ( Hrlow.H. d Welch J.E. Nmercl clclo of me-deede scos comressle flow of fld wh free srfce Phys. lds 8 ( ( Amsde A.A. d Hrlow.H. The SMAC mehod: A mercl echqe for clclg comressle fld flows Los Almos scefc lorory of he U. of Clfor Re. No. LA-47 (97. ( Hr C.W. d Cook J.L. Clclg hree-dmesol flows rod srcres d oer rogh err J. Com. Phys. ( (4 Km J. d Mo P. Alco of frcol-se mehod o comressle Ner-Sokes eqos J. Com. Phys. 5.9 ( (5 Pero J.B. A lyss of he frcol se mehod J. Com. Phys. 8 ( (6 Dkowcz J.K. d Dsky A.S. Aroxme fcorzo s hgh order slg for he mlc comressle flow eqos J. Com. Phys. ( (7 Pker S.V. d Sldg D.B. A clclo rocedre for he mss d momem rsfer hree-dmesol rolc flows I. J. He Mss Trsfer 5. ( (8 Pkr S.V. Nmercl he rsfer d fld flow (McGrw-Hll New York (98 ( 水谷 香月訳 コンピュータによる熱移動と流れの数値解析 ( 森北出版 (985.. (9 V Doorml J.P. d Rhy G.D. Ehcemes of he SIMPLE mehod for redcg comressle fld flows Nmercl He Trsfer 7 ( ( 日本機械学会 計算力学ハンドブック ( 第 Ⅱ 巻 (6 第. 節 (