第 1４回数値流体力学シンポジウム E03-2 直交格子を使用した３次元の任意形状物体まわりの流体シミュレーション Computation of the Flow Field around Arbitrary Three-Dimensional Body Geometry Using Cartesian Grid 市川 治, 日本 IBM(株), 〒242-8502 大和市下鶴間 1623-14, E-mail: [email protected] 藤井 孝藏, 宇宙科学研究所, 〒229-8510 相模原市由野台 3-1-1, E-mail: [email protected] Osamu ICHIKAWA, IBM Japan Ltd. 1623-14 Shimotsuruma, Yamato, Kanagawa, Japan Kozo FUJII, The Institute of Space and Astronautical Science, 3-1-1 Yoshinodaii, Sagamihara, Kanagawa, Japan A finite difference method for the simulation using Cartesian grid is developed. As the body boundary is not necessarily located on the grid points, finite difference formulation near the boundary is critically important. The distance from the adjacent stencils to the body boundary is integrated into the finite difference formulation to satisfy the body boundary conditions. From the stability analysis, it is proved that the mid-term in the dissipation term should be treated implicitly to avoid the severe diffusion number condition. The artificial dissipation term based on the reconstruction that uses a virtual center point gives the most accurate solution near the body boundary. As for the pressure boundary condition, the internal pressure is defined by the least squares method to satisfy all the differential calculus along the normal vector to the body boundary should be zero as much as possible around the stencil. 1. β= d / Δx ν= 粘性係数 φ= 流速を想定した関数 Φ= 圧力の境界条件の２乗残差合計 ψ= セル内境界プロファイルの２乗残差合計 Ψ= セル内圧力プロファイルの２乗残差合計 ω= 境界点隣接セル番号 ai = 圧力プロファイルを与える係数 ki = 境界プロファイルを与える係数 H H n = 法線ベクトル u = 流速ベクトル 序論 任意形状の物体まわりの流れを単純な直交格子を使用して 解こうという試みが近年いくつか提案されている。物体境界 の処理方法としては、セルが物体境界によって横切られる全 てのケースを評価する方法や、物体境界を格子に沿った階段 状の境界として近似する方法などがある。前者ではカットセ ルによる方法(1)や仮想境界(2)による方法などにより物体境界 が厳密に扱われるが、その分計算処理が複雑である。後者で は計算コストは低いが、階段状近似に由来する精度不足が懸 念される。 直交格子を使用した計算では、その精度がしばしば議論の 対象となる。現在のところ、境界層の遷移まで正確に計算す るまでには至っていない。それでもなお、直交格子を使用す るメリットは、境界適合型格子を生成するコストがかからな いという点にあり、自動車のエンジン・ルーム内の流れ(3)(4) や電子機器の内部流などの実用問題で有用となる。 格子生成が不要であっても、物体境界の処理に多大な計算 コストがかかってしまっては、メリットが生きてこない。そ こで、ある程度の精度と安定性を有しながら境界処理にコス トがあまりかからない計算スキームを考える。 ここで紹介する方法は、格子点間に位置する物体境界まで の距離のみを差分スキームに取り込む。カットセルのような 有限体積評価をしないので、処理が低コストである。安定性 については、対流項、拡散項のスキームを注意深く選択する ことによって、それが改善されることが示される。圧力の境 界条件については、３次元における物体法線方向の傾きゼロ の条件のシンプルな与え方を提案する。 なお、ここではスキームの単純化のために、格子間隔以下 で入り込んでくる薄い物体のケースを考えない。 H H H a , b , c , M ： 圧力境界条件の設定で使用 1, 0, C, B, i ： 添え字または点名で使用 θ,V, G, D, A, B ： Von Neumann の安定解析で使用 3. 差分スキーム 3.1. 境界までの距離を含む差分スキーム u u-1 Body −1 u1 Body uB -1 u0 uC u1 uC 0 B e C +1 -1 f Δx (e<0, f=(Δx+e)/2 ) 0 C B e ll +1 x f Δx ll Fig. 1 Right side boundary 2. 記号 ここでは以下の記号を使用する。 c = 移流速度 d, e = 隣接セルから境界までの距離 u = 流速 f = 隣接セルから仮想中点までの距離 P = 圧力 Rc = セル・レイノルズ数 t = 時間 Re = レイノルズ数 x, y, z = 座標 Δx, Δy, Δz = 格子間隔 I= u0 ここでは、物体右側の境界における x 方向の微分を考える。 Fig. 1 に示すように境界上の点 B と 点(+1)との間に仮想的 な中点 C を考える。点 C での値 uc は、それを挟む２つの格 子点 u1 ,u0 の線形内挿により推定する。点(0)から物体境界 までの距離を e (ただし e < 0)とすると 点(0)から点 C まで の距離 f は、f = (Δx + e) / 2 で表される。 3.1.1. α= 人工粘性の係数 ２階の差分スキーム まず、テイラー展開の定義に忠実に 点 B、点(0)、点(+1) の 1 Copyright © 2000 by JSCFD é ∂ 2 φ ù ∆x − d é ∂ 2 φ ù æ ∆x + d ö é ∂φ ù c⋅ê ú +c⋅ê 2 ú⋅ −α ⋅ c ⋅ ê 2 ú ⋅ç ÷ ∂x û c 2 ∂x û ∂x û è 2 ø ë ë ë " " ! "" " """ ! """" """" ! ３点を使用して２階微分を求めると次のスキームを得る。 æ ∂ 2u ö (e ⋅ u1 + (∆x − e ) ⋅ u 0 − ∆x ⋅ u B ) çç 2 ÷÷ ≅ 2 ⋅ ∆x ⋅ e ⋅ (∆x − e ) è ∂x ø 0 --- (1) １次精度部分 æ ∂ 2u ö u − 2 ⋅ uc + u B çç 2 ÷÷ ≅ 4 ⋅ 1 (∆x − e)2 è ∂x ø c f ⋅ u1 + (∆x − f ) ⋅ u 0 ただし、 u B = 0 , uc = ∆x --- (3) １階の差分スキーム 点(0)での１階の微分項については、次のように中点 C での 微分項に、点(0)が中点からずれていることへの補正項を加え たものと理解することができる。 2 æ ∂u ö æ ∂u ö æ ∂ u ö ç ÷ = ç ÷ + çç 2 ÷÷ ⋅ (− f ) x ø è ∂x ø 0 è ∂x ø c è∂" " "" ! --- (4) 中心補正項 ここでは、 æç ∂u ö÷ ≅ u1 − u B ∆x − e è ∂x ø c --- (5) と差分化される。 これは、中心補正項を加えない時に１次精度、中心補正項 として式(1)を使用した時に２次精度、式(2)を使用した時に 疑似２次精度の差分となる。 また、対流項の場合、３点しかとれない境界部では、差分 化の際、次のように２次の人工粘性項が追加される。 （以下、 [ ] は差分による演算を表わす。 ） æ ∂ 2φ ö φ − 2 ⋅ φ c + φ −1 çç 2 ÷÷ ≅ 4 ⋅ B = 4⋅ (∆x + d )2 è ∂x ø c 中心補正項 人工粘性項 --- (6) ただし、境界の無い場所では５点を使用し、４次の人工粘 性による３次精度の風上差分を構成する。 3.2. 対流項のスキームの選択 ここでは、対流項について、式(6)の中心補正項と人工粘性項 にどの精度のスキームを使用すべきかを検討する。 １次元定常移流拡散方程式の安定性解析の手法を、以下の ように格子間に境界がある場合に拡張し、それぞれの精度の スキームを使用して適性を検査する。 Fig. 2 のように物体の左側境界が点(+1)と点(0)の間にある 場合、Φに関する１次元の定常移流拡散方程式 c⋅ ∂ 2φ ∂φ =ν⋅ 2 ∂x ∂x --- (7) の解析解は、 φ 0 − φ −1 exp( Rc) − 1 = φ B − φ −1 exp( Rc ⋅ (1 + d )) − 1 ∆x --- (8) Case Ⅰ Ⅱ Ⅲ Ⅳ Ⅴ Ⅵ と与えられる。ここで Rc はセル・レイノルズ数 Rc = c ⋅ ∆x ν --- (10) æ φ − φ 0 φ 0 − φ −1 ö d ⋅ç B − ÷ ∆x ø è d (∆x + d )2 --- (11) 検証のため、ケースⅡとⅢ、およびケースⅤとⅥを比較す ると、人工粘性として疑似２次精度のスキームを使用した方 が、２次精度スキームを使用するよりも、特に d が小さい 領域で解析解の傾向をより良く表わしていることがわかる。 中心補正項については、ケースⅣの結果より、２次精度ス キームは適していないことがわかる。 以上より、人工粘性項については疑似２次精度スキームを、 中心補正項については疑似２次精度スキームまたは補正項 なしを選択すべきであることがわかった。これはケースⅢま たはケースⅤである。今回は、解析解との一致の度合いがよ り良いケースⅤを採用した。ケースⅤは、グラフを見る限り においては、セル・レイノルズ数が大きい時に中程度の d の 領域でオーバーシュートする可能性を含んでいる。しかし、 境界部では流速が低くセル・レイノルズ数はあまり大きくは ならないこと、また、多少のオーバーシュートは連続の式を 満たすように決まる圧力のポアソン方程式により抑え込ま れることが期待できるため、大丈夫と判断した。実際にこれ を用いて計算させると、かなり粗い格子でもオーバーシュー トすることなく、きれいな解が得られている。 é ∂ 2u ù é ∂ 2 u ù æ ∆x − e ö é ∂u ù æ ∂u ö c ⋅ç ÷ ≅ c ⋅ ê ú − c ⋅ ê 2 ú ⋅ f −α ⋅ c ⋅ ê 2 ú ⋅ç ÷ ∂x û c ∂x û è 2 ø è ∂x ø 0 ë ∂x "! û "ë " ! " ""ë" " """" ! １次精度部分 人工粘性項 のように、対流項を分解し、差分化に関して Table 1 に示 すケースⅠ∼Ⅵを考える。どのケースも右辺の拡散項には、 ２次精度のスキームを使用する。 c < 0 の場合には中心補正項がスキームを風上化し、安定性 が増大してしまうと予想されるので、より厳しいケースであ る c > 0 の場合のみを考える。人工粘性の係数αを 0.5 にと り、ケースⅠ∼Ⅵについて差分式によって得られるφ0 の値 をプロットしたものを Fig. 4 に示す。 ケースＩの結果は、中程度以上の d の領域で、セル・レイ ノルズ数が２を越える場合に、φ0 が非物理的な値をとりう ることを示している。これは、この領域では何らかの人工粘 性が必須であることを表わしている。一方で注目すべきは、 d が小さい（境界までの距離が近い）領域では、セル・レイ ノルズ数にかかわらず、人工粘性なしでも解析解に近い妥当 な値を示していることである。これは、人工粘性としては d が小さい時には効果が小さくなるようなものを選択すべき であることを示唆している。ところで、疑似２次精度の２階 の差分スキームは、次のように d によりスケーリングされ ていると見ることができるので、d が小さい時には効果が自 動的に小さくなることが期待できる。したがって、この差分 スキームは、今回の目的に適していると予想される。 --- (2) これは、点(0)での２階微分として使用した場合、 e ≈ −∆x の時に２次精度を有するが、その他の時は内挿によって精度 が落ちるので、ここでは疑似２次精度のスキームと呼ぶこと にする。 3.1.2. 中心補正項 é ∂ 2φ ù =ν⋅ ê 2 ú ë ∂x û これは２次精度のスキームである。 一方で、文献 5 にあるように、仮想的な中点 C、点 B、点(+1) の３点を使用して２階微分を求めると次のスキームを得る。 --- (9) である。φB=0, φ-1= 1 として、この解析解を、点(0)から物 体境界点 B までの距離 d を横軸に、その時のφ0 の値を縦 軸にプロットしたものを Fig. 3 に示す。 次に、この定常移流拡散方程式を差分化する。 Table 1 Cases of Scheme Selection Center Corrective Term Artificial Dissipation Term None None None 2nd order None Pseudo 2nd order nd 2 order 2nd order nd Pseudo 2 order Pseudo 2nd order nd Pseudo 2 order 2nd order 2 Copyright © 2000 by JSCFD ∂φ ∂ 2φ = ν⋅ 2 ∂t ∂x φ φ-1 φB (C) (-1) 0 B d f Δx ll ２次精度の差分スキームを使用し、時間方向には陽解法を 適用すると、次の差分式を得る。 Body φ0 φ in +1 − φ in 2 −ν⋅ {φin_ d − (1 + β ) ⋅ φ in + β ⋅ φin−1 }= 0 ∆t d ⋅ (d + ∆x ) --- (13) x (+1) ここでは、 β = d ∆x ll である。 φ in_ d = V n ⋅ e I ⋅(i⋅θ +θ ⋅β ) φ Analytical Solution 0.8 0.6 0.4 n =V ⋅e --- (16) I ⋅(i ⋅θ −θ ) --- (17) V n +1 = G ⋅ V n ここで、 G = 1 − D ⋅ 1 + β − e I ⋅i⋅ β − β ⋅ e − I ⋅θ 2 ⋅ ∆t D =ν⋅ d ⋅ (d + ∆x) ( Rc=2 Rc=4 0.2 0 n i −1 ここでは、I は純虚数である。 これらを式(13)に代入すると次式を得る。 Φ Rc=6 0 --- (14) この差分式について Von Neumann の安定性解析を行う。 --- (15) φ in = V n ⋅ e I ⋅i⋅θ Fig. 2 Left side boundary 1.2 1 --- (12) 0.5 1 d G 2 = (1 − D ⋅ A) − (D ⋅ B ) 2 Φ Case I 1 1 0.8 0.8 0.6 0.6 Rc=2 Rc=4 Rc=6 0.4 0.2 0 Φ 1.2 0.5 d Case III ( Φ 1.2 1 0.8 0.6 0.5 0.2 Rc=2 Rc=4 Rc=6 0.4 d 0 0 0 Φ 1.2 0.5 1 Case V Rc=2 Rc=4 Rc=6 0 Φ 1.2 1 1 0.8 0.8 0.5 d Rc=2 Rc=4 Rc=6 0.4 0.2 0 D≤ Case VI 0 0.5 ν⋅ Rc=2 Rc=4 Rc=6 0.2 1 0 0 0.5 --- (28) 1 ∆t ≤ d ⋅ ∆x 2 --- (29) 通常知られている拡散数制限は、 d 1 ν⋅ Fig. 4 CaseⅠ∼Ⅵ 3.3. 1 1+ β であれば十分である。 式(20)と式(28)より、次式が得られる。 0.4 d 2 しい。よって、 1 0.6 0.6 --- (25) --- (26) 0 ≤ B ≤ (1 + β ) 2 2 であるので、 A ≤ B の時、式(24)は常に成立する。 A 2 ≥ B 2 の時は、 2 --- (27) D≤ B2 A− A となり、この条件は、B = 0 かつ A = 2 ⋅ (1 + β ) の時に最も厳 2 Case IV 0.2 --- (24) 0 ≤ A ≤ 2 ⋅ (1 + β ) d 1 0.6 0.4 --- (23) であれば安定である。 A, B のとりうる値の範囲は、 Rc=2 Rc=4 Rc=6 0 1 ) D ⋅ A2 − B 2 ≤ 2 ⋅ A 0 0.8 --- (22) G ≤ 1 となることが安定条件であるので、 0.2 1 B = sin(θ ⋅ β ) − β ⋅ sin(θ ) --- (21) 2 0.4 0 2 ここで、 A = (1 + β ) − cos(θ ⋅ β ) − β ⋅ cos(θ ) Case II Φ 1.2 --- (20) 増幅係数は、 Fig. 3 Analytical Solution 1.2 ) --- (18) --- (19) 1 ∆t ≤ 2 2 ∆x --- (30) であるので、式(29)は式(30)を格子間に物体境界がある場合 のために修正されたものと見ることができる。ここでは、式 (29)を修正拡散数制限と呼ぶことにする。 一般に高レイノルズ数流体計算の場合にはクーラン数の制 限が厳しくなり、拡散数の制限は問題になることは少ないが、 この修正拡散数制限は、格子点が物体境界に近い場合にクー ラン数の制限よりもはるかに厳しい条件となる。このままで は、時間刻みΔt を実用的でないほど小さくとらなくてはな らないので、式(13)を次のように変更する。 拡散項のスキームの選択 前項の議論により、疑似２次精度の差分スキームは、境界 までの距離が近い時には効果が小さ目に表れることがわか っているため、拡散項には、２次精度の差分スキームを使用 すべきである。ただし、２次精度の差分スキームを陽解法で 使用した時には、以下に述べるように拡散数の制限が非常に 厳しくなることがあるので、対策が必要である。 Fig. 5 のように物体の左側境界が点(+1)と点(0)の間にある 場合、φに関する次の１次元拡散方程式を考える。 3 Copyright © 2000 by JSCFD φ in +1 − φ in 2 −ν⋅ {φ in_ d − (1 + β ) ⋅ φ in+1 + β ⋅ φ in−1 }= 0 d ⋅ (d + ∆x ) ∆t (a) --- (31) ここでは、拡散項の中心項を陰的に扱っている。 この差分式について Von Neumann の安定性解析を行うと、 次の増幅係数を得る。 G2 = (1 + D ⋅ A)2 − (D ⋅ B )2 (1 + D ⋅ C )2 ここで、 A = cos(θ ⋅ β ) + β ⋅ cos(θ ) Pi,j+1 (b) --- (32) (c) --- (33) B = sin(θ ⋅ β ) − β ⋅ sin(θ ) C =1+ β Fig. 6 Pressure at boundary point --- (34) --- (35) 3.5. 2 G ≤ 1 となることが安定条件であるので、 ( ) D ⋅ A 2 − B 2 − C 2 ≤ 2 ⋅ (C − A) − (1 + β ) ≤ A ≤ (1 + β ) --- (36) --- (37) であるので、 C−A≥0 C+ A≥0 --- (38) --- (39) l よって、 A 2 − B 2 − C 2 = −(C + A) ⋅ (C − A) − B 2 ≤ 0 差分スキームのまとめ 以上の議論から、今回使用する差分スキームをまとめる。 スキームの単純化のため、フラグを用いて境界の有無・方向 に関わらず共通に使用できるように構成する。ただし、Fig. 7 に示すように、予め物体内部点での流速をゼロに設定してお くことと、物体境界までの距離 d, e を物体境界の無いとこ ろでは格子間隔の値に設定しておくことを前提とする。 であれば安定である。D は式(20)より正である。 A のとりうる値の範囲は、 対流項およびポアソン式右辺用の１階微分 部分) (u i +1 − u i −1 ) é ∂u ù ê ∂x ú = (d − e ) ë ûc _i --- (40) 式(38)、式(40)より、式(36)は常に成立する。よってこのスキ ームは無条件安定である。修正拡散数の制限を受けない。 l Body φi-1 i d --- (41) 境界での人工粘性項および中心補正項用の２階微分 é (d + e ) ⋅ (u − u )ù − ( flag _ right _ boundary ) ⋅ ê4 ⋅ i +1 i ú 2 ë (d − e) ⋅ d û é ù ( d + e) − ( flag _ left _ boundary ) ⋅ ê4 ⋅ ⋅ (ui −1 − ui )ú 2 ë (d − e ) ⋅ e û i+1 Δx Ø Fig. 5 Left side boundary 3.4. (1 次精度 é ∂ 2u ù (u i +1 − 2 ⋅ u i + u i −1 ) ê 2 ú = 4⋅ (d − e )2 ë ∂x û c _ i φi φi_d i-1 Pi-1,j Pi,j 圧力ポアソン式 圧力は、Fig. 6 の P(i,j)のような物体境界に隣接した物体内 部の格子点で値を持つ。以下、これらの点を境界点と呼ぶこ ととする。ただし、いくら物体境界に隣接していても、Fig. 6 (c)のように斜め方向で隣接する点は差分で使用されないの で境界点に含まれない。 境界点での圧力は複数の方向からの計算に共有されること がある。例えば、Fig. 6 の P(i,j)は、P(i,j+1)中心の差分計算の ために(a)方向で使用され、P(i-1,j)中心の差分計算のために(b) 方向で使用される。このように境界点の圧力を定めることに より、圧力のポアソン式の左辺（圧力に関する差分）は、境 界部でも通常の中心差分で計算することができる。 また、物体内部で定義した仮想的な圧力を複数方向で共有 することは、物体表面上で計算される Cp 分布を滑らかにす る効果がある。 境界点の圧力は、物体境界での境界条件により決められる。 詳しくは次章で説明する。 --- (42) ここで、flag_right_boundary は物体の右側の境界 部、flag_left_boundary は物体の左側の境界部であ れば 1 を、そうでなければ、0 を与えるフラグで ある。 これらを用いて、対流項は次のように構成される。 é ∂ 2u ù d +e æ ∂u ö é ∂u ù u ⋅ç ÷ ≅ u ⋅ ê ú −u ⋅ ê 2 ú ⋅ 2 è ∂x ø i ë ∂x û c _ i ë ∂x û c _ i é ∂ 2u ù æd −eö 3 −α ⋅ u ⋅ ê 2 ú ⋅ç ÷⋅ ∂ x ë ûc_i è 2 ø 4 æ − u i + 2 + 2 ⋅ u i +1 − 2 ⋅ u i −1 + u i −2 + ( flag _ 5 pt _ ok ) ⋅ ç u ⋅ 12 ⋅ ∆x è u − u i +1 − u i −1 + u i − 2 ö --- (43) + α ⋅ u ⋅ i +2 ÷ 4 ⋅ ∆x ø Ø ここで、flag_5pt_ok は物体境界に邪魔されずに 5 点の中心差分を採ることができれば 1 を、そうで なければ、0 を与えるフラグである。 粘性項は次のように構成される。陽解法の場合 ui は、次の 時間ステップの値を使用するために、uin+1 として左辺に移動 する。 l 粘性項用の２階微分 4 Copyright © 2000 by JSCFD é ∂ 2u ù (e ⋅ u i +1 + (d − e) ⋅ u i − d ⋅ u i −1 ) ê 2 ú = 2⋅ d ⋅ e ⋅ (d − e ) ∂ x ë ûi H セルごとの法線ベクトル nω は既知である。あとは、セルご --- (44) との圧力勾配ベクトル (∇p )ω を、セル内の既知の圧力と P0 Body Body を使用して記述すれば、式(48)を解くことができる。 １つのセルに注目した時に、圧力の分布を次式のように仮 定する。ここでは、Fig. 11 のように、注目する境界点を中心 にローカルなｘ､y､ｚ座標を設定している。 --- (49) p ( x, y, z ) = a1 ⋅ x + a 2 ⋅ y + a3 ⋅ z + P0 ui ui-1 = 0 i i-1 ui-1 ui+1 i+1 e ui i i-1 Δx d ui+1= 0 セル内の３点で既知の圧力があれば、式(49)の a1, a2, a3 を 求めることができる。そのような点とは、外部点（物体外の 点）と境界点（物体内であるが圧力を持つ点）である。これ らの圧力は１ステップ前の値を使用することになるが、ポア ソン式の反復解法の中で正しい値に収束する。 セル内にはローカルの原点以外に格子点が７点あるので、 既知の圧力が最大で７点与えられる場合がある。また、最低 ３点は常に利用可能なので、条件過多に備え、最小２乗法を 使用する。既知の圧力を Pi、その圧力を与える点の座標を(xi, yi, zi) とし、次の式を最小にする a1, a2, a3 を求める。 i+1 Δx (e<0, d=Δx ) (d>0, e= -Δx ) Fig. 7 Boundary of right side and left Side 4. 物体境界での境界条件 物体境界での流速の境界条件については、差分スキームに 取り込まれている。ここでは、圧力の境界条件について述べ る。 4.1. 圧力境界条件 Ψ = å ( p (xi , y i , z i ) − Pi ) = å (a1 ⋅ xi + a 2 ⋅ y i + a 3 ⋅ z i + P0 − Pi ) H 物体境界上での圧力の境界条件は、 n を法線ベクトルとし て、次式で与えられる。 H 1 H H ∇p ⋅ n = ∆u ⋅ n Re ∂Ψ = å 2 ⋅ xi (a1 ⋅ xi + a 2 ⋅ y i + a 3 ⋅ z i + P0 − Pi ) = 0 ∂a1 i --- (45) --- (51) a2, a3 についても同様に行い、次のようにまとめる。 H H H Μ ⋅ a = b − P0 ⋅ c æ H 2 çç å Pi ⋅ z i ÷÷ è i ø æ ç å x i ⋅ xi ç i Μ = ç å y i ⋅ xi ç i çç å z i ⋅ xi è i ö åx ⋅y åy ⋅y åz ⋅y i i i i i i çç å z i ÷÷ è i ø ö åi xi ⋅ zi ÷ ÷ åi yi ⋅ zi ÷÷ åi zi ⋅ zi ÷÷ø i i i --- (53) である。 これにより、a1, a2, a3 が求まる。 H H H a = Μ −1 ⋅ b − P0 ⋅ c ( ) --- (54) また、式(49)より、 H ∇p = a --- (55) なので、 H H ∂∇ p ∂a = = − Μ −1 ⋅ c ∂P0 ∂P0 --- (56) 式(54),(55),(56)を、式(48)に代入する。 å 2 ⋅ flag (ω ) ⋅ (− (Μ ω ⑧ −1 =① H H H H H ⋅ c )⋅ n )ω ⋅ Μ −1 ⋅ b − P0 ⋅ c ⋅ n ω = 0 (( ( )) ) --- (57) --- (47) これを、P0 について解く。 =① å flag (ω ) ⋅ ((Μ ω ⑧ 式(47)を最小にする P0 を求めるため、P0 で微分し 0 とおく。 ⑧ æ ∂∇p H ö H ∂Φ = å 2 ⋅ flag (ω ) ⋅ çç ⋅ n ÷÷ ⋅ (∇p ⋅ n )ω = 0 ∂P0 ω =① è ∂P0 øω æ ö ç å Pi ⋅ xi ÷ ç å xi ÷ æ a1 ö ÷, H ç i ÷ H ç ÷, H ç i a = ç a2 ÷ b = ç P ⋅ y ÷ c = ç y ÷ , å å i i i ça ÷ ç i ç i ÷ ÷ è 3ø うに Pi,j を決めたとしても、それが、セル②で式(46)を良く 満たすとは限らない。境界点での圧力を設定する際には、そ れが共有されるすべての境界セルで、式(46)が総合的になる べく良く満たされることを考えなければならない。 ２次元の場合には Fig. 9 のように最大で４セルが、３次元 の場合には Fig. 10 のように最大で８セルが、圧力の境界条 件の下に境界点を共有するセルとなりうる。ここで、セルと は、格子点を中心とするコントロール・ボリュームではなく、 格子点を頂点に持つ四角形（直方体）のことである。それら のセルは、物体境界線（面）と交差した時に、それぞれのセ ルを代表する法線ベクトルを持つ。それらは、圧力の境界条 件の計算に先立ち、予め計算して用意しておくこととする。 以下、３次元の場合を考える。flag(ω)は、①から⑧のセル ω に関して、境界を持つ時に１を、持たない時に 0 を与え るフラグであるとする。セル①を例にとれば、Fig. 11 で点(Δ x,0,0)、点(0,Δy,0)、点(0,0,Δz)のどれかが物体外部点であれ ば、図中のローカルなｘ軸、ｙ軸、ｚ軸上のどれかに物体境 界線が刻まれる。このとき、flag(①)は１となり、それ以外の 場合は 0 となる。 境界点の圧力 P0 は、式(46)を満たさねばならない全てのセ ルでの式(46)の最小２乗法により求められる。すなわち、次 の関数を最小にするように定められる。 ⑧ --- (52) ここで、 を使用して、式(46)が満たされるよ å flag (ω ) ⋅ (∇p ⋅ n )ω ω --- (50) これを a1 で微分して、それを 0 とおく。 これを満たすように境界点の圧力を設定するのが、基本的 な考え方である。 ところが、Fig. 8 を例にとると、セル①において、法線ベク Φ= 2 i ここで、レイノルズ数が極めて小さい場合を除き、実用上 は右辺をゼロと置いても問題ないと言われている。 H --- (46) ∇p ⋅ n = 0 H トル n① と Pi+1,j+1, Pi,j+1 2 i P0 = --- (48) =① ⑧ å ω =① −1 ) ) ⋅ ((Μ H H ⋅c ⋅n ω −1 H H flag (ω ) ⋅ Μ ⋅ c ⋅ n (( −1 H H ⋅b ⋅n ω ) ) --- (58) ) ) 2 ω 5 Copyright © 2000 by JSCFD Fig. 11 Boundary cell No.1 in 3D case これが、境界点における圧力値を与える式である。２次元 の場合も、ほぼ同様に求まる。 4.2. → ｎ① → ｎ② Pi,j+1 Pi+1,j+1 ① ② 薄い物体の場合の２価処理 ２次元で考えると、４つのセルが境界を持つセルであって も、Fig. 9 のように同じ空間に開いている場合には、問題な い。このケースは薄い物体だけでなく、厚い物体の端点でも しばしば表れる。しかし、Fig. 12 のように、上下（または左 右）２つの異なる空間に対して開いている場合には、境界点 の圧力を２つの方向で共有するのは適切ではない。この場合、 １つの境界点について２つの値を持たせるような例外処理 が必要である。 Pi,j ② ① P0 Fig. 8 Pressure boundary condition ③ ② ④ ① Fig. 12 Exceptional Case P0 ③ 4.3. ④ Fig. 9 Boundary cells in 2D case ③ これは、ローカルの原点で -1 の値を与える。そして、グ リッドと物体境界の交点で 0 を与えるように、k1, k2, k3 を 決める。グリッドと物体境界の交点が３点以上利用可能なら ばそれは求められる。境界を含むセルなら、通常は交点が３ 点以上あるはずだが、まれに交点が格子点の真上に来た場合 に、複数の交点が１点に集まったような形になり、３点の利 用ができないことがある。この場合には、微少な値だけ、交 点をずらし、複数の交点に戻す例外処理を行う。そうするこ とで、常に３点以上の交点が利用可能となり、以下を使用で きる。 グリッドと物体境界の交点(i)の座標を(xi, yi, zi) とし、次の 式を最小にする k1, k2, k3 を求める。(i=1∼imax, imax≧3) ② ④ ⑦ ① ⑥ ⑧ ϕ = å (φ (xi , y i , z i ))2 ⑤ i = å (k1 ⋅ xi + k 2 ⋅ y i + k 3 ⋅ z i − 1) Fig. 10 Boundary cells in 3D case 2 --- (60) i ∂ϕ ∂ϕ ∂ϕ = 0, = 0, = 0 より、次式が求まる。 ∂k1 ∂k 2 ∂k 3 z æ ç xi ⋅ xi æ k1 ö ç å i ç ÷ ç k 2 ÷ = ç å y i ⋅ xi çk ÷ ç i è 3 ø ç å z i ⋅ xi ç è i (0,0,Δz) å xi ⋅ y i i åy ⋅y åz ⋅y i i i i i i −1 ö ö æ åi xi ⋅ zi ÷÷ çç åi xi ÷÷ åi yi ⋅ z i ÷÷ ⋅ çç åi y i ÷÷ åi zi ⋅ zi ÷÷ø ççè åi z i ÷÷ø --- (61) H これを用いて、法線ベクトル n は次のように求められる。 P0 (0,Δy,0) x 法線ベクトル 今回の方法では、境界を含むセルごとに法線ベクトルを用 意する必要がある。これを求めるのは一見煩雑そうであるが、 実は、それぞれのセルについてグリッドと物体境界の交点が わかっていれば、次のように容易に算出することができる。 Fig. 11 のようにセルごとにローカルな座標軸を用意し、次 のポテンシャル関数を考える。 --- (59) φ ( x, y , z ) = k1 ⋅ x + k 2 ⋅ y + k 3 ⋅ z − 1 y H ∇φ 1 = n= ∇φ k12 + k 22 + k 32 (Δx,0,0) æ k1 ö ç ÷ ⋅ ç k2 ÷ çk ÷ è 3ø --- (62) 6 Copyright © 2000 by JSCFD 5. 有効性の評価 5.1. 翼型（NACA0012）まわりの流れ（２次元） 迎角 0 度の NACA0012 まわりの流れを、今回の直交座標 による方法と従来の物体境界適合座標による方法とで計算 を行った。Fig. 13 に NACA0012 を配置した直交格子を、 Fig. 14 に物体境界適合格子を示す。両者の最小格子間隔は、 ほぼ同じに取っている。レイノルズ数は 1000 である。 Fig. 15, 16 に両者の等圧線分布を示す。Fig. 17 に直交格子 計算での流速ベクトルを示す。概ね妥当な計算結果が得られ ている。 Fig. 18 に Cp 分布を示すが、両者は前縁で少し異なってい る。この部分は格子点にして数点しかないにもかかわらず、 流れが大きく変化するので、物体形状を計算の中にどう取り 込むかの違いが結果に大きく表れる。 Cd 値を Table 2 に示 す。Fig. 18 と Table 2 には、比較のために、物体から離れ る方向への格子間隔を半分にした物体境界適合格子の計算 結果を合わせて表示した。 Cp 分布や Cd 値の相違は、物体 境界適合格子の解像度が上がるにしたがって、小さくなって いる。これは、通常の物体境界適合格子を使用した計算の場 合、Fig. 19 に示すように、法線方向の勾配ゼロという圧力 の境界条件が、物体境界上ではなく、物体から半格子離れた 点で満たされるため、法線方向の解像度が精度に大きく影響 することを示唆している。Fig. 20 に直交格子計算による前 縁圧力分布を示す。結果的に、同程度の格子間隔で比較した 場合、今回の直交格子計算の方がより良い Cp 分布が得られ ている。 Fig. 15 Pressure contour around NACA0012 (Cartesian grid system) Fig. 16 Pressure contour around NACA0012 (boundary fitted coordinate grid system) Fig. 13 Cartesian grid system around NACA0012 Fig. 17 Velocity vector plot around NACA0012 (Cartesian grid system) 1.2 1 Fig. 14 Boundary fitted coordinate grid system around NACA0012 0.8 BFC Cartesian Cp 0.6 0.4 0.2 code 0 -0.2 0 0.2 0.4 0.6 0.8 1 -0.4 -0.6 7 Copyright © 2000 by JSCFD 1.2 1 0.8 BFC(fine) Cartesian Cp 0.6 0.4 0.2 code 0 -0.2 0 0.2 0.4 0.6 0.8 1 -0.4 -0.6 Fig. 18 Cp distribution around NACA0012 Table 2 Cd around NACA0012 Cd Cdv Boundary fitted 0.1262 0.0760 coordinate grid Boundary fitted 0.1161 0.0795 coordinate grid (fine mesh case) Cartesian grid 0.0968 0.0625 (current) Fig. 21 3D Cartesian grid system around sphere Cdp 0.0502 0.0366 0.0343 Fig. 22 Pressure contour around sphere Fig. 19 Pressure contour near leading edge (boundary fitted coordinate grid system) Fig. 23 Velocity vector plot around sphere Fig. 20 Pressure contour near leading edge (Cartesian grid system) 5.2. 球まわりの流れ（３次元） Fig. 21 に球を配置した直交格子を示す。球の直径に対して、 格子はわずかに 16 点しかとっていないが、Fig. 22,23 に示 すように滑らかな計算結果が得られている。レイノルズ数は 1000 である。Fig. 24,25 に T=80 の時の流線を示す。後流は、 ３次元的な渦構造になっていることがわかる。Fig. 26 に Cd 値について今回直交格子により計算された値と他の計算に よる値および実験による値を比較した。今回は粗い格子を使 用したにもかかわらず、１０%程度の誤差に収まっている。 Fig. 24 Stream line plot around sphere (Side view) 8 Copyright © 2000 by JSCFD Fig. 28 Pressure contour around the complicated object Fig. 25 Stream line plot around sphere (Back view) 1.8 Current Shirayama (Ref.6) Rimon (Ref.7) LeClair (Ref.8) Liebster Horii (Ref.2) 1.6 1.4 Cd 1.2 1 0.8 0.6 0.4 0.2 Re 0 10 100 1000 10000 Fig. 29 Velocity vector plot around the complicated object Fig. 26 Cd plot around sphere 5.3. 複雑形状まわりの流れ（３次元） Fig. 27 に示す複雑形状まわりの流れを計算した。レイノル ズ数は 1000 である。部品である球の直径に対して、格子は わずかに 16 点しかとっていない。Fig. 28 に圧力分布を、Fig. 29 に速度ベクトルを、Fig. 30,31,32 に流線を示す。格子が 粗いにもかかわらず、形状を良く反映した滑らかな流れが得 られている。 Fig. 30 Stream line plot around the complicated object Fig. 27 3D Cartesian grid system around the complicated object 9 Copyright © 2000 by JSCFD 今回の方法は、実用的な精度と安定性を有しながら境 界処理にコストがあまりかからない計算方法となって いる。 7. 参考文献 (1) T.Ye, R.Mittal, H.S. Udaykumar and W.Shvy, “A Cartesian Grid Method For Viscous Incompressible Flows With Complex Immersed Boundaries”, AIAA-99-3312, pp. 547557. (2) 堀井, 西田, 里深, “仮想境界デカルト格子法による物体 まわりの非圧縮性 DNS”, 第 13 回数値流体力学シンポジ ウム (1999) (3) 小野, 姫野, 冨田, 藤谷, “車のエンジンルーム内の実用 的な流れ解析手法の開発”, 第 11 回数値流体力学シンポ ジウム (1997) (4) 橋口, 桑原, “複雑形状をすぎる非定常流の数値計算”, 第 10 回数値流体力学シンポジウム (1996) (5) 市川, 藤井, “直交格子を利用して任意物体形状まわりの 流れを解く場合の境界処理”, 第 13 回数値流体力学シン ポジウム (1999) (6) S.Shirayama, “Flow Past a Sphere: Topological Transitions of the Vorticity Field”, AIAA Journal Vol.30, No.2, pp. 349-358 (1982) (7) Y.Rimon, S.I.Cheng, “Numerical Solution of a Uniform Flow over a Sphere at Intermediate Reynolds Numbers”, Physics of Fluids, Vol.12, No.5, pp. 949-959 (1969) (8) B.P.Le Clair, A.E.Hamielec, “A Numerical Study of the Drag on a Sphere at Low and Intermediate Reynolds Numbers”, Journal of the Atmospheric Sciences, Vol.27, pp. 308-315 (1970) Fig. 31 Stream line plot around the complicated object (Side view) Fig. 32 Stream line plot around the complicated object (Top view) 5.4. ² 物体境界適合座標計算との計算負荷の比較 従来の物体境界適合格子を使用した計算と今回の直交格子 を使用した計算とでＣＰＵの消費時間を計測した。Table 3 に結果を示す。値は 1 ステップあたりの時間を、計算に使用 した格子点の数で除したものである。計測に使用したＣＰＵ は、Intel Pentium II 450MHz である。両者のコンパイラ が異なるため厳密な比較はできないが、圧力のポアソン式の 反復では、直交格子計算の方が約 2 倍速くなっている。これ は、メトリックが無いためと考えられる。一方で、NavierStokes の運動量式の方は、直交格子計算の方が若干遅いと いう結果になった。これは、式(41),(42),(43)で見られるよう に、境界までの距離を差分に取り込むために、スキーム自身 が長くなっているためと推定される。計算全体としては、両 者の計算負荷はおおよそ同じと見積もることができる。 Table 3 CPU time per stencil per step (sec.) Poisson Momentum Compiler equation equation Boundary fitted 1.08×10-6 3.78×10-6 Visual coordinate grid Fortran Cartesian grid 4.61×10-7 4.68×10-6 Visual C++ (current) 6. ² ² ² 結論 境界部での人工粘性項については疑似２次精度スキー ムを、中心補正項については疑似２次精度スキームま たは補正項なし（１次精度）を選択すべきである。 境界部での拡散項には、２次精度スキームを使用し、 その中心項は陰的に扱うべきである。 今回の方法は、計算を単純に行えるという直交格子の メリットを損なわずに、３次元物体の境界を処理でき る。 10 Copyright © 2000 by JSCFD