将来行ってみたい研究

(1)

1 中部CAE懇話会：基礎講座第７回(Ｈ21年度)「流体伝熱基礎講座」 H21年度・４日目（２００９年８月２９日（土））

流体数値計算法

講師：名古屋工業大学・大学院工学研究科

教授森西洋平

(2)

１．時間進行法１・１線形多段解法 1・1・1 アダムス・バシュフォース公式 1・1・2 アダムス・モルトン公式 1・1・3 後退差分公式１・２ルンゲ・クッタ法 1・2・1 陽的ルンゲ・クッタ法 1・2・2 陰的ルンゲ・クッタ法１・３時間進行法の安定性解析 1・3・1 線形多段解法の絶対安定領域 1・3・2 ルンゲ・クッタ法の絶対安定領域２．非圧縮性流れの計算アルゴリズム２・１差分格子系とMAC法系統の解法 2・1・1 レギュラー格子を用いる差分法の問題点 2・1・2 MAC法とスタガード格子 2・1・3 SMAC法 2・1・4 SOLA法（HSMAC法）２・２フラクショナル・ステップ法 2・2・1フラクショナル・ステップ法と分離誤差 2・2・2 DDのフラクショナル・ステップ法２・３ SIMPLE法系統の解法 2・3・1 SIMPLE法 2・3・2 SIMPLER法 2・3・3 SIMPLEC法講義内容午前 ( 9:30-12:30）午後 (13:30-16:30）

(3)

3 １．時間進行法１・１線形多段解法 1・1・1 アダムス・バシュフォース公式 1・1・2 アダムス・モルトン公式 1・1・3 後退差分公式１・２ルンゲ・クッタ法 1・2・1 陽的ルンゲ・クッタ法 1・2・2 陰的ルンゲ・クッタ法１・３時間進行法の安定性解析 1・3・1 線形多段解法の絶対安定領域 1・3・2 ルンゲ・クッタ法の絶対安定領域２．非圧縮性流れの計算アルゴリズム２・１差分格子系とMAC法系統の解法 2・1・1 レギュラー格子を用いる差分法の問題点 2・1・2 MAC法とスタガード格子 2・1・3 SMAC法 2・1・4 SOLA法（HSMAC法）２・２フラクショナル・ステップ法 2・2・1フラクショナル・ステップ法と分離誤差 2・2・2 DDのフラクショナル・ステップ法２・３ SIMPLE法系統の解法 2・3・1 SIMPLE法 2・3・2 SIMPLER法 2・3・3 SIMPLEC法講義内容午前 ( 9:30-12:30）午後 (13:30-16:30）

(4)

( )

t

₀

u

₀

u

=

( )

₍

_{( )}

_{) (}

₎

0

,

u

t

F

dt

t

du

₌

_>

目的：１階常微分方程式に対する時間進行法を理解する．本章では変数 _{u(t) に関する次の}１階常微分方程式の初期値問題の時間進行法（数値時間積分法）を考える． (1.1) １階常微分方程式 (1.2) 初期条件つまり，時間刻み幅毎の時間離散点上の値，

u(t₀+

Δ

_{t), u(t}₀₊₂

Δ

_{t), u(t}₀₊₃

Δ

_{t), ・・・}

を，式(1.2)を初期条件として，式(1.1)を満足させながら求めたい．

(5)

5 時間進行法を構築するため，まず_u(t+

Δ

_{t)をu(t)まわりに} テイラー展開して式(1.1)を使用すると次式を得る． ( + Δ ) ( )= + Δ ( )+ Δ ( )+ Δ ( )+_L+ Δ n n _n( )+_L dt t u d n t dt t u d t dt t u d t dt t du t t u t t u ! ! 3 2 3 3 3 2 2 2 ( )+Δ ( ( ))+ Δ ( ( ))+ Δ ( ( ))+_L = 2 2 3 2 , ! 3 , 2 , dt t u t F d t dt t u t dF t t u t tF t u ( ) ₍ _{( )}₎ ( ) L L+ Δ ₋ + − 1 1 , ! n n n dt t u t F d n t 時間に関する離散点上の値を表現するため， (1.3)

( )

n n _u _t u =

( )

(

n n

)

n _F _t _u _t F = , t t t t n tn = Δ , n+1 = n + Δ の表記を使用する．以降も同様．

(6)

t u(t) t n Δt Euler explicit F n t n+1 u(t) u n u n+1 n n n _F 図1.1 陽的オイラー法による時間進行最も単純な時間進行法はテイラー展開の右辺第３項以降を打ち切って構成される陽的オイラー法である． t u u +1 = + Δ ⋅

(

n n

)

_Fn t u u ₌ Δ − +1 あるいは，陽的オイラー法では，_tn_における勾配_Fn_と_un_から_un+1_{を計算する．} 陽的オイラー法 F dt du = _に対する

(7)

7 式(1.1)と陽的オイラー法をテイラー展開を通して比べると，打ち切誤差の初項が_O(

Δ

t1_{)となり，}_{陽的オイラー法}_は_１次精度しか持たないことがわかる．時間進行法の精度を向上させる方法として，１）線形多段階法２）１段階法がしばしは採用される．

(8)

1)線形多段解法・_un+1_{を計算する際に，時刻}_tn_および_tn+1_{での情報のみならず，} さらに過去の情報（時刻tn-1_，_tn-2_{，… ）も利用して高次精度} 化を行う．・計算の出発点近傍で別の時間進行法を使用する必要があり，また時間刻み幅の変更が容易でないことに注意が必要．・代表は，アダムス・バシュフォース（AB）公式アダムス・モルトン（AM）公式後退差分（BD）公式 2) １段解法：・_un+1_{を計算する際に，基本的に時刻}_tn_{での情報（およびそれ} から計算される値）のみを利用して高次精度化を行うもの．・代表は，ルンゲ・クッタ（RK）法以下本章ではこれら時間進行法を取り上げ紹介する．さらに，これら時間進行法の数値安定性についても示す．

(9)

9 補足：式(1.1)の１階常微分方程式を扱う理由流体運動の支配方程式は一般に偏微分方程式である輸送方程式として記述されるが，なぜ１階常微分方程式を扱うのかを説明する．例として外力を伴う１次元移流拡散方程式を考える． f x u x u U t u ₊ ∂ ∂ + ∂ ∂ − = ∂ ∂ 2 2 ν 変数 u=u(t,x) が独立変数 t とx の双方に依存するので，この方程式は偏微分方程式である．またより一般的に，移流速度 U=U(t,x)，拡散係 数ν=ν(t,x) ，および外力 f=f(t,x) も時間と空間の双方に依存するもの とする．

(10)

ここで，空間に等間隔格子（x_i=iΔx, 0≦i≦N+1）を使用し，離散値を

u(t,x_i)=u_i(t)=u_i のように表現し，空間２階微分を２次精度中心差分で近似する．

(

i N

)

f x u u u x u u U dt du i i i i i i i i i ₊ _≤ _≤ Δ + − + Δ − − = + − + 2 − , 1 2 2 1 1 1 1 ν 空間離散化後の変数はu_i=u_i(t)なので偏微分が常微分に変わる．さらに変数を

u=[u₁,u₂,…,u_N]Tとベクトル表記し，端点i=0,Nに対する境界条件も含めて上式 を書き直すと次の離散システム方程式を得る． ( ) _A_{( ) ( ) ( )}_t _t _t _F₍_t _{( )}_t ₎ dt t d u b u u , = + = A(t)は差分の係数を要素とする(N×N)の行列，b(t)は外力と境界条件を含むベ クトルである．このように流体運動の支配方程式を空間離散化した後の時間発展方程式は基本的に式(1.1)の形となる．空間多次元の場合も同様である．なお，離散システム方程式の変数はベクトル，式(1.1)の変数はスカラーであるが，時間進行法の構成に違いは生じない．

(11)

11 １．時間進行法１・１線形多段解法 1・1・1 アダムス・バシュフォース公式 1・1・2 アダムス・モルトン公式 1・1・3 後退差分公式１・２ルンゲ・クッタ法 1・2・1 陽的ルンゲ・クッタ法 1・2・2 陰的ルンゲ・クッタ法１・３時間進行法の安定性解析 1・3・1 線形多段解法の絶対安定領域 1・3・2 ルンゲ・クッタ法の絶対安定領域２．非圧縮性流れの計算アルゴリズム２・１差分格子系とMAC法系統の解法 2・1・1 レギュラー格子を用いる差分法の問題点 2・1・2 MAC法とスタガード格子 2・1・3 SMAC法 2・1・4 SOLA法（HSMAC法）２・２フラクショナル・ステップ法 2・2・1フラクショナル・ステップ法と分離誤差 2・2・2 DDのフラクショナル・ステップ法２・３ SIMPLE法系統の解法 2・3・1 SIMPLE法 2・3・2 SIMPLER法 2・3・3 SIMPLEC法講義内容午前 ( 9:30-12:30）午後 (13:30-16:30）

(12)

1・1 線形多段階法

(N+1)段の線形多段階法（Linear multi-step methods）は一般に次式で記述される．ただし

α

₁

≠

_0． ∑ ∑ − = =− + + ₌ Δ 1 1 1 N j j N j n j j n ju F t α β 線形多段階法の一般形 F dt du = に対するもの特に，

α

₁_=1,

α

₀

=-1

_,

β

₁_=0,

α

_j (j

≦-

1)

＝０

とすれば, アダムス・バシュフォース(Adams-Bashforth)公式

α

₁_=1,

α

₀₌

-

_1,

β

₁

≠

_0,

α

_j (j

≦-

1)

＝０

とすれば, アダムス・モルトン(Adams-Moulton)公式

α

₁

≠

1,

β

₁=1,

β

_j (j

≦０

)

＝０

とすれば, 後退差分(Backward-difference)公式

(13)

13 t u(t) t n Δt Adams–Bashforth F n t n+1 u(t) u n u n+1 t n–1 F n–1 t n–2 F n–2 Δt Δt

(

_{) (}

₌ ₊ ₊ ₊_L

₎

Δ − ₋ − − − + 2 2 1 1 0 1 n n n n n F F F t u u _β _β _β 図1.2 アダムス・バシュフォース公式による時間進行アダムス・バシュフォース公式陽解法(

β

₁₌₀₎

(14)

図1.3 アダムス・モルトン公式による時間進行アダムス・モルトン公式

(

_{) (}

₌ ₊ ₊ ₊ ₊_L

₎

Δ − ₋ − − − + + 2 2 1 1 0 1 1 1 n n n n n n F F F F t u u _β _β _β _β t u(t) t n Δt Adams–Moulton F n t n+1 u(t) u n u n+1 t n–1 F n–1 t n–2 F n–2 Δt Δt F n+1 陰解法(

β

₁

≠

₀₎

(15)

15 図1.4 後退差分公式による時間進行後退差分公式

(

1 +1 0 −1 −1

)

₌ +1 Δ + + + n n _n n F t u u u α α _L α t u(t) t n Δt Backward–difference t n+1 u(t) u n u n+1 t n–1 t n–2 Δt Δt F n+1 u n–1 u n–2 陰解法(

β

₁

≠

₀₎ ※ 以下でそれぞれに対する４次精度までの公式を示す．

(16)

1・1・1 アダムス・バシュフォース公式

時間１次精度から４次精度までのアダムス・バシュフォース（AB）公式は以下のとおりである．

(

n n

)

_Fn t u u ₌ Δ − +1

(

1

_{) (}

1

₎

3 2 1 ₋ + − = Δ − n _n _n n F F t u u

(

1

_{) (}

1 2

₎

5 16 23 12 1 ₋ ₋ + + − = Δ − n _n _n _n n F F F t u u

(

1

_{) (}

1 2 3

₎

9 37 59 55 24 1 ₋ ₋ ₋ + − + − = Δ − n _n _n _n _n n F F F F t u u １次精度の公式（AB1）は陽的オイラー法である．特に精度が明示されていない場合，流体計算では２次精度の公式（AB2）をアダムス・バシュフォース法と称する場合が多い．陽的オイラー法(AB1） AB2 (アダムス・バシュフォース法） AB3 AB4

(17)

17 補足：AB公式の導出（１／３）

(

_{) (}

2

₎

₀ 2 1 1 0 1 = + + + − Δ − − − − − + L n n n n n F F F t u u _β _β _β とおき，時刻tn_まわりの_{テイラー展開}_，

(

)

L + Δ + Δ + Δ + = Δ − + 3 3 3 2 2 2 1 24 1 6 1 2 1 t dt F d t dt F d t dt dF F t u un n _n n n n ( ) ( ) ( ) _L, 0,1,2,_L 6 1 2 1 ₃ 3 3 2 2 2 = + Δ − Δ + Δ − = − _q _t _q dt F d t q dt F d t q dt dF F F n n n n q n から，AB公式が指定した時間精度（２次精度ならばO(Δt2_{)）で近似さ} れるように係数β₀, β_－1, β_－2，・・・を定める． AB公式を，

(18)

補足：AB公式の導出（２／３）例として２次精度アダムス・バシュフォース法（AB２）を取り上げる．

(

)

1 ₀ 1 0 1 = − − Δ − ₋ − + n n n n F F t u u _β _β AB２の係数

β

₀

,β

_－1 を定めるためテイラー表を利用する．

(

)

1 1 1 1 1 0 0 1 2 2 2 2 1 0 0 6 1 2 1 1 / − − − − − − − + − − − − + + + Δ − Δ Δ β β β β β β n n n n n n n F F t u u t dt F d t dt dF F O(Δt０_)とO(_Δ_t1_{)の列それぞれの係数の合計を０とする２つの条件から，} 係数β₀ ,β_－1 に対する連立方程式が得られる．

(19)

19 補足：AB公式の導出（３／３）係数β₀ ,β_－1 に対する連立方程式 ⎪⎩ ⎪ ⎨ ⎧ − = = + − − 2 1 1 1 1 0 β β β これを解きβ₀=3/2, β_－1=－1/2を得る．これはAB2公式を与える．

(

1

_{) (}

1

₎

3 2 1 ₋ + − = Δ − n _n _n n F F t u u なお，テイラー表の３列目（O(Δt2_{) ）の合計はAB2公式の}_{打切り誤差} の初項εを与える． AB2 (アダムス・バシュフォース法） ) ( 12 5 2 1 6 1 ₂ ₂ 2 2 2 2 2 1 t O t dt F d t dt F d n_Δ ₌ n_Δ ₌ _Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = β₋ ε

(20)

補足：陽解法による計算陽解法は既知量の単純代入計算のみで時間進行が行えるので計算コードの作成が容易であり，また陰解法に比べて時間ステップあたりの計算量が非常に少なくて済む．例えばAB2は実際には次の形で使用される．

(

1

)

1 3 2 − + ₌ n ₊ Δ n ₋ n n _u t _F _F u これは既知量（時刻tn_{以前の量）の単純代入計算のみで未知量}_un+1 が得られることを意味し，１時間ステップの進行に繰り返し計算を必要としない．このような時間進行法は陽解法と呼ばれる． AB公式陽的RK法

(21)

21

1・1・２アダムス・モルトン公式

時間１次精度から４次精度までのアダムス・モルトン（AM）公式は以下のとおりである．１次精度の公式（AM1）は陰的オイラー法であり，２次精度の公式（AM2）は通常クランク・ニコルソン（Crank-Nicolson， CN）法と呼ばれる．AM公式はun+1_{の値を得るために}_Fn+1_=F(un+1_{)を必要と} するが，このような時間進行法は陰解法と呼ばれる．陰的オイラー法(AM1）クランク・ニコルソン法(AB2) AM3 AM4

(

n 1 n

)

_Fn 1 t u u + − ₌ ₊ Δ

(

n 1 n

_{) (}

_Fn 1 _F n

₎

2 1 t u u + − ₌ ₊ ₊ Δ

(

n 1 n

_{) (}

₅_Fn 1 ₈_Fn _Fn 1

₎

12 1 t u u + ₊ ₋ − + = − Δ

(

n 1 n

_{) (}

₉_Fn 1 ₁₉_Fn ₅_Fn 1 _Fn 2

₎

24 1 t u u + − ₌ ₊ ₊ ₋ ₋ ₊ ₋ Δ

(22)

補足：AM公式の導出（１／３）

(

_{) (}

1

₎

₀ 1 0 1 1 1 = + + + − Δ − − − + + L n n n n n F F F t u u _β _β _β とおき，時刻tn_まわりの_{テイラー展開}_，

(

)

L + Δ + Δ + Δ + = Δ − + 3 3 3 2 2 2 1 24 1 6 1 2 1 t dt F d t dt F d t dt dF F t u un n _n n n n ( ) ( ) ( ) _L, 0,1,2,_L 6 1 2 1 ₃ 3 3 2 2 2 = + Δ − Δ + Δ − = − _q _t _q dt F d t q dt F d t q dt dF F F n n n n q n から，AM公式が指定した時間精度（２次精度ならばO(Δt2_{)）で近似さ} れるように係数

β

₁,β₀, β_－1, ・・・を定める． AM公式を， L + Δ + Δ + Δ + = + 3 3 3 2 2 2 1 6 1 2 1 t dt F d t dt F d t dt dF F F n n n n n

(23)

23 補足：AM公式の導出（２／３）例としてクランク・ニコルソン法（AM２）を取り上げる． AM２の係数

β

₁

,β

₀ を定めるためテイラー表を利用する． O(Δt０_)とO(_Δ_t1_{)の列それぞれの係数の合計を０とする２つの条件から，} 係数β₁,β₀ に対する連立方程式が得られる．

(

)

0 0 1 1 1 = − − Δ − ₊ + n n n n F F t u u _β _β

(

)

0 0 2 1 6 1 2 1 1 / 0 0 1 1 1 1 1 1 2 2 2 β β β β β β − − − − − − + + + Δ − Δ Δ + − n n n n n n n F F t u u t dt F d t dt dF F

(24)

補足：AM公式の導出（３／３）係数β₁ ,β₀ に対する連立方程式 ⎪⎩ ⎪ ⎨ ⎧ = = + − − 2 1 1 1 1 0 β β β これを解きβ₁=β₀=1/2を得る．これはクランク・ニコルソン法(AM2) を与える．

(

n n

_{) (}

_Fn _Fn

₎

t u u ₌ ₊ Δ − + + 1 1 2 1 なお，テイラー表の３列目（O(Δt2_{) ）の合計はクランク・ニコルソ} ン法の打切り誤差の初項εを与える．クランク・ニコルソン法(AB2) ) ( 12 1 2 1 6 1 ₂ ₂ 2 2 2 2 2 1 t O t dt F d t dt F d n_Δ ₌ ₋ n_Δ ₌ _Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = β ε

(25)

25 補足：陰解法による計算例えばクランク・ニコルソン法（AM2）は実際には上式を変形した次の形で使用される右辺は既知量であるが左辺は未知量un+1_{を含む関数であり，上式を} 満足するun+1を特定するために通常繰り返し計算が必要となる．このような時間進行法は陰解法と呼ばれる． AM公式後退差分公式陰的RK法

(

n 1 n 1

)

n

(

n n

)

1 n _F _t _u 2 t u u t F 2 t u + − Δ + , + = + Δ ,

(26)

1・1・3 後退差分公式

時間１次精度から４次精度までの後退差分(BD)公式は以下のとおりである．１次精度の公式（BD1）は陰的オイラー法である．これら後退差分公式はギアー（Gear）法と呼ばれることもある． BD公式も陰解法である．陰的オイラー法(BD1） BD2 BD3 BD4

(

n 1 n

)

_Fn 1 t u u + − ₌ ₊ Δ

(

n 1 n n 1

)

_Fn 1 t 2 u u 4 u 3 + − + − ₌ ₊ Δ

(

n 1 n n 1 n 2

)

_Fn 1 t 6 u 2 u 9 u 18 u 11 + − + − − − ₌ ₊ Δ

(

n 1 n n 1 n 2 n 3

)

_Fn 1 t 12 u 3 u 16 u 36 u 48 u 25 + − + − − − + − ₌ ₊ Δ

(27)

27 補足：BD公式の導出（１／３）とおき，時刻tn+1_まわりの_{テイラー展開}_，から，BD公式が指定した時間精度（２次精度ならばO(Δt2_{)）で近似さ} れるように係数α_１,

α

₀,

α

_-1,

α

_-2,・・・を定める． BD公式を，

(

)

0 1 2 2 1 1 0 1 1 − = Δ + + + + ₋ − ₋ − ₊ + n n n n n F t u u u u α α α _L α ( + )Δ + [( + )Δ ] − [( + )Δ ] +L − = + + + + − 3 1 3 3 2 1 2 2 1 1 1 6 1 1 2 1 1 q t dt u d t q dt u d t q dt du u u n n n n q n ( + )Δ + [( + )Δ ] − [( + )Δ ] +_L − = + + + + 3 1 2 2 2 1 1 1 1 6 1 1 2 1 1 q t dt F d t q dt dF t q F u n n n n L , 2 , 1 , 0 = q

(28)

補足：BD公式の導出（２／３）例として２次精度後退差分公式（BD２）を取り上げる． BD２の係数

α

_１,

α

₀,

α

_-1を定めるためテイラー表を利用する． O(Δt－1_{), O(}_Δ_t０_{),およびO(}_Δ_t1_{)の列それぞれの係数の合計を０とする３} つの条件から，係数α_１,α₀,α_-1 に対する連立方程式が得られる．

(

)

0 1 1 1 0 1 1 − = Δ + + ₋ − ₊ + n n n n F t u u u α α α 0 0 1 0 6 8 2 2 / 6 1 2 1 / 0 0 0 / 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 2 1 2 2 1 1 1 − − − − + Δ − + − + Δ + Δ Δ Δ Δ + − − − − − − + + + + + n n n n n n n n F t u t u t u t dt F d t dt dF F t u α α α α α α α α α α α α

(29)

29 補足：BD公式の導出（３／３）係数α_１,α₀,α_-1 に対する連立方程式これを解きα₁=3/2，α₀=-2，α_-1=1/2 を得る．これはBD2公式を与える．なお，テイラー表の４列目（O(Δt2_{) ）の合計はBD2公式の}_{打切り誤差} の初項εを与える． BD2 ⎪ ⎩ ⎪ ⎨ ⎧ = + − = + = + + − − − 0 4 1 2 0 1 0 1 0 1 0 1 α α α α α α α

(

)

1 1 1 2 4 3 + − ₊ = Δ + − n n _n n F t u u u ( ) ) ( 3 1 6 8 1 ₂ ₂ 2 2 2 1 2 2 1 0 _t _O _t dt F d t dt F d n _Δ ₌ ₋ n _Δ ₌ _Δ + − = + + − α α ε

(30)

(31)

31

1・2 ルンゲ・クッタ法

式(1.1)を時間進行するための s 段ルンゲ・クッタ法 (Runge-Kutta method) の一般公式は次式で与えられる．

(

)

( )

(

)

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ Δ + Δ + = ≤ ≤ Δ + Δ + = ∑ ∑ = + = s i i i n i n n j j n s j ij n i u t c t F b t u u s i u t c t F a t u u 1 ) ( 1 ) ( 1 ) ( , 1 , , ただし，係数に対し

(

i s

)

a c s j ij i =

∑

≤ ≤ =1 , 1 (1.2.1) の条件を与える． (1.2.2) 係数_b_i, _c_i,および _a_ij は要求精度を満たすように式(1.1)のテイラー展開との比較から決定される． ( ) ₍ _{( )}_{) (} 0 ,u t t t t F dt ) t du ₌ _>

( )

t₀ u₀ u = (1.1) (1.2)

(32)

RK法の係数は次のブッチャー配列（Butcher tableau）で表現すると便利である． s ss s s s s s b b b a a a c a a a c a a a c L L M O M M M L L 2 1 2 1 2 22 21 2 1 12 11 1 (1.2.3) s段RK法のブッチャー配列 T A b c

(33)

33

( )

1 2 3 4 4 5 6 6 7 7

( )

2 10 10 9 8 7 6 5 4 3 2 1 − ≤ ≥ s s p s p s s なお，j≧i に対してa_ij=0とするものを(陽的)RK法と呼ぶ．s 段(陽的)RK法 の到達可能精度p(s)は次表となることが知られている． 精度と計算量の兼ね合いからは４段４次精度公式が最も効率が良いことがわかる．また，j>iに対してa_ij=0とするものを半陰的RK法と呼ぶ．ｓ段半陰的RK法の到達可能精度はp(s)=(s+1)となることが知られている． さらに，あるj>iについてa_ij≠0のものを陰的RK法と呼ぶが，ｓ段陰的RK法の到達可能精度はp(s)=2sである．

(34)

補足：係数に対する条件式(1.2.2)に対しての必然性はないが，このように設定すると扱い易い形となるので，通常この条件を加える．また，常微分方程式の数値解法の教科書におけるRK法の一般公式は式 (1.2.1)よりも次式によるものが多い．しかしここでは，計算流体力学への応用における明快性から式(1.2.1)を一般公式として採用する． ( ) ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ Δ + = ≤ ≤ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ + Δ + = ∑ ∑ = + = s i i i n n s j ij j n i n i K b t u u s i K a t u t c t F K 1 1 1 1 , ,

(35)

35

1・2・1 陽的ルンゲ・クッタ法

陽的RK法はj>iに対してa_ij=0と設定するものであり，s段陽的ＲＫ法に対するブッチャー配列は以下となる．陽的RK法のブッチャー配列 s s ss s s s b b b b a a a c a a c a c 1 2 1 1 2 1 32 31 3 21 2 0 0 0 0 0 0 0 0 − − L L M O O M M M L L L L L L

(36)

この場合（_c₁=0として）， ( ) _t _c _t

₍

_i _s

₎

t i = n + _iΔ , 1≤ ≤ ( ) _F

(

_t( ) _u( )

)

₍

_i _s

₎

F i = i , i , 1≤ ≤ と記述すると，ｓ段陽的RK法は次のように表現できる（更新方式）. ( ) ( ) ( ) ( ) ( ) ₍ ₎ ( ) ⎪ ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎪ ⎨ ⎧ Δ + = ≤ ≤ Δ + = + = = =

∑

= + − = s i i i n n i n i 1 i j j ij n i n n F b t u u s i t c t t F a Δt u u t t u u 1 1 1 1 1 2 , , , なおこの公式から，_c_iは_{i段目の計算時点でのtからt+}

Δ

_tの間 の到達割合いを意味することがわかる．

(37)

37 以下４次精度までの陽的ＲＫ法の公式の具体例を挙げる．まず，陽的オイラー法は１段１次精度陽的RK法(RK1)と見なすことができ，対応するブッチャー配列と計算式は 1 0 0 ⎪⎩ ⎪ ⎨ ⎧ ⋅ Δ + = = = +1 (1) ) 1 ( ) 1 ( , F t u u t t u u n n n n 1 11 1 b a c となる．これが陽的オイラー法 n n n _u _t _F u +1 = + Δ ⋅ と一致することは容易に理解できる．

(38)

２段２次精度陽的RK法（RK2）２段陽的RK法の一般形は次式で与えられる. ( ) ( ) ( ) ( ) ( ) ( ) ( )

(

)

⎪ ⎩ ⎪ ⎨ ⎧ + + = Δ + = ⋅ + = = = + 2 2 1 1 1 2 2 1 21 2 1 1 , , F b F b Δt u u t c t t F a Δt u u t t u u n n n n n n ここで_b₂=

α

を自由パラメータに選ぶと，２段２次精度陽的ＲＫ法(RK2)に対する係数のブッチャー配列は以下で与えられる． α α α α − = 1 0 ) 2 /( 1 ) 2 /( 1 0 0 0 0 0 0 0 2 1 21 2 b b a c RK2の一般係数(αは自由パラメータ）

(39)

39 特に，α=1の場合を修正オイラー法，α=1/2の場合をHeun法(改良オイラー法)と呼ぶ．またα=3/4の場合はRalstonの最適公式である． 1 0 0 2 / 1 2 / 1 0 0 0 2 / 1 2 / 1 0 1 1 0 0 0 4 / 3 4 / 1 0 3 / 2 3 / 2 0 0 0 修正オイラー法(α=1）改良オイラー法(α=1/2） Ralstonの最適公式(α=3/4） t u(t) t n t n+1 RK2 un+1 u(2) = u(1) un t(2) t(1) = F (1) F (2) Δt/2 u(t) (α=1, Modified Euler) Δt/2 図1.5 RK2(修正オイラー法)による時間進行

(40)

補足：RK2の係数（１／３）と書ける．ただし，F=F(t,u(t))，F_u=∂F/∂u，およびとする．これより，u(t+Δt)のテイラー展開は次式で記述される． まずF=F(t, u(t))の時間微分は， F dt du = DF dt u d 2 2 = 3₃ D2F ( )DF F_u dt u d + =

( )

( ) 2 ( ) _u u u 2 3 4 4 DF DF 3 F DF F F D F D dt u d ₌ ₊ ₊ ₊ ，，， u F t D ∂ ∂ + ∂ ∂ = 2 2 2 2 2 2 2 u F u t F 2 t D ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = 3 3 3 2 3 2 2 3 3 3 3 3 3 u F u t F u t F t D ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = ，，

(

t t

) ( )

u t tF t DF t

[

D F

( )

DF F_u

]

u + Δ = + Δ + Δ + Δ 2 + 3 2 ! 3 2

( )

( ) ( )

[

+ + +

]

+L Δ + t D F D F F_u DF F_u 3 DF DF ! 4 2 2 3 4

(41)

41 補足：RK2の係数（２／３）これをＲＫ２の時間進行の式，un+1_=un₊_Δ_t(b 1F(1)+b2F(2))，に代入すると，ＲＫ２は次に，ＲＫ２の一般公式におけるF(1)，F(2)をtn_=tおよびun_{=u(t)まわりに} テイラー展開すると次式を得る． F F(1) =

( )

(

(1)

)

21 2 ) 2 ( , u t ta F t c t F F = + Δ + Δ L + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ Δ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ Δ + = F u F a t c t F u F a t c t F 2 21 2 2 21 2 2

(

t + Δt

) ( )

= u t + Δt

(

b + b

)

F + Δt

(

b c F_t + b a FF_u

)

+_L u 2 ₂ ₂ ₂ ₂₁ 2 1 の計算を行っていることに対応する．ただしF_t=∂F/∂tである． これをu(t+Δt)のテイラー展開と比較してO(Δt2_{)まで係数を一致させ，} さらに式(1.2.2)の条件も加えると，RK2の係数に対する条件式を得る．

(42)

補足：RK2の係数（３／３）ＲＫ２の係数に対する条件式これは自由度１の連立方程式（式３つ，変数４つ）であるので，b₂=α を自由パラメータとして係数の解を表現すると c₂=a₂₁=1/(2α) ， b₁=1-α を得る． 21 2 a c = 1 2 1 + b = b 2 / 1 2 2c = b

(43)

43 ３段３次精度陽的RK法（RK3）３段陽的RK法の一般形は次式で与えられる. ３段陽的RK法に対する係数のブッチャー配列は以下となる． 3段陽的RK法のブッチャー配列 ( ) ( ) ( ) ( ) ( ) ( )

(

( ) ( )

)

( ) ( ) ( ) ( )

(

)

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + + + = Δ + = + + = Δ + = ⋅ + = = = + , , , , 3 3 2 2 1 1 1 3 3 2 32 1 31 3 2 2 1 21 2 1 1 F b F b F b Δt u u t c t t F a F a Δt u u t c t t F a Δt u u t t u u n n n n n n n n 3 2 1 32 31 3 21 2 0 0 0 0 0 0 0 b b b a a c a c

(44)

代表的な３段３次精度RK法(ＲＫ３)に対する係数のブッチャー配列を挙げておく． 古典的RK3(c₂=1/2,c₃=1） HeunのRK3(c₂=1/4,c₃=2/3） RalstonのRK3(c₂=1/2,c₃=3/4） 6 / 1 3 / 2 6 / 1 0 2 1 1 0 0 2 / 1 2 / 1 0 0 0 0 − 4 / 3 0 4 / 1 0 9 / 8 9 / 2 3 / 2 0 0 4 / 1 4 / 1 0 0 0 0 − 9 / 4 3 / 1 9 / 2 0 4 / 3 0 4 / 3 0 0 2 / 1 2 / 1 0 0 0 0 4 / 3 0 4 / 1 0 12 / 5 4 / 1 3 / 2 0 0 15 / 8 15 / 8 0 0 0 0 7 / 4 7 / 3 0 0 4 / 7 1 4 / 3 0 0 6 / 1 6 / 1 0 0 0 0 − 15 / 8 10 / 3 6 / 1 0 16 / 15 16 / 3 4 / 3 0 0 3 / 1 3 / 1 0 0 0 0 − WrayのRK3(c₂=8/15,c₃=2/3） b₁=0となるRK3(c₂=1/6,c₃=3/4） WilliamsonのRK3(c₂=1/3,c₃=3/4）これらのうち，b₁=0となるＲＫ３は更新方式における低容量型ＲＫ３で ある．つまり

b

₁

=0

によってun+1_の計算に_F(1)が必要なくなるので，更新計算において他のＲＫ３公式よりも配列を１つ減らせる．WrayのＲＫ３は積み上げ方式における低容量型ＲＫ３である．WilliamsonのＲＫ３は， Williamsonの計算法に従えば低容量（２Ｎ容量）型となるＲＫ３である．

(45)

45 補足：RK3の係数（１／３）これは自由度２の連立方程式（式６つ，変数８つ）である． RK2と同様にしてRK3に対する係数の条件式も得られる． 21 2 a c = 32 31 3 a a c = + 1 3 2 1 + b +b = b 2 / 1 3 3 2 2c + cb = b 3 / 1 2 3 3 2 2 2c + cb = b 6 / 1 2 32 3a c = b

(46)

補足：RK3の係数（２／２）を得る．また，b₃=α≠0を自由パラメータとして，以下の解が知られている． c₂とc₃を自由パラメータとして係数の解を表現すると，c₂≠c₃，c₂≠0， c₂≠2/3 に対して 2 21 c a = ( ) 2 3 2 2 2 3 31 3 2 1 3 c c c c c c a − − − × = 2 2 3 2 3 32 3 2 c c c c c a − − × = ( ) 3 2 3 2 1 6 2 3 1 c c c c b = − + − ( ₃ ₂) 2 3 2 6 2 3 c c c c b − − = ₍ ₎ 2 3 3 2 3 6 3 2 c c c c b − − = ，， _， α α α α − − 4 / 3 4 / 1 0 ) 4 /( 1 ) 4 /( 1 3 / 2 3 / 2 0 0 3 / 2 3 / 2 0 0 0 0 α α α α 4 / 3 4 / 1 0 ) 4 /( 1 ) 4 /( 1 0 0 0 3 / 2 3 / 2 0 0 0 0 − −

(47)

47 ４段４次精度陽的RK法（RK4）４段陽的RK法の一般形は次式で与えられる. ４段陽的RK法に対する係数のブッチャー配列は以下となる． 4段陽的RK法のブッチャー配列 ( ) ( ) ( ) ( ) ( ) ( )

(

( ) ( )

)

( ) ( )

(

( ) ( ) ( )

)

( ) ( ) ( ) ( ) ( )

(

)

⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ + + + + = Δ + = + + + = Δ + = + + = Δ + = ⋅ + = = = + 4 4 3 3 2 2 1 1 1 4 4 3 43 2 42 1 41 4 3 3 2 32 1 31 3 2 2 1 21 2 1 1 , , , , F b F b F b F b Δt u u t c t t F a F a F a Δt u u t c t t F a F a Δt u u t c t t F a Δt u u t t u u n n n n n n n n n n 4 3 2 1 43 42 41 4 32 31 3 21 2 0 0 0 0 0 0 0 0 0 0 0 b b b b a a a c a a c a c

(48)

４次精度を満足する４段陽的ＲＫ法の係数の一般解はいくつか知られているが，ここでは代表的な４段４次精度RK法(ＲＫ４)に対する係数のブッチャー配列を挙げておく．古典的RK4 KuttaのRK4(Kuttaの3/8公式） GillのRK4 6 / 1 3 / 1 3 / 1 6 / 1 0 1 0 0 1 0 0 2 / 1 0 2 / 1 0 0 0 2 / 1 2 / 1 0 0 0 0 0 8 / 1 8 / 3 8 / 3 8 / 1 0 1 1 1 1 0 0 1 3 / 1 3 / 2 0 0 0 3 / 1 3 / 1 0 0 0 0 0 − −

(

) (

)

(

)

(

2 2

) (

/6 2 2

)

/6 1/6 6 / 1 0 2 / 2 2 2 / 2 0 1 0 0 2 / 2 2 2 / 1 2 2 / 1 0 0 0 2 / 1 2 / 1 0 0 0 0 0 + − + − − −

(49)

49 補足：RK4の係数（１／３）ＲＫ２あるいはＲＫ３と同様にしてＲＫ４に対する係数の条件式も得られる． 21 2 a c これは自由度２の連立方程式（式１１，変数１３）である． = 32 31 3 a a c = + 43 42 41 4 a a a c = + + 1 4 3 2 1 +b +b +b = b 2 / 1 4 4 3 3 2 2c +b c +b c = b 3 / 1 2 4 4 2 3 3 2 2 2c + b c + b c = b 4 / 1 3 4 4 3 3 3 3 2 2c +b c +b c = b ( ₄₂ ₂ ₄₃ ₃) 1/6 4 2 32 3a c +b a c + a c = b

(

₄₂ ₂2 ₄₃ ₃2

)

1/12 4 2 2 32 3a c +b a c + a c = b ( ₄₂ ₂ ₄₃ ₃) ₄ 1/8 4 3 2 32 3a c c +b a c + a c c = b 24 / 1 2 32 43 4a a c = b

(50)

補足：RK4の係数（２／３）を得る．この解はc₂=1/3，c₃=2/3でKuttaのRK4を与える． まず c₂とc₃を自由パラメータとし，c₂c₃(c₂-1)(c₃-1)(c₂-c₃)≠0の場合につい てRK4の係数の解を表現すると， 1 4 = c a₂₁ = c₂

(

)

( ₂) 2 2 2 3 2 3 31 2 1 2 4 3 c c c c c c a − − − = (₍ )₎ 2 2 2 3 3 32 2 1 2c c c c c a − − = ( )

[

( )

]

( )[6 4( ) 3] 2 1 2 1 1 3 2 3 2 2 3 2 2 3 3 2 2 42 ₋ ₋ ₊ ₊ − − − + − = c c c c c c c c c c c a ( )( )( ) ( )[6 4( ) 3] 1 1 2 1 3 2 3 2 2 3 3 3 2 2 43 ₋ ₋ ₊ ₊ − − − = c c c c c c c c c c a ( ) 3 2 3 2 1 12 2 1 2 1 c c c c b = + − + ( ₃ ₂)( ₂) 2 3 2 1 12 1 2 c c c c c b − − − = ( ₃ ₂)( ₃) 3 2 3 1 12 2 1 c c c c c b − − − = ( ₄₂ ₄₃) 41 1 a a a = − + b₄ =1−(b₁ +b₂ +b₃) また，，，，，，，，

(51)

51 補足：RK4の係数（３／３）この解はα=1/3で古典的ＲＫ４，α=(2+√2)/6でGillのＲＫ４を与える．他にも４つの１パラメータ解が知られているが，そのうちの２つを示しておく（ただしα≠0）． c₂=c₃=1/2に対してはb₃=α≠0を自由パラメータとする次の１パラメータ 解が知られている． 6 / 1 3 / 2 6 / 1 0 3 3 1 0 1 0 0 ) 6 /( 1 ) 6 /( 1 2 / 1 2 / 1 0 0 0 2 / 1 2 / 1 0 0 0 0 0 α α α α α α − − − α α α α α 3 / 2 6 / 1 6 / 1 0 ) 3 /( 1 ) 12 /( 1 ) 4 /( 1 1 1 0 0 8 / 1 8 / 3 2 / 1 0 0 0 1 1 0 0 0 0 0 − − − 6 / 1 3 / 2 6 / 1 0 6 2 / 3 6 2 / 1 1 0 0 ) 12 /( 1 ) 12 /( 1 0 0 0 0 2 / 1 2 / 1 0 0 0 0 0 α α α α α α − − − −

(52)

1・2・2 陰的ルンゲ・クッタ法

陰的ＲＫ法はあるj>iについてa_ij

≠

0なので，基本的に式 (1.2.1)をそのまま適用する必要がある．例えば２段陰的ＲＫ法に対するブッチャー配列と計算式は 2 1 22 21 2 12 11 1 b b a a c a a c

(

)

(

)

{

}

(

)

(

)

{

}

⎪⎩ ⎪ ⎨ ⎧ Δ + + Δ + Δ + = Δ + + Δ + Δ + = ) 2 ( 2 22 ) 1 ( 1 21 ) 2 ( ) 2 ( 2 12 ) 1 ( 1 11 ) 1 ( , , , , u t c t F a u t c t F a t u u u t c t F a u t c t F a t u u n n n n n n

(

)

(

)

{

(2)

}

2 2 ) 1 ( 1 1 1 , , u b F t c t u t c t F b t u un+ = n + Δ n + Δ + n + Δ およびであり，u(1)とu(2)の連立方程式を解いた後にun+1を計算する．一般にs段陰的ＲＫ法では_u(1)_{, u}(2)_{, …, u}(s)の連立方程式をニュートン法等の繰り返し法で時間ステップ毎に解く必要があり計算量が膨大となる．

(53)

53 半陰的ＲＫ法は_{j>iについてa}_ij_{=0と設定するものであり，例} えば２段半陰的ＲＫ法に対するブッチャー配列と計算式は

(

)

(

)

{

(2)

}

2 2 ) 1 ( 1 1 1 , , u b F t c t u t c t F b t u un+ = n +Δ n + Δ + n + Δ およびとなる．この場合はu(1)およびu(2)に対する単独の陰的方程式を順番に解いた後に_un+1_{を計算すれば良いので，１段あたり} の計算負荷はAM法と同程度となる．陰的RK法はs段で最大2s次精度まで，半陰的RK法ではs段で最大(s+1)次精度まで到達可能である． 2 1 22 21 2 11 1 0 b b a a c a c

(

(1)

)

1 11 ) 1 ( , u t c t F a t u u = n + Δ ⋅ n + Δ

(

)

(

)

{

(2)

}

2 22 ) 1 ( 1 21 ) 2 ( , , u a F t c t u t c t F a t u u = n + Δ n + Δ + n + Δ

(54)

陰的RK法（IRK法）および半陰的RK法（SIRK法）の例を挙げておく．

陰的オイラー法(1段1次IRK法）陰的中点則(1段2次精度IRK法） CN法(2段2次SIRK法）

2段4次IRK法 2段3次SIRK法(ノルセット法） 1 1 1 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 1 0 0 0

(

)

(

)

(

) (

)

2 / 1 2 / 1 4 / 1 12 / 3 2 3 6 / 3 3 12 / 3 2 3 4 / 1 6 / 3 3 + + − −

(

) (

)

(

)

(

)

2 / 1 2 / 1 6 / 3 3 3 / 3 6 / 3 3 0 6 / 3 3 6 / 3 3 + − − + + 3段6次IRK法

(

)

(

)

18 / 5 9 / 4 18 / 5 36 / 5 15 / 15 9 / 2 30 / 15 36 / 5 10 / 15 5 24 / 15 36 / 5 9 / 2 24 / 15 36 / 5 2 / 1 30 / 15 36 / 5 15 / 15 9 / 2 36 / 5 10 / 15 5 + + + − + − − −

(55)

55 ここで，陰的オイラー法，陰的中点法，クランク・ニコルソン法（CN法）について補足しておく．陰的オイラー法は１段１次精度陰的RK法とみなせ，ブッチャー配列と計算式はそれぞれ，と記述できる．この場合の２つの計算式は等価でu(1)_=un+1であり，これを用いるとAM公式での陰的オイラー法， 1 1 + + ₌ n ₊ _Δ _⋅ n n _u _t _F u が得られる． 1 1 1

(

(1)

)

) 1 ( , u t t F t u u = n + Δ ⋅ n + Δ

(

(1)

)

1 , u t t F t u un+ = n + Δ ⋅ n +Δ および

(56)

陰的中点則は１段２次精度陰的RK法とみなせ，ブッチャー配列と計算式はそれぞれ，と記述できる．２つの計算式を比べるとu(1)_=(un+1_+un_{)/2となっており，} これが陰的中点則，であることがわかる．および 1 2 / 1 2 / 1 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ₊ Δ Δ + = (1) ) 1 ( , 2 2 u t t F t u u n n ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ₊ Δ ⋅ Δ + = +1 (1) , 2 u t t F t u un n n ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ _Δ ₊ + ⋅ Δ + = + + 2 , 2 1 1 n n n n n _u _t _F _t t u u u

(57)

57 クランク・ニコルソン法（CN法）は２段２次精度陰的RK法とみなせ，ブッチャー配列と計算式はそれぞれ，と記述できる．この場合 u(1)_=un，u(2)_=un+1であり，これがAM公式でのクランク・ニコルソン法，であることがわかる．および 1 2 / 1 2 / 1

(

n n

)

n n _u t _F _F u +1 = + Δ +1 + 2 n u u(1) =

(

(1)

)

(

(2)

)

) 2 ( , 2 , 2 F t t u t u t F t u u = n + Δ n + Δ n +Δ

(

(1)

)

(

(2)

)

1 , 2 , 2 F t t u t u t F t u un+ = n + Δ n + Δ n + Δ

(58)

(59)

59

1・3 時間進行法の安定性解析

線形テスト問題に対する時間進行法の数値安定性を考える．まずこの問題の解析解はで与えられ，

λ

_R

≦

0ならばt>0でu(t) は有界となる．

u

dt

du

₌

_λ

( )

0 u

₀

u

=

，

λ

=

λ

_R

+

i

λ

_I

( )

_t _u _e t _u _e _R t_ei _I t u = ₀ λ⋅ = ₀ λ ⋅ λ ⋅ –2 –1 0 1 2 –2 –1 0 1 2 λ R λ I Stable Unstable 図1.6 線形テスト問題の安定領域 (複素平面)

(60)

補足：なぜ線形テスト問題を考えるのか（１/２）流体運動の支配方程式は輸送方程式であるので，その本質を表現するモデル方程式（１次元の移流・拡散方程式の空間離散式）を考える． 2 2 x u x u U dt du δ δ ν δ δ ₊ − = ここで，U：移流速度 ν：拡散係数 δ/δx：∂/∂x に対する空間離散オペレータ 解析的な扱いを容易にするため，x方向に周期的な場を仮定し，離散フー リエ変換を導入する． ( )

( )

[

( )

]

( ) k L N x L k N N k x k i t x u t u _x x j N j j k = ∑ − = − − = = ⋅Δ − = , 2 1 2 / 2 / , exp , ˆ 1 0 π ω ω ～，

( )

, ˆ

( )

exp

[

( )

]

, 0 1 1 2 / 2 / − = + =

∑

− − = N j x k i t u t x u N _j N k k j ω ～

(61)

61 補足：なぜ線形テスト問題を考えるのか（２/２）ここで，ω’(k)とω”(k) は空間１階および２階微分に対する離散オペレー タの修正波数で，空間離散化手法に応じて定まる．例えば，離散フーリエ変換によりモデル方程式を対角化する． ( ) _[ _{( )} _{( )}_] _{( )} 1 2 / 2 / , ˆ ˆ − − = ′ − ′′ − = k iU k u k k N N dt t u d k k _ν_ω _ω _～ ( ) ( ) ( ) ( )2 , k k k k ω ω ω ω′ = ′′ = ( ) ( ( ) ) ( ) ( ( ) ) ( )2 2 2 / 2 / sin , sin x x k k x x k k Δ Δ = ′′ Δ Δ = ′ ω ω ω ω フーリエ・スペクトル法：２次精度中心差分：以上より，モデル方程式は線形テスト問題（固有値問題）に帰着され，固有値の実部は拡散項に，固有値の虚部は移流項に関係することが理解される． ( ) _{( )} _{( )} _{( )} k U k i t u dt t u d k I k R k I k R k k k k ₌ _λ _ˆ _, _λ ₌ _λ ₊ _λ _, _λ ₌ ₋_ν _ω_′′ _, _λ ₌ ₋ _ω_′ ˆ なお，純粋移流問題の安定性はω’(k)の最大値，純粋拡散問題の安定性 はω”(k)の最大値に対応して定まる．

(62)

線形テスト問題の解の有界性が数値計算においてどのような条件で満足されるかを検討するため，まず陽的オイラー法を線形テスト問題に適用してみる．

( )

n

(

)

n n n _u _t _u _t _u u +1 = + Δ ⋅ λ = 1+ λΔ この操作は，初期値u0_=u 0，数値解の複素増幅率をξとして 0 u un =

ξ

n と記述するとき，陽的オイラー法の複素増幅率が

(

+ z

)

= 1 ξ ，ただし z = λΔt となることを意味する．ここでz=λΔt=λ_RΔｔ+iλ_IΔtとしている．

(63)

63 これから明らかなとおり，陽的オイラー法による数値解が有界であるための条件は 1 1+ ≤ = z

ξ

で与えられる．これを満足するzは複素 平面において(-1,0i)を中心とする半径１の円の内側に存在する．この領域を陽的オイラー法の絶対安定領域と呼ぶ．これは，陽的オイラー法を用いて安定な数値解を得るためには，時間刻み幅の設定に制限があることを示している．具体的には，λを負の実数とすると λ − ≤ Δt 2 の範囲に時間刻み幅Δtを設定する必要 がある． –3 –2 –1 0 1 –2 –1 0 1 2 λ R_Δ t λ I Δ t

Euler explicit

Unstable Stable 図1.7 陽的オイラー法の絶対安定領域

(64)

次に，陰的オイラー法を線形テスト問題に適用し絶対安定領域を調べる．

( )

1 1 + + ₌ n ₊ _Δ _⋅ n n _u _t _u u λ これより，陰的オイラー法の複素増幅率ξは次式で与えられる．，ただし z = λΔt ( − z) = 1 1 ξ これに対応する絶対安定領域は次式で与えられる． 1 1 1 ≤ − = z ξ これを満足するzは複素平面において (-1,0i)を中心とする半径１の円の外側の全領域である．この領域は λ_R≦0となる複素平面の左半分を全 て含むので，陰的オイラー法は λ_R≦0に対し時間刻み幅Δtによらず 無条件安定な数値解を与える． –1 0 1 2 3 –2 –1 0 1 2 λ RΔ t λ I Δ t

Euler implicit

Unstable Stable 図1.8 陰的オイラー法の絶対安定領域

(65)

65 なお，絶対安定の緩和表現として，絶対安定領域が複素平面上で無限扇形領域 {z; -α < (π-arg z) < α} (0<α<π/2) を含む場合，A(α) 安定と呼ぶ． Re(z) α α O Im(z) Stable Unstable 陰的オイラー法のように絶対安定領域が複素平面の左半分(λ_R≦0)全 てを含む時間進行法はＡ安定（絶対安定）という．図1.9 A(α)安定 ところで，陽的オイラー法および陰的オイラー法は１次精度の時間進行法であるため打ち切り誤差の影響が懸念される

(66)

そこで２次精度の時間進行法である中点蛙跳び法（leap-frog法）を線形テスト問題に適用してみる．これに un₌_ξn_u0 _および _z=_λΔ_{t を代入して整理すると，}_{中点蛙跳び法}_の複素増幅率ξはの解として次式で得られる．ここでは複素増幅率ξの解を求めるのではなく，絶対安定領域の境界|ξ|=1を与えるｚの範囲を調べてみる．そのため上式を変形しξ=e iθ_{を代入すると，}

( )

n n n _u _t _u u +1 = −1 + 2Δ ⋅ λ 0 1 2 2 ₋ ξ ₋ ₌ ξ z

( )

_θ ξ ξ θ θ θ θ cos 2 2 1 2 1 2 2 i i e e i e e z i i i i = − = − = − = − を得る．これは，zが虚軸上の(-i, +i)の範囲のみで安定，つまり中点蛙跳び法はλ_R=0以外では数値解が無条件に不安定となることを示している．以上のように，時間進行法は必ずしも精度が高いものが適切なわけではなく，数値安定性も含めて検討する必要がある．

(67)

67

1・3・1 線形多段階法の絶対安定領域

線形多段階法の一般形（ただし

α

₁≠0），を線形テスト問題に適用し，_un₌

_ξ

n_u0_および

_z=λΔt

_を代入すると次式を得る．この式による_{複素増幅率の根が全て |}

ξ

|≦1 を満足する場合 に線形多段階法は絶対安定である．

∑

− = =− + + ₌ 1 N j 1 N j j n j j n ju F t 1 _α _β Δ ∑ ∑ − = =− + + ₋ ₌ 1 1 0 N j j N j n j j n jξ z β ξ α

(68)

線形多段階法の安定限界曲線（絶対安定領域の境界），つまり |

ξ

|=1 を与える z は，LM法の安定限界関数，

( )

に対して

ξ

_=eiθ _{を代入し，０から２πの範囲の}

_θ

_に対する_{zを描けば良い．} 以下に，AM公式，AM公式およびBD公式に対する安定限界関数と安定限界曲線を示しておく． j n j N j N j j n j z + − = − = + ∑ ∑ = = ξ β ξ α ξ ϕ 1 1 LM法の安定限界関数

(69)

69 AB公式陽的オイラー法(AB1) 1 − = ξ z

(

)

1 3 1 2 − − = ξ ξ ξ z

(

)

5 16 23 1 12 2 2 + − − =

ξ

z

(

)

9 37 59 55 1 24 2 3 3 − + − − = ξ ξ ξ ξ ξ z AB2 AB3 AB4 –3 –2 –1 0 1 –2 –1 0 1 2 λ R_Δ t λ I Δ t AB3 AB4 AB2 Euler explicit us u s u s us u u 図1.10 AB法の安定限界曲線

(70)

AM公式陰的オイラー法(AM1) クランク・ニコルソン法(AM2) AM3 AM4 ξ ξ −1 = z

(

)

1 1 2 + − = ξ ξ z

(

)

1 8 5 1 12 2 ₊ ₋ − = ξ ξ ξ ξ z

(

)

1 5 19 9 1 24 2 3 2 + − + − = ξ ξ ξ ξ ξ z –6 –4 –2 0 2 –4 –2 0 2 4 λ R_Δ_t λ I Δ t _AM3 AM4 AM2 AM1 Euler explicit us u s us u s u s 図1.11 AM法の安定限界曲線

(71)

71 BD公式陰的オイラー法(BD1) BD2 BD3 BD4 ξ ξ −1 = z 2 2 2 1 4 3

ξ

− + = z 3 2 3 6 2 9 18 11 ξ ξ ξ ξ − + − = z 4 2 3 4 12 3 16 36 48 25 ξ ξ ξ ξ ξ − + − + = z –4 0 4 8 12 –8 –4 0 4 8 λ R_Δ_t λ I Δ t BD3 BD4 BD2 BD1 Euler explicit u s u s us u s us 図1.12 BD法の安定限界曲線 –6 –3 0 3 6 –6 –3 0 3 6 λ R_{Δ t} λ I Δ t BD3 BD4 u s u s 73o 86o 図1.13 BD3とBD4の原点付近の安定限界曲線

(72)

以上より，線形多段解法（AB法，AM法，BD法）では，精度が上がると絶対安定領域が狭くなることがわかる．また，これらの中でA安定な時間進行法は，陰的オイラー法（AM1），クランク・ニコルソン法（AM2），および２次精度後退差分法（BD2）のみである．なお，BD3はA(86°_)安定_，_BD4_は_{A(73°)安定}_である．

将来行ってみたい研究

流体数値計算法

講師：名古屋工業大学・大学院工学研究科

教授 森西洋平

( )

t

u

u

=

( )

(

( )

) (

)

,

u

t

t

t

t

F

dt

t

du

=

>

Δ

Δ

Δ

Δ

( )

( )

(

)

(

)

Δ

(

)

1・1 線形多段階法

α

≠

α

α

=-1

β

α

≦-

＝０

α

α

-

β

≠

α

≦-

＝０

α

≠

β

β

≦０

＝０

(

) (

)

β

(

) (

)

β

≠

(

)

β

≠

1・1・1 アダムス・バシュフォース公式

(

)

(

教授森西洋平

₍

_{( )}

_{) (}

₎

₌

_>

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎

1・1・２アダムス・モルトン公式

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎

_{) (}

₎