第 27 回数値流体力学シンポジウム
講演番号 E10-2

埋め込み境界法と SA 乱流モデルに基づく壁関数を用いた定常流解析
Steady Flow Simulation using Immersed Boundary Method
and Wall Function based on SA Turbulence model○ 高橋悠一, 東大院, 東京都文京区本郷 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
The purpose of this study is to develop a method which enables to simulate fully steady turbulent flow around complex
geometries using Cartesian grid. An immersed boundary method is coupled with a wall function, so that highly
clustered grid in wall normal direction which is necessary for turbulent flow simulation can be alleviated. This isessential when Cartesian grid is used because the grid can only be refined equally in all directions. The results fromflow simulations around the flat plate show that turbulent boundary layer profile can be resolved with y+≈120 at thefirst cell above the wall when the plate is oriented tangential to the grid. Also, the profile can be resolved with y+≈30
when the plate is oriented 45° to the grid.

１． 序論
航空機設計における数値流体力学（CFD）の利用頻度が高まる
中で，計算時間の短縮化が大きな課題となっている．そのため格子生成にかかる時間は CFD の主要な問題の一つであり，直交格子法が近年再注目されている．直交格子法は，複雑な物体形状に対しても完全自動で格子生成が可能であるという利点を持ち，生成時間を大幅に短縮できる．
しかし直交格子法を乱流計算に利用する際には，格子数の大幅な増大が問題である．乱流境界層の解像には物体壁面近傍で格子を十分に細分化する必要があるが，壁面垂直方向のみに細かい格子を寄せられる物体適合格子に対し，直交格子は格子が等方的に分割される．その結果，物体適合格子の壁面垂直方向の格子幅と同等の格子幅を有する格子を直交格子で生成すると，例えば二次
元平板乱流境界層の場合は物体適合格子と比較して10倍から100倍程度の格子数が必要である．この課題を克服するために，物体近傍にのみ物体適合格子を生成し直交格子と組み合わせる方法がある(1)．この方法は物体近傍に細かな格子を寄せられるが，物体形状が複雑になると格子生成が困難になり，格子生成の容易さという直交格子の最大の利点を損ねる．
直交格子のみを用いる方法で格子増大の問題を避けるために，
境界層解像に必要な格子数を削減する手法が考案されている．壁関数（wall model, wall function）を用いる手法(2)(3)，偏微分方程式に基づく２層モデルによる手法(4)である．前者の手法では，理論的な境界層プロファイルを壁近傍で定義し，背面の直交格子とカップリングをする．後者の手法では，境界層方程式などの一次元の偏微分方程式を壁近傍で解く．これらの手法を文献(2)ではカット
セル，文献(3)，(4)では埋め込み境界法と組み合わせて用いており，平板を用いた検証計算で境界層解像に必要な格子数を減らすことに成功している．しかし，摩擦抗力係数の推算精度や格子に対して傾けた平板における境界層の解像に関して課題があり，直交格子法に適した新しい乱流計算手法の開発が望まれている．
本研究では、複雑な物体形状周りにおける定常乱流の計算を可能とするために、直交格子に適した定常乱流計算手法を開発することを目的とする．乱流モデルには Spalart-Allmaras (SA) 乱流モ
デルを用い，SA 乱流モデルから導かれた壁関数(2)を壁近傍で定義する．壁関数とカットセルを組み合わせる文献(2)とは異なり，埋め込み境界法と組み合わせることで境界層解像に必要な格子数の増大を抑える．平板乱流境界層の速度分布や摩擦抗力係数の推算値を理論解と比較し，新しい乱流計算手法を検証する．また，壁関数の導入による必要格子数軽減の効果，物体壁面の格子に対する位置・角度の計算結果への影響を調べる．
２． 数値計算法
2.1 圧縮性流体ソルバー
現在研究室で開発中の二次元圧縮性流体ソルバー(5)(6)に乱流モデルを組み込む．Tab. 1 に数値計算法をまとめる．
格子生成は四分木法を用いた直交格子法によって行う．
ここでは Fig. 1 のように，格子を wall cell・body cell・fluid cell の 3 つの種類に分類する．wall cell は物体壁面に交差する格子，body cell・
fluid cell はそれぞれ物体側・流体側の格子である．流体計算は
fluid cell に対して行う．
Tab. 1 Numerical calculation method
Grid generation
Cartesian grid method
Governing equation
2-D compressible RANS equation
Turbulence model
Spalart-Allmaras turbulence modelSpecial discretization
Cell-centered finite volume method
Physical quantity gradient valuation Green - Gauss Method
Inviscid flux valuation
SLAU scheme
Viscous flux valuation
α-Scheme
Time integration
LU-SGS method

