幾何的マルチグリッド法の導出

離散化された(C.8)^式をFig. C.1^左の1^〜4^{の点で求め},^{足し合わせる}.

∑4 i=1

(∇ · ®u^∗_i)

−

∑4 i=1

(A_ieδp_ie+A_iwδp_iw+ A_inδp_in+A_isδp_is+ A_icδp_ic)=0 (C.9) この式を0とするように各差圧の修正量∆δp^を加える.

∑4 i=1

(∇ · ®u^∗_i)

−

∑4 i=1

©­

«

Aie(δpie+∆δpie)+Aiw(δpiw+∆δpiw) +A_in(δp_in+∆δp_in)+A_is(δp_is+∆δp_is)

+Aic(δpic+∆δpic) ª®

¬

=0 (C.10)

さらに,∆付きの部分を分離すると

∑4 i=1

(∇ · ®u^∗_i)

−

∑4 i=1

©­

«

A_ieδp_ie+ A_iwδp_iw +Ainδpin+Aisδpis

+A_icδp_ic ª®

¬

−

∑4 i=1

©­

«

A_ie∆δp_ie+A_iw∆δp_iw +Ain∆δpin+ Ais∆δpis

+A_ic∆δp_ic

ª®

¬

=0 (C.11)

(C.11)^式の前2^項はn+1ステップでの発散値の和と等価なので,

∑4 i=1

(®∇ · ®u_iⁿ⁺¹ )

とおく.

∑4 i=1

(®∇ · ®u_iⁿ⁺¹ )−

∑4 i=1

©­

«

Aie∆δpie+Aiw∆δpiw

+A_in∆δp_in+A_is∆δp_is +Aic∆δpic

ª®

¬

=0 (C.12)

と書き表すことができる. ^ここで,最も簡単なものとして以下のような近似を施す.

∆δp_ic =δP_1C(i=1,· · · ,4) (C.13)

Fig. C.1: Approximation for four grids

(C.13)^{式を用いて}, (C.12)^{式を書き替えると},

∑4 i=1

(®∇ · ®u_iⁿ⁺¹ )−



(A3e+A4e)δPE+(A1w+A2w)δPW

+(A_2n+A_4n)δP_N +(A_1s+ A_3s)δP_S

− {(A_3e+A_4e)+(A_1w+ A_2w)+(A_2n+A_4n)+(A_1s+A_3s)}δP_C





=0 (C.14) となる. ^ここで,^{粗格子での各差圧}δP^{の係数はそれぞれ} A_1E,A_1W,A_1N,A_1S,A_1C ^{と置くことができ}

るので, (C.14)式は以下のように表すことができる.

∑4 i=1

(®∇ · ®uⁿ⁺¹_i

)− {A_1EδP_E + A_1WδP_W +A_1NδP_N +A_1SδP_S− A_1CδP_C}=0 (C.15)

発散値をR^∗_j ^として, (C.15)^式は

R^∗_j− {A1EδPE +A1WδPW + A1NδPN + A1SδPS−A1CδPC}=0 (C.16)

となる. (C.8)^式と(C.16)^{式を比較してみると},同じ形の方程式となっていることが分かる. ^このよう

に密格子と粗格子を補間することで,プログラムコードを作成する際に同じアルゴリズムで解くこと が可能となり,コーディングが比較的容易となることがメリットとなる.

C.4 2 段グリッド法

ここでは,マルチグリッド法の基礎となる2段グリッド法の概要について述べる. 2^{段グリッド法で} は計算格子より粗い格子を1^つ生成し,^この2つの格子を交互に解くことで効率的に解を求めること ができる. 細かい格子（計算格子）での方程式には下付き添え字 f ^を用い,粗い格子での方程式を下 付き添え字c^として,^以下に2段グリッド法の概要を示す.

(1) ^{細かい格子の方程式} A_fx®_f =b®_f に定常反復計算を数回適用する. (2) ^{細かい格子の残差}r®f =b®f −Afx®f を求める.

(3) 細かい格子から粗い格子への補間（制限補間）r®c = R®rf を行う. (R^{：制限補間行列}) (4) ®bc =r®c,x®c =0®として粗い格子の方程式 Acx®c =b®cに定常反復計算を数回適用する. (5) 粗い格子から細かい格子への補間（延長補間）e®_f = Pe®_c^を行う. (P^{：延長補間行列}) (6) 細かい格子の解を修正する. x®_f =x®_f +e®_f

