第 26 回数値流体力学シンポジウム
講演番号

直交格子法を用いた粘性計算における力計算と物体壁面境界の取扱いについて
Force Calculation and Wall Boundary Treatment
for Viscous Flow Simulation using Cartesian Grid Methods○ 高橋 悠一, 東大院, 東京都文京区本郷 7-3-1, 2008136238@mail.ecc.u-tokyo.ac.jp
今村 太郎, 東大工, 東京都文京区本郷 7-3-1, timat@ mail.ecc.u-tokyo.ac.jp
Yuichi Takahashi, University of Tokyo, 7-3-1, Hongo, Bunkyo 113-8656, Japan
Taro Imamura,
University of Tokyo, 7-3-1, Hongo, Bunkyo 113-8656, Japan
A Cartesian grid method is useful for Computational Fluid Dynamics (CFD) because the method can generate grids
easily and rapidly around complex solid boundaries. When a Cartesian grid method is used to compute viscous flow, it
is difficult to calculate velocity gradient on the wall surface for viscous force calculation because grids cannot alignalong the wall surface. When using wall boundary condition and force calculation method which consider the locationof real wall boundary by velocity interpolation, viscous force on the wall surface can be calculated accurately.
It is essential to improve both wall boundary condition and force calculation method in order to calculate viscous force
accurately using Cartesian grid methods.
１． 序論
航空機の設計現場における数値流体力学（CFD）の利用において，短時間で計算が可能であることが望まれる．特に複雑な物体形状の解析においては流体計算のための格子生成にかかる時間と人手が問題である．
近年注目されている手法に直交格子法がある．
直交格子法は複雑物体形状に対しても完全自動で格子生成が可能であり，短時間での生成が可能である．
しかし粘性計算に直交格子法を用いる場合，境界層の解像が困難になる．直交格子は一般的に格子が物体壁面境界と一致しないため，構造格子や非構造格子と比較して物体壁面近傍での計算精度が問題となる．部分境界適合直交格子法を用いると境界層の解像は可能であるが，物体形状が複雑になると格子生成が困難になる場合もあり直交格子法が本来持つ利点が損なわれる．
そこで本研究は，直交格子法のみによる粘性計算を行った場合に問題となる，境界層速度分布や壁面上の速度勾配の推算精度，また力計算法に対する壁境界条件の影響について体系的に調べた結果を報告する．平板と円柱を用いた検証計算を行い，理論解や過去の計算結果との比較を行う．

課す．
Fluid CellがWall Cellに隣接するとき，Fluid Cell，Wall Cellをそれぞれセル i，セル j とする．
(1) 階段状の壁境界条件（壁境界条件(1)）
Fig. 1 のように，セル i とセル j の境界面を物体壁面とみなして壁境界条件を与える．セル j について，密度・圧力はセル i での値を代入し，速度はセル i での速度を符号を反転させて代入する．

時間積分
時間刻み法

(1)

u j  ui , v j  vi(2)(2) 速度の一次内挿による壁境界条件（壁境界条件(2)）
Fig. 2 のように物体壁面の法線上に Image Point (IP) を定義
し，物体壁面の位置を考慮した境界条件を与える．ここで，d0，dIP は物体壁面からセル j の中心，IP までの距離を表す．セルj について，密度・圧力は IP での値を代入し，速度は IP と壁面での値から一次内挿を行う(3)(4)．ただし，dIP は最小格子幅をΔxmin，物体壁面の単位法線ベクトルを n(nx，ny)として，式(6)のように決定する(4)．

２． 数値計算法
現在研究室で開発中の圧縮性流体ソルバー(1)(2)を用いて計算を行う．Tab. 1 に数値計算法をまとめる．
Tab. 1 Numerical calculation method
格子生成
四分木法を用いた直交格子法
支配方程式
空間の離散化
物理量勾配の評価
非粘性流束の評価
粘性流束の評価 j   i , p j  pi

二次元圧縮性ナビエ・ストークス方程式
セル中心有限体積法
Green - Gauss Method
SLAU スキーム
左右のセルの内分でセル境界面での物理量と物理量勾配を評価
LU-SGS 法
局所時間刻み法 j   IP , p j  p IP(3)uj d0u IPd IP(4)vj d0v IPd IP(5)d IP xmin