Fig. 1 Cell definition
支配方程式は，式(1)で表される二次元圧縮性 RANS 方程式を用いる．

Q F  FV  G  G V 


0
t
x
y
1

(1)

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2
Q は保存量ベクトル Q =(ρ, ρu, ρv, e)T で，ρ は密度，(u, v) はそれ
ぞれ x 方向，y 方向の速度成分，e は全エネルギーを表しており，それぞれレイノルズ平均が施されている．F，G はそれぞれ x 方
向，y 方向の非粘性流束，Fv，Gv はそれぞれ x 方向，y 方向の粘性流束で，式(2)，式(3)のように表される．

 v  uv G 2
v  p  (e  p )v  u 
 2
 u  p ,F
uv  (e  p)u 

Cb1  0.1355,   2 / 3, Cb 2  0.622,   0.41
Cw1  Cb1  2  (1  Cb 2 )  , Cw2  0.3, Cw3  2,
Cv1  7.1

離散化はセル中心有限体積法を用い，非粘性流束，粘性流束の評価はそれぞれ SLAU scheme(8)，α-scheme(9)(10)によって行う．時
間積分には LU-SGS 法(11)(12)を用いる．
物体壁面は埋め込み境界法によって取り扱う．壁関数とのカ
ップリングを行うことで乱流境界層解像に必要な格子数の軽減を行う．詳細は次節に示す．

(2)

0
0








 xx
 yx




Fv  
 , Gv  
 (3)


xy
yy




 u xx  v xy  q x  u xy  v yy  q y 

2.2 壁関数のカップリング方法
有限体積法では，各セル境界面を通過する流束をその境界面を共有するセルの物理量から評価を行う．wall cell とそれに隣接する fluid cell の境界面での流束を求める場合は，wall cell または境界面における物理量を境界条件として与える必要がある．
本研究では，埋め込み境界法と壁関数を用いて境界面の中心における物
理量を与える．以下にその手順を示す．
まず Fig. 2 のように，wall cell と fluid cell の境界面中心（点 0）を通るような物体壁面の法線上に，Forcing point（点 F）を定義する．物体壁面から点 0，点 F までの距離をそれぞれ y0 ，yF とし，yF = 1.5 Δxmin (Δxminは最小格子幅)とする．
yFは点Fがfluid cell内に位置するように，最小格子幅の対角線長さより長く定義した(2)．点 0 における各物理量は，点 F の物理量を用いて以下に示す方法で与えられる．なお点 F の物理量 qF (ρF, uF, vF, pF) は，点 Fが存在する格子中心における物理量 qFc とその勾配∇qFc から式(13)ように計算する．ここで，d はセル中心から点 F に向かうベ
クトルを表す．

ここで，粘性応力とレイノルズ応力の和である全応力と，乱流熱流束は以下のように表される．ただし，μ は分子粘性係数，μt は渦粘性係数，Pr はプラントル数，Prt は乱流プラントル数，M∞は一様流マッハ数，Re はレイノルズ数である．

2 M
(   t )(2u x  v y )3 ReM xy   yx (   t )(u y  vx )
Re
2 M
 yy

(   t )(2v y  u x )
3 Re

(4)

 c p t c p Tq  
Prt 
 Pr

(5)

 xx



q F  q Fc  q Fc  d乱流モデルには，式(6)で表される SA 乱流モデル(Spalart-Allmaras
turbulence model)(7)を用いる．変数~ は，渦粘性係数 μt と式(7)のような関係が成り立つ．