(7) ^{細かい格子の方程式} Afx®f =b®f に定常反復計算を数回適用する.

2^{段グリッド}(^{マルチグリッド})法で重要となるのは粗い格子と細かい格子の補間である,^制限補間 と延長補間の扱いである. この補間方法によって粗い格子で解くべき方程式が決まる.

細かい格子の解の真値を x®^t_f,^近似値をx®f,^誤差をe®f とすると,

®

x^t_f = x®_f +e®_f (C.17)

であり,^{このとき残差は}

®

r_f = ®b_f − A_fx®_f

= ®b_f − A_f

(x®^t_f − ®e_f

) (C.18)

®

x^t_f ^{は方程式の真値なので}0^となり,残差は以下のように表すことができる.

®

r_f = A_fe®_f (C.19)

また,制限補間と延長補間の関係式は以下のとおりである.

® r_c = Rr®_f

®

e_f =Pe®_c (C.20)

ここで, R^とPは制限補間行列と延長補間行列を表す. (C.19)^式と(C.20)式を粗い格子での方程式 Acx®c = ®bcに代入して整理すると,以下の関係式を導くことができる.

A_c = RA_fP (C.21)

このことから,粗い格子での係数行列は細かい格子での係数行列と制限補間行列,^{延長補間行列から} 生成されることが分かった. ^また,圧力のポアソン方程式では残差がn+1ステップでの速度の発散値 を表すため,そのまま方程式の収束判定に用いることができる. Fig. C.2^{のような格子の場合},^制限補 間は以下の式により補間を行う.

r_c = 1 4 (

r¹_f +r²_f +r_f³+r⁴_f )

(C.22) また,延長補間は以下に従って補間を行う.

Fig. C.2:Restriction and prolongation

eⁱ_f = 1

2ec (i=1· · ·4) (C.23)

(C.22)^式と(C.23)^{式中に現れる}1/4^と1/2は制限補間行列と延長補間行列の係数である. ^これらの

数字は数学的には必要のないものであるが,^発散値は"単位体積当たりの流出量"^{として定義される}.^そ のため, 2^{次元計算では}1/4, 3^{次元計算では}1/8^{を掛ける必要がある}.

また,速度の修正には以下の式,

®

uⁿ⁺¹=u®ⁿ+ ∆t

ρ ®∇ (δp) (C.24)

という圧力勾配を用いて行っているため, "粗い格子での圧力修正量と細かい格子での圧力修正量は同 じである"^{必要があるので}1/2^を掛ける.

∂p

∂x

f

= pf1−pf2

∆x

∂p

∂x

c

= p_c1−p_c2 2∆x

 →p_f1 = 1 2p_c1 p_f2 = 1

2p_c2

(C.25)

Fig. C.3: Relationship between coarse and fine grid widths

(C.22) ^式と(C.23) ^{式を用いると}, 粗い格子での係数行列も大規模な対角疎行列となり,^その成分 A_cc,A_ce,A_cw, A_cn, A_csは以下のように細かい格子での係数行列により求まる.

A_cci,j =−(

A_cei,j+ A_cwi,j+ A_cni,j+A_csi,j) Acei,j = 1

8

(Af e2i,2j−1+Af e2i,2j

)

Acwi,j = 1 8

(Af w2i−1,2j−1+Af w2i−1,2j

)

Acni,j = 1 8

(Af n2i−1,2j+Af n2i,2j

)

Acsi,j = 1 8

(Af s2i−1,2j−1+Af s2i,2j−1

)

(C.26)

以下,4×4格子における具体例を元に, (C.26)式が成り立つことを示す. ^格子のindex^を1^次元配 列0^〜15^とすると,^計算格子4×4の係数行列は以下のとおりである.

Fig. C.4: Index of4×4grids

また,制限補間行列と延長補間行列は以下のようになる.

これらを用いることで粗い格子での係数行列が以下のように求まる. (C.7)^{式を考慮して},^すべてのAc に対して

Aci =−(Aei+Awi+Ani+Asi) (i=0^〜15) (C.27)