max n x , n y(6)(3) 速度の二次内挿による壁境界条件（壁境界条件(3)）
Fig. 3 のように物体壁面の法線上に IP1 と IP2 を定義し，(2)と同様に物体壁面の位置を考慮した境界条件を与える．ここで，d0，d1，d2 は物体壁面からセル j の中心，IP1，IP2 までの距離を表す．セル j について，密度・圧力は IP1 での値を代入し，速度は IP1，IP2，壁面での値から二次内挿を行う(5)．
物体壁面に沿う座標系を XY とし，X，Y 方向の速度をそれぞれ U，V とすると，

３． 壁境界条件の設定方法と力の計算方法
３.１ 壁境界条件の設定方法
物体壁面と交差するセルを Wall Cell，流体側のセルを Fluid
Cell とする．壁境界条件は Wall Cell に持たせる物理量によって決まり，粘性計算では物体壁面に滑りなし条件と断熱壁条件を
1

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 j  1 , p j  p1

(7)

U j  d 0 Au  d 0 (d o  d1 ) Bu(8)V j  d 0 Av  d 0 (d o  d1 ) Bv

(9)

３.２ 力の計算方法
粘性計算における力は，物体壁面上の圧力と粘性力を物体周りで積分することで計算できる．Fluid Cell が Wall Cell に隣接するとき，Fluid Cell，Wall Cell をそれぞれセル i，セル j とする．
(1)セル境界面に沿った階段状の力計算（力計算法(1)）
Fig. 4 のように，階段状壁に沿った計算を行う．セル i とセル j の境界面 k にセル i が持つ密度，圧力，速度勾配を持たせ，境界面 k に対して垂直方向に圧力 pk，平行方向に粘性力 τk を与える．境界面 k の面積を sk として，物体周りの力 f は以下のように計算できる．

U 2  U1
 Au UddA  1 , B  2
1u u d1d2V2  V1 Av 
Vd  d1
 Av  1 , Bv  2

dd21



f   p k  τ k   s k(12)kとなる．
Uj，
Vj は x 方向，y 方向の速度になるように回転する．
ただし， d1，d2 は最小格子幅をΔxmin として以下のように決定する．

d1  1.5xm i n(10)d 2  3.0xm i n(11)(2)物体形状データに射影した力計算（力計算法(2)）
Fig. 5 のように実際の物体壁面に沿った計算を行う．物体形状データの各要素 k に周囲のセルが持つ物理量を射影し，要素k に対して垂直方向に圧力 pk，平行方向に粘性力 τｋを与える．
要素 k の面積を sk として，物体周りの力 f は以下のように計算できる．

f   p k  τ k   s k

(13)

k

物体形状データの要素への射影方法を示す．ここで，物体壁面に沿った座標系を XY とし，X 方向の速度を U とする．Fig. 6のように物体壁の各要素の中心を通る法線上に IP1，IP2 を定義し，これら２点から壁面上での密度，圧力，速度勾配を求める．
密度，圧力は IP1 での値を代入する．速度勾配は式(15)のように２点の IP と壁面を通る二次関数を求め，式(16)のように勾配を計算する(5)．

 0  1 ,Fig. 1 Illustration of wall boundary condition (1)p0  p1(14)U (Y )  YA  Y (Y  d1 ) B A  U1 ,
d1

(15)

U 2  U1 Ad 2  d1

B

d2



 U 
  U (0)  A  d1 B
 Y  0(16)Fig. 2 Illustration of wall boundary condition (2)Fig. 4 Illustration of force calculation method (1)Fig. 3 Illustration of wall boundary condition (3)2Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号Fig. 5 Illustration of force calculation method (2)(a)

Fig. 6 Illustration of Image Point for method (2)４． 検証計算１：平板周り流れの検証計算
４.１ 目的３章で述べた壁境界条件(1)から(3)において平板周り流れの検証計算を行い，各壁境界条件における層流境界層の速度分布と壁面上の摩擦抗力係数を求め，理論解と比較する．平板と格子のなす角θ(Fig. 7)を変化させて，直線状に表現した場合と階段状に表現した場合で計算する．
摩擦抗力係数を求める際，壁面上の速度勾配の計算は３．２
節の力計算法(2)と同様に壁面，IP1，IP2 の速度から二次関数を求め，式(16)により求める．