(13)

2

 ~   ~

1 M   
~ ~~   
u j (   ) cb 2  x j  tx j Re  x j x j 

 


 

+Production + Destruction

t  ~f v1 ~  0 , f v1 3

  C v13
3

(6)

, 

~


(7)

Fig. 2 Definition of point F and point 0
物体壁面に平行な速度成分 u//は，壁関数より求める．計算方法の概略を Fig. 3 に示す．壁関数には SA 乱流モデルから導かれた SA wall model(2)を用いる．SA wall model は以下の式(14)で表され，壁法則に一致する関数である．

また，Production と Destruction は式(8)から(12)で表される．d は
壁面からの距離を表す．

~
Production Cb1S ~(8)M
Destruction  Cw1  f w 2
Re
d





M  ~rg  r  Cw2 r  r ,~Re S  2d 2
6
 C Arctany  a , b  C Arctany  a , b 





3

+

1

1

4



2

2

+

ここで，y = y u τ / ν ，u = u// / uτ である．y は物体壁面からの距
離，ν は動粘性係数，u//は物体壁面に平行な速度成分，uτ は摩擦速度を表す．また，各定数は以下の通りである．

(10)

1

 1 C 6  6w3

f w  g g6  C 6 
w3 




2
2




 B  C1 log y   a1  b12   C2 log y   a2  b2 2  (14)





(9)

v u ~M ~f
S   2v 22 , S ,S  S Sx yRe  df v2  1 ,
1  f v1

 

u y

~ 2

B  5.0333908790
505579(11)(12)

a1  8.1482215800
24245

a2  6.9287093849
022945b1  7.4600876082
527945b2  7.4681457904
01841

c1  2.5496773539
754747c2  1.3301651588
535228c3  3.5994591093
32379

c4  3.6397531868
684494

任意の点において，
適切な摩擦速度 uτ を与えれば y+と u+は式(14)各パラメータは以下の通りである．
2

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2
みに格子が存在するため，
この式を用いて壁面近傍の渦粘性を表す．~ と νt は式(7)を変形した式(23)のような関係が成り立つ．
そこで，式(23)を χ の νt / ν についての近似関数（式(24)）に書き改めることで，混合長モデルによって得られた νt から直接~ を計算する．

を満たす．またそのような uτ は一意に定まる．このことを利用して，点 F での壁面に平行な速度成分 u//F と壁面・点 F 間の距離
yF が式(14)を満たすように uτ を決めることができる．uτ に任意の初期値を代入して点 F での u+，y+（uF+，yF+）（式(15)）を計算し，Newton Iteration により uF+，yF+が(14)式を満たすような uτ を求める．

y F  y F u  F ， u F  u // F u~ t  (y) 2 S D,

(15)

求めた摩擦速度 uτ と壁関数を用いて，点 0 での速度 u//0 は以下のように評価できる．

 





u // 0  u 0  u  u  y 0  u

t
4

,

 3  C v13

(16)


0

また，
求めた摩擦速度 uτ は u//の法線方向の勾配(∂u// /∂y)0 の評価にも使うことができる．まず点 0 における壁関数の勾配を式(17)のように計算し，uτ を用いて速度勾配(∂u// /∂y)0 を式(18)のように
求める．









 u  
2 y0  a1
2 y0  a2
    C1

C
2
2
2
2
2
y0  a1  b1
y0  a2  b2
 y 0
b1
b2
 C3
 C4
2
2
2
2
y0  a1  b1
y0  a2  b2







 u u    y  0

2









(17)



 u    y  0(18)

壁面に垂直な速度成分 u⊥0 は，点 F での値 u⊥F の線形内挿で与
える．
(19)

圧力 p は点 F と等しい値を与える．密度は Crocco-Busemannの式から温度を求めた後に，状態方程式から求める．

p0  p F
T0  TF 

2  1u F2  u //2 0  ，

 0  p0 T0

~


(23)



 0  t  0.2 







 0.2  t  5.0 