を代入して整理すると次の式を得ることができる.

これにより粗い格子での係数行列を導出することができた. 細かい格子での係数行列と比較をして

みると, 1/8^{という係数は違うが},成分は同じ形式をとっていることが分かる.

マルチグリッド法ではこの2段グリッド法を再帰的に用いることにより, 2^{つの格子を繰り返して} 解くのではなく,様々な格子サイズを解く. ^{再帰の方法によって}V^サイクル, W^{サイクルなどがある} が,^ここではVサイクルについてのみ述べることとする.^以下のFig. C.5^とFig. C.6^{にマルチグリッ} ド法の概要を示す. Fig. C.5^でのsmootherとは定常反復計算を数回適用する作業のことである.

Fig. C.5:Schematic of Multigrid method (V-Cycle)

Fig. C.6:Restriction and prolongation of Multigrid method

参考文献

[1] ^姫野武洋: "ロケット・宇宙機における推進薬管理",^{日本混相流学会誌}(^混相流), Vol.27, No.4, pp.385-392 (2013).

[2] ^姫野武洋: "液体ロケットの推進薬挙動予測",日本流体力学学会誌「ながれ」, Vol.32, No.3 (2013).

[3] ^前村孝志,^後藤智彦,^秋山勝彦,^二村幸基,^{渡邉篤太郎}: "H-IIAロケットの新技術と初号機打ち上 げ結果",^{三菱重工技報}, Vol.39, No.1 (2002-1).

[4] Gen, M., Ogasawara, K., Kitayama, O., Ochiai, T.: "Conceptual Study of Space Environmental Experiment Platform using H-IIA 2nd Stage", 2008-g-08, ISTS 2008 (2008).

[5] Kitayama, O., Tokunaga, T., Igarashi, I., Okita, K., Fujita, M.: "CRYOGENIC PROPELLANT MANAGEMENT OF H-IIA UPPER STAGE", 2006-a-08, ISTS 2006 (2006).

[6] Himeno, T., Watanabe, T., Nonaka, S., Naruo, Y., Inatani, Y. and Aoki, H.: "Numerical and Experimental Investigation on Sloshing in Rocket Tanks with Damping Devices", AIAA 2007-5557.

[7] ^大橋昭文: "バッフル板が極低温流体のスロッシングと圧力変動に及ぼす影響",^{東京大学学士論}

文(2015).

[8] ^{古市侑太郎},^大橋昭文,^亀山頌互,^幅大地,^{佐久間康典},^鵜沢聖治,^姫野武洋,^渡辺紀徳: "^液温成層 厚みが極低温流体のスロッシングと圧力変動に及ぼす影響",^第57回航空原動機・宇宙推進講演 会, 1B15 (2017).

[9] NASA@SC15: "Predicting Baffled Propellant Tank Slosh for Spacecraft" (Last visited Jan. 5, 2020).

https://www.nas.nasa.gov/SC15/demos/demo12.html

[10] ^姫野武洋,^渡辺紀徳: "微小重力環境における気液界面挙動の数値解析",^{日本機械学会集}(B^編), 65^巻, 635^号, pp.2333-2340 (1999).

[11] Yokoi, K.: "Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm",Journal of Computational Physics, Vol.226 , pp.1985-2002 (2007).

[12] Yokoi, K.: "A practice numerical framework for free surface flows based on CLSVOF method, multi moment meth od and density scaled CSF model: Numerical simulation of droplet splashing", Journal of Computational Physics, Vol.232, pp.252-271 (2013).

[13] Yokoi, K.: "A density-scaled continuum surface force model within a balanced force formulation", Journal of Computational Physics, Vol.278, pp.221-228 (2014).

[14] Brackbill, J.U., Kothe, D.B., and Zemach, C.: "A Continuum Method for Modeling Surface Tension",Journal of Computational Physics, Vol.100, pp.335-354 (1992).

[15] ^{關 輝彦}: "A−ϕ法による電磁浮遊技術を利用した表面張力測定法の数値解析",^{首都大学東京大}

学院修士論文(2017).

[16] 桜庭雅明，弘崎聡，樫山和男: "自由表面流れ解析のためのCIVA/VOF法に基づく高精度界面捕

ドキュメント内 2020 年 1 月 24 日 (ページ 47-54)