(b)
Fig. 8 Grids generated around the plate

４.３ 計算結果と考察
Fig. 9 から Fig. 11 は，壁境界条件(1)から(3)のもとで各レイノルズ数における x  1.0 での層流境界層の速度分布を各回転角について計算し，理論解と比較したものである．U，Y はそれぞれ物体壁面に対して平行方向の速度と垂直方向の壁面からの距離を表す．Fig. 12 は平板周りのマッハ数分布の例を示す．
どの壁境界条件においても，すべてのレイノルズ数で回転角によらず定性的に理論解に一致する．Re = 104 までは従来と同程度の格子幅を用いれば，壁境界条件(1)のような階段状の壁境界でも層流境界層の速度分布は計算可能である．
同一の条件で各壁境界条件における x  1.0 での摩擦抗力係数 Cf の値を計算し，理論解と比較したものを Tab. 3 から Tab. 5に示す．各壁境界条件で二乗平均誤差を求め，比較を行う．壁
境界条件(1)における理論解との二乗平均誤差は７５％である．
一方壁境界条件(2)における理論解との二乗平均誤差は８％，壁境界条件(3)では１８％である．Fig. 9 から Fig. 11 で示したように，壁境界条件(1)，(2)，(3)で境界層全体のプロファイルには明確な違いは観察されなかったが，Cf については壁境界条件(1)では計算誤差が大きく，壁境界条件(2)，(3)では改善される．これは壁面での速度勾配の違いが要因と考えられる．Fig. 13 は，Re=103 ，θ=45°において式(14)で与えられる速度の二次関数分布を，各壁境界条件について示したものである．壁境界条件(1)と比較して壁境界条件(2)と(3)で壁面における速度勾配がより理論解に近いことが確認できる．

θFig. 7 Illustration of parameter“θ”
４.２ 計算条件
計算条件を Tab. 2 にまとめる．最小格子幅は，構造格子で境界層を解像する際に用いられる格子幅と同程度にする．θ=0°，45°の場合に生成される格子を Fig. 8 に示す．
Tab. 2 Computational condition for case 1
レイノルズ数 Re
10 2 103 104
マッハ数 M
最小格子幅Δxmin
平板の回転角θ（deg）
CFL 数

0.20.1/ Re

0 15 30 45
1003Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号Tab. 3 Skin friction coefficient on the plateusing wall boundary condition (1)ReCfCfθ理論解との計算値
誤差(％)
理論解00.032-52
10 20.0664150.037-44
300.029-56
45-0.016-124
300.011-48
100.0210150.0033-84300.0046-78
450.0020-91
400.0020-70
100.00664
150.0030-53
300.0020-70450.00048-9375
二乗平均(a)Tab. 4 Skin friction coefficient on the plateusing wall boundary condition (2)Re
CfCfθ
理論解との計算値
誤差(％)
理論解00.070510 20.0664
150.0694300.070
5450.080
2100.021
01030.0210
150.0225300.022
5450.02310400.0066-1
100.00664
150.0087
31300.0075
13450.007288
二乗平均(b)Tab. 5 Skin friction coefficient on the plateusing wall boundary condition (3)Re
CfCfθ
理論解との計算値
誤差(％)
理論解00.070510 20.0664
150.080
21300.070
5450.0923900.021
0103
0.0210150.0210300.022
5450.026
2400.0076
15104
0.00664150.0085
28300.0079
19450.0078
1818
二乗平均(c)Fig. 9 Boundary layer profile
using wall boundary condition (1)((a): Re=102,(b):Re=103,(c):Re=104)

4

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号

(a)

(a)

(b)

(b)

(c)
Fig. 10 Boundary layer profile
using wall boundary condition (2)((a): Re=102,(b):Re=103,(c):Re=104)(c)Fig. 11 Boundary layer profile
using wall boundary condition (3)((a): Re=102,(b):Re=103,(c):Re=104)

5Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号Fig. 12 Mach number distribution around the plate