(24)




 5.0  t  20



t 

 20  
３． 検証計算
Figure 4 に示される長さ 2 の平板周りの流れ場を解析する．乱流境界層の速度分布・摩擦抗力係数分布を理論解と比較し，乱流境界層の解像に必要な格子数を調べる．
さらに物体壁面の格子に対する位置・角度を変化させ，境界層プロファイル・摩擦抗力係数の計算結果への影響を調べる．
計算条件を Tab. 2 に示す．レイノルズ数を 1.0×106，マッハ数を 0.2，迎角を 0 (deg)とする．レイノルズ数の代表長さは，コード長の半分の 1 とする．
物体壁面の格子に対する位置は，Fig. 5 に示されるように fluid
cell・wall cell の境界面と物体壁面の y 方向の距離 d で表す．角
度は，Fig. 6 の θ で表す．d と θ がともに 0 のとき，fluid cell・wallcell の境界面と物体壁面は完全に一致し，物体壁面に沿った格子となる．Fig. 7 に計算格子の一例を示す．

(20)





以上のように，wall cell の境界面では壁関数と混合長モデルを用いて速度と渦粘性をそれぞれ評価し，fluid cell における RANS
計算とカップリングをする．このとき摩擦速度 uτ を正しく計算できるように点 F は対数領域内に位置させる必要があり，これが格子幅を設定する際の制約となる(2)．また，壁関数は物体外部(y＞0)でのみ使用可能であるが，本研究では境界条件の定義点をfluid cell と wall cell の境界面中心にすることでこの条件を満たしている．

u0   yo / y F  uF

0.2509 
4.39 t 

 

4
3
2

 t 
 t   t 
 0.0168   0.2 13        


 
 
 2.91 t   2.52
 


2






0.0075 t   0.651 t   4.79
 
 

 
 

 t

 

（ y  y0u  0 ， u  u // 0 u ）

0



3

D  1  e x 
p y  A   (22)



(21)

さらに wall cell に隣接する fluid cell について，u//の法線方向の勾配を Green-Gauss Method の計算結果から壁関数で与えられる
勾配に置き換える．セル中心と壁面の距離，セルでの壁に平行な速度成分から，(17)式，(18)式を用いて速度勾配を求めればよい．

Tab.2 Computational condition for flow simulation
Reynolds number
1.0×106
Mach number
0.2
Angle of attack (deg)0Fig. 3 Illustration of velocity computation at point 0 using wall model
渦動粘性係数 νt についての混合長モデル(13)(14)（式(22)）を用いて~ を評価する．ここで，κ = 0.4，A+ = 25， S~ は歪み速度の大きさを表す．この式は LES 解析で壁面近傍を表現するのに多く用いられる．
本研究も粗い格子の計算では対数領域以上の領域のFig. 4 Illustration of flow simulation around the flat plate3

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2

Fig. 5 Definition of “d”(a) Velocity profile at x=1.0
Fig. 6 Definition of “θ”(a) θ = 0 (deg)(b) Skin friction coefficientFig. 8 Flow simulation around the flat plate
without wall model, compared with theoretical value
3.2 壁関数を導入した計算
壁関数を導入して計算を行う．
最小格子幅を壁関数を用いない場合に対して 2 倍，4 倍，16 倍，64 倍に変化させ，格子数をどこまで抑えられるかを調べる．物体壁面の格子に対する位置 d
は 0，0.5Δxmin，Δxmin (Δxmin は最小格子幅 )の 3 通り，角度 θ は0，45(deg)の 2 通りに変化させる．
(1) 格子を壁面に一致させた場合の計算
各最小格子幅での格子数，細分化レベル，壁から一点目における y+を Tab. 3 にまとめる．
Tab. 3 Characteristics of the computational grid around the flat plate
Minimum
Number of Refinement
y+ at the first cell
grid size
cells
level
above the wall
1.95×10-3
8320
10
45.4
4.88×10-4
42880
12
10.6
1.22×10-4
13504014
2.99-5
6.10×1027020815
1.72(b) θ = 45(deg)Fig. 7 Computational grid3.1 すべりなし条件のみによる計算
壁関数を使わずに通常のすべりなし条件で計算を行い，RANS
計算の妥当性を検証する．
物体適合格子で乱流計算を行う場合には，一般に物体壁面から一点目のセルにおける y+は 1 に設定する．そこで本計算でも一点目における y+が物体適合格子と同程度の 1.05 になるように，最小格子幅を 3.05×10-5（格子細分化レ
ベル 16）とする．d=θ=0 とし，格子と壁面を一致させる．
Figure 8 に，x=1.0 における乱流境界層速度分布と 0<x<2 における摩擦抗力係数を理論解と比較した結果を示す．速度分布・摩擦抗力係数分布ともに理論解と一致する．
よって物体適合格子と同程度まで格子細分化を行えば壁関数なしで定常乱流計算が可能であり，RANS 計算は妥当である．ただしこのときの格子数は54 万点であり，計算容量・計算時間の観点から格子数を抑える
必要性がある．
4

計算結果を理論解とすべりなし条件での結果と比較したもの
を Fig. 9 に示す．凡例のy+min は壁から一点目におけるy+を表す．
速度分布，摩擦抗力係数分布ともに最も粗い格子でも理論解・すCopyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2
Tab. 4 y+ at the first cell above the wallin each minimum grid size, where d is 0.5Δxmin
Minimumy+ at the first cell
grid size
above the wall
1.95×10-3
124
4.88×10-4
32.1
1.22×10-4
7.88
6.10×10-5
3.93

べりなし条件での結果と一致する．
物体適合格子の計算での一点目の y＋は 1 程度であるが，壁関数を導入することで y+=45 でも
乱流境界層を解像可能である．このとき摩擦速度を評価する点 F
での y＋は 128 で，粘性低層の外でも対数領域内であれば正しく
摩擦速度を評価できる．細分化のレベルは 6 下げられており，格子数を 1 万点以下まで抑えられる．

Tab. 5 y+ at the first cell above the wallin each minimum grid size, where d isΔxmin
Minimumy+ at the first cell
grid size
above the wall
1.95×10-3
123
4.88×10-4
31.4
1.22×10-4
7.79
6.10×10-5
3.86(a) Velocity profile at x=1.0(a) Velocity profile at x=1.0(b) Skin friction coefficientFig. 9 Mesh convergence study for flow simulation
around aligned flat plate using wall model, where d=0.0,
compared with theoretical value and result without wall model(2) 壁面を y 方向に動かした場合の計算
物体壁面の格子に対する位置を 0.5Δxmin，Δxmin と変化させて計算を行う．
最小格子幅は壁面と格子が一致した場合と同様に設定する．
各最小格子幅での壁から一点目におけるy+をd=0.5ΔxminとΔxmin それぞれの場合でまとめたものを Tab. 4，Tab. 5 に示す．
各格子幅について同様の計算を行った結果を Fig.10，Fig11 に示す．速度分布，摩擦抗力係数分布ともに最も粗い格子でも理論
解・すべりなし条件での結果と一致する．一点目の y+が 120 程度であっても境界層を解像可能である．
以上の θ=0(deg)の結果から，物体壁面と格子が同一方向の場合，格子が物体壁面から 1 セル分離れた場合でも乱流境界層を解像可能である．このとき壁面から一点目の y+は 120 程度でよく，格子数軽減の効果が現れている．

(b) Skin friction coefficientFig. 10 Mesh convergence study for flow simulation
around aligned flat plate using wall model, where d=0.5Δxmin,
compared with theoretical value and result without wall model

5

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2(a) Velocity profile at x=1.0(a) Velocity profile at x=1.0(b) Skin friction coefficientFig. 11 Mesh convergence study for flow simulation
around aligned flat plate using wall model, where d=Δxmin,
compared with theoretical value and result without wall model(b) Skin friction coefficientFig. 12 Mesh convergence study for flow simulation
around non-aligned flat plate using wall model,
compared with theoretical value and result without wall model(3) 壁面を格子に対して傾けた場合の計算
物体壁面を格子に対して 45 (deg) 傾けて同様の計算を行う．
各最小格子幅での壁から一点目における y+を Tab. 6 にまとめる．
計算結果を Fig. 12 に示す．一点目の y+が 29 程度になるように格子を細分化すると，乱流境界層を解像できる．