Fig. 14 Grids generated around the cylinder
５.３ 計算結果と考察 1
円柱周りの流れ場の再現性を検証する．各レイノルズ数における円柱周りのマッハ数分布と流線を Fig. 15 に示す．後方で双子渦ができており，レイノルズ数の増加に伴い大きさを増すことが確認できる．各レイノルズ数におけるよどみ点圧力係数
Cp,s，再循環領域の長さ L と円柱の半径の比 L/0.5D の計算結果
を過去の研究(6)(7)(8)と比較を行うと，Tab. 7 のようになる．いずれのレイノルズ数においても壁境界条件(1)と(2)共に過去の研究結果の範囲内に収まっており，流れ場は正しく計算できている．
抗力係数 Cd を力計算法(1)と(2)により求めた結果を過去の研究(6)(7)(8)と比較すると Tab. 8 のようになる．力計算法(1)では，レイノルズ数によらず壁境界条件(1)，(2)ともに１０％以上値が小さい．力計算法(2)では，壁境界条件(1)の場合は力計算法の改善にもかかわらず過去の研究との差が３０％に増大する．一方壁境界条件(2)の場合は力計算法を(1)から(2)に改善することで過去の研究との差は１％程度となる．つまり，流体計算に壁境界条件(2)を，力計算に力計算法(2)を用いた場合に限り，抗力係数が過去の計算例と一致することが分かる．
そこで，流体計算と力計算の方法の組み合わせと抗力係数の計算結果を詳しく調べる．Tab. 9 は Re = 10 における力計算の結果を，圧力抗力係数 Cdp と摩擦抗力係数 Cdv に分解して示したものである．これから各計算結果の差は摩擦抗力係数のみにみられることが分かる．平板の例で示した通り，壁境界条件(1)での流体計算では Cf が平均で１/４程度の値に見積もられるため，力計算法(2)を用いても Cdv が小さく見積もられる．一方，壁境界条件(2)では Cf が８％程度の誤差で見積もられるため，力計算法(2)を用いることで摩擦抗力の計算精度が向上されると考えられる．
よって直交格子法を用いた粘性計算は，壁境界条件と力計算共に改善を行うことで抗力係数の推算精度の向上が可能となる．

Fig. 13 Velocity distribution near the wall boundary
using each wall boundary condition
５． 検証計算２：円柱周り流れの検証計算
５.１ 目的
直径１の円柱周り流れの検証計算を行い，流れ場の再現性と力計算の精度について検証する．流れ場は円柱後方に発生する双子渦の大きさやよどみ点での圧力係数について比較する．力計算は４章で述べた力計算法(1)と(2)で計算を行い，計算精度の違いを調べる．このとき，壁境界条件を(1)と(2)に設定した場合について計算を行う．
５.２ 計算条件
計算条件を Tab. 6 にまとめる．レイノルズ数は定常流れが保たれる範囲内で設定する．円柱周りに生成される格子を Fig. 14に示す．
Tab. 6 Computational condition for case 210 20 40
レイノルズ数 Re0.2
マッハ数 M0.01
最小格子幅Δxmin
100
CFL 数

6Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号Tab. 7 Pressure coefficient at stagnation point and length ofrecirculation zone around the cylinder
ReCp,sL/0.5D
著者
10Dennis and Change
1.489
0.53
Nieuwstadt and Keller
1.500
0.434
1.403
0.478
Imamura ら
1.426
0.501
本研究(境界条件(1))
1.445
0.477
本研究(境界条件(2))
20
Dennis and Change
1.269
1.88
Nieuwstadt and Keller
1.274
1.786
1.233
1.852
Imamura ら
1.225
1.860
本研究(境界条件(1))
1.242
1.805
本研究(境界条件(2))
40
Dennis and Change
1.144
4.69
Nieuwstadt and Keller
1.117
4.357
1.158
4.454
Imamura ら
1.114
4.458
本研究(境界条件(1))
1.128
4.360
本研究(境界条件(2))

Re102040(a)Tab. 8 Drag coefficient of the cylinder
Cd
著者
Dennis and Change
2.846
Nieuwstadt and Keller
2.828
2.848
Imamura ら
力計算法(1)
力計算法(2)2.503
2.158
本研究（壁境界条件(1)）
2.491
2.803
本研究 (壁境界条件(2)）
Dennis and Change
2.045
Nieuwstadt and Keller
2.053
2.051
Imamura ら
力計算法(1)
力計算法(1)1.826
1.591
本研究（壁境界条件(1)）
1.823
2.039
本研究（壁境界条件(2)）
Dennis and Change
1.522
Nieuwstadt and Keller
1.550
1.538
Imamura ら
力計算法(1)
力計算法(2)
1.379
1.221
本研究（壁境界条件(1)）
1.380
1.538
本研究（壁境界条件(2)）(b)

Tab. 9 Drag coefficient of the cylinder at Re=10 divided into
pressure drag coefficient and viscous drag coefficient
力計算
力計算法(1)
力計算法(2)
壁境界条件
Cd = 2.506
Cd = 2.178
壁境界条件(1)Cdp = 1.587
Cdp = 1.544
Cdv = 0.919
Cdv = 0.633
Cd = 2.489
Cd = 2.803
壁境界条件(2)Cdp = 1.582
Cdp = 1.543
Cdv = 0.908
Cdv = 1.260(c)Fig. 15 Mach number distribution and flow line
around the cylinder((a):Re = 10, (b):Re = 20, (c):Re = 40)

7

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号
６． 結論
直交格子法による粘性計算を行った場合に問題となる境界層速度分布や壁面上の速度勾配の推算精度について，３つの異なる壁境界条件で比較を行った．また，２つの力計算法に対する壁境界条件の影響について体系的に調べた．
層流境界層速度分布は，レイノルズ数 104 の計算までは階段状の境界条件でも再現可能で，理論解に一致する．一方で壁面上の速度勾配は壁面のごく近傍の速度分布をとらえることが必要で，速度の一次内挿や二次内挿の計算による壁面の位置を考慮した境界条件に改善することで理論解と同程度の値となる．
抗力係数のうち，圧力抗力係数の計算結果は壁境界条件と力計算法によらず一定値となるが，摩擦抗力係数の計算結果には差が見られる．セル境界面に沿った階段状の力計算では，壁境界条件によらず摩擦抗力係数が正しく計算できない．物体壁面に沿った力計算法に改善すれば，壁面の位置を考慮した境界条件を用いて壁面上の速度勾配を正しくとらえられる場合に限り，摩擦抗力係数が計算可能となり抗力係数の値が過去の計算結果と一致する．
直交格子法による粘性計算は，壁境界の計算方法と力計算法を共に改善することで抗力係数の推算精度の向上が可能となる．

謝辞
本研究は科学研究費助成事業，若手研究（B）（基金）
・23760767に基づき実施された．ここに感謝の意を表す．
参考文献(1) 今村太郎，神園仁志，高橋悠一，“航空機空力性能評価を目的とした直交格子法ベース流体ソルバーの開発，”平成 23 年度航空宇宙空力班シンポジウム，2012 年 1 月.27 日．
(2) 神園仁志，今村太郎，“埋め込み境界法を用いた直交格子オイラーソルバーの改良，”第 44 回流体力学講演会・航空宇宙数値シミュレーション技術シンポジウム 2012．
(3) Mittal, R. and Iaccarino, G., “Immersed Boundary Methods,”
Annual Review of Fluid Mechanics, Vol. 37, 2005, pp.239-261.
(4) 中橋和博, “直交格子による圧縮性非粘性流れ計算のための壁境界の取り扱い,”第 25 回数値流体力学シンポジウム講演論文集,2011．
(5) Aftosmis, M. and Berger, M., “Progress Towards a Cartesian Cut –Cell Method for Viscous Compressible Flow,” AIAA Journal, January
2012.
(6) Dennis, S.C.R. and Change, G.Z., “Numerical solutions for steadyflow past a circular cylinder at Reynolds number up to 100,” Journal ofFluid Mechanics., Vol. 42, 1980, pp471.
(7) Nieuwstadt, F. and Keller, H.B., “Viscous flow past circularcylinders,” Computers and Fluids, Vol. 1, 1973, pp59-71.
(8) Imamura, T., Suzuki, K., Nakamura, T. and Yoshida, M.,“Acceleration of steady-state lattice Bolzmann simulations on
non-uniform mesh using local time step method,” Journal ofComputational Physics, Vol. 202, 2005, pp.645-663.

8Copyright © 2012 by JSFM