４．結論
複雑な物体形状周りにおける定常乱流の計算を可能とするために、直交格子に適した定常乱流計算手法を開発した．埋め込み境
界法に SA 乱流モデルから導かれた壁関数を組み合わせ，境界層解像に必要な格子数の増大を抑えた．
検証のために平板を用いた計算を行った．物体壁面と格子が同一方向の場合は壁から一点目のy+が120 程度で乱流境界層が解像可能であり，壁面を格子に対して傾けた場合でも一点目の y+が 30程度であれば乱流境界層を解像可能であった．

Tab. 6 y+ at the first cell above the wallin each minimum grid size, where θ is 45
Minimumy+ at the first cell
grid size
above the wall
1.95×10-3
107
4.88×10-4
28.9
1.22×10-4
7.08
6.10×10-5
3.54

謝辞
本研究は科学研究費助成事業，若手研究（B）（基金）
・23760767に基づき実施された．ここに感謝の意を表す．
参考文献(1) Delanaye, M., Aftosmis, M., Berger, M., Lin, Y., and Pulliam, T.,“Automatic Hybrid-Cartesian Grid Generation for High-Reynolds
Number Flows around Complex Geometries,” AIAA-99-0777, Jan,
1999.
6

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
講演番号 E10-2
(2)

(3)

(4)
(5)

(6)

(7)
(8)

(9)

(10)
(11)

(12)

(13)

(14)

Berger, Marsha J., and Michael J. Aftosmis, “Progress Towards a
Cartesian Cut-Cell Method for Viscous Compressible Flow,”
AIAA-2012-1301, 2012.
Lee, J.-D. and Ruffin, S. M., “Development of a Turbulent
Wall-Function based Viscous Cartesian-Grid Methodology,”
AIAA-2007-1326, 2007.
Capizzano, F., “A turbulent Wall Model for Immersed Boundary
Methods,” AIAA-2010-712. 2010.
今村太郎，神園仁志，高橋悠一，“航空機空力性能評価を目的とした直交格子法ベース流体ソルバーの開発，”平成 23年度航空宇宙空力班シンポジウム，2012 年 1 月.27 日．
神園仁志，今村太郎，“埋め込み境界法を用いた直交格子オ
イラーソルバーの改良，”第 44 回流体力学講演会・航空宇
宙数値シミュレーション技術シンポジウム 2012．
Spalart, P. and Allmaras, S., “A one equation turbulence model for
aerodynamic flows,” AIAA-92-0439, 1992.
Shima, E., "Parameter-free Simple Low-Dissipation
AUSM-Family Scheme for All Speeds," AIAA Journal Vol. 49,
No.8, August 2011, pp.1693-1709.
Nishikawa, H., “Two Way to Extend Diffusion Schemes toNavier-Stokes Schemes: Gradient Formula or Upwind Flux,”
AIAA -2-11-3044, 2011.
Nishikawa, H., “Robust and Accurate Viscous Discretization via
Upwind Scheme, I: Basic Principle,” Computers and Fluids, 2011.
Jameson, A., Turkel, E., “Implicit Schemes and LU
Decompositions,” Mathematics of Computation, 37, 156, pp.
385-397, 1981.
Chakravarthy, S., and Osher, S., " Numerical experiments with theOsher upwind scheme for the Euler equations, AIAAPaper82-0975,
1982.
Balaras, E., Benocci, C. and Piomelli, U., Two-layer approximateboundary conditions for large-eddy simulations. AIAA J. 34 (1996)
1111–1119.
W. Cabot and P. Moin, "Approximate wall boundary conditions in
the large-eddy simulation of high Reynolds number flow," Flow,
Turbul. Combust. 63, 269 (2000).

7

Copyright © 2013 by JSFM