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

LES のダイナミック壁面モデル：
剥離・遷移・再付着を伴う高レイノルズ数翼流れ
Dynamic wall modeling in large-eddy simulation:
Transitional separated flow over an airfoil at high Reynolds number

浅田 健吾，東大・院, 〒 252-5210 神奈川県相模原市中央区由野台 3-1-1, E-mail：asada@flab.isas.jaxa.jp
河合 宗司，ISAS/JAXA, 〒 252-5210 神奈川県相模原市中央区由野台 3-1-1, E-mail：kawai@flab.isas.jaxa.jp
藤井 孝藏，ISAS/JAXA, 〒 252-5210 神奈川県相模原市中央区由野台 3-1-1, E-mail：fujii@flab.isas.jaxa.jp

Kengo Asada, University of Tokyo, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210, Japan
Soshi Kawai, ISAS/JAXA, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210, Japan
Kozo Fujii, ISAS/JAXA, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210, Japan
Wall-modeled large-eddy simulation (LES) of subsonic flow over an A-Airfoil near stall condition at high Reynolds
number (chord based Rec = 2.1 × 106 ) is conducted to investigate the predictability of the wall-modeled LES forthe complex flow features including a laminar separation, turbulent transition and turbulent reattachment. By
incorporating the pressure and convective (non-equilibrium) effects and transition treatment in the wall model,
the wall-modeled LES well predicts the mean and turbulence statistics in the region of laminar separation, turbulent transition, turbulent reattachment, and attached turbulent boundary layer development. The present
LES with non-equilibrium wall-model well captures the physical processes of two-dimensional intermittent laminar separation, hairpin-like vortices, and breakdown of the laminar two-dimensional separation vortex to thethree-dimensional turbulent structures that allow boundary layer to re-attach. Comparisons between the nonequilibrium and equilibrium wall-models highlight the importance of including the non-equilibrium effects in themodel.

1.

組みが行われてきた．Chapman は翼型の例において外層域のみを解像した場合，必要な総格子点数は劇的に削減可能であり，Ntotal ∼ Rec0.4 であると見積もっている．

Introduction

航空機や宇宙機の設計を考える場合に乱流境界層は空力や空力加熱といった非常に重要な設計課題を含む分野の一つである．これらの工学的応用対象の多くは高レイ
ノルズ数流れであり乱流境界層における剥離や遷移，再付着，熱輸送といった現象を正確に予測するには壁面近傍の微細な渦構造が重要となる．
近年 Large-eddy simulation(LES) は乱流を解析する手法として非常に有効な手段であると注目されている．この手法は様々な分野において多くの成果を残しているが，計算資源の制約から工学的応用対象の多くが高レイノル
ズ数流れ（Re & 106 ）であるにもかかわらず，適用可能であるレイノルズ数の範囲は比較的低いレイノルズ数に限られてしまっている．これは乱流境界層が複数の長さスケールからなる現象である事が主な理由であり，特に境界層壁近傍 10%程度の内層域の渦スケールはレイノ
ルズ数の増加と共に急激に小さくなる．従って，高レイ
ノルズ数流れにおいて壁面近傍の微細な乱流構造を解像しようとすると，非常に膨大な計算資源が必要となる．
Chapman(1) は翼型を例に内層域まで解像する LES を行うのに必要な総格子点数 Ntotal を Ntotal ∼ Re1.8
c と見積もったが（Rec はコード長基準のレイノルズ数），これはDirect numerical simulation (DNS) に迫る計算コストであり，壁境界を持つ現実的な高レイノルズ数流れへ LESを適用しようとする際の大きな妨げとなっている．この
壁面近傍の問題（near-wall proble）を解決するためには内層域は直接解像せずにモデル化する必要性は不可避であり(2, 3) ，長い間多くの研究者によってモデル化の取り

1

LES における内層域のモデルは過去多くの提案がなされてきており，主要なものとして，2 つの手法が提案されている：1) 内層域で RANS タイプの乱流渦粘性にスイッチする手法 (例：DES や LES/RANS ハイブリッド
手法，review by Spalart(4) )，2) 壁面摩擦応力 τw を直接モデル化する手法 (review by Piomelli & Balaras(5) )．
手法 1(DES や LES/RANS ハイブリッド手法等) は複雑な剥離を伴う複雑形状にまで幅広く多くの研究者に用いられている．しかし Nikitin 等(6) が示したように，DESに代表されるような内層域で LES タイプの subgrid-scale
渦粘性 µt,LES を RANS タイプの乱流渦粘性 µt,RANS にス
イッチする手法を壁面モデルとして用いると，そのスイッ
チ部付近 (DES buffer layer とも呼ばれる) で “log-layer
mismatch” が発生し，最も基本的な付着乱流境界層においても壁面摩擦を 15%程度も小さく予測してしまうことが知られている．これはこのスイッチ部付近で非物理的に大きなスケールの渦が発生してしまい，解像している
乱れ成分が過小となってしまうのが原因である(7) ．従って現状，これらの手法はスイッチ部付近で，人工的にランダムな擾乱を付加する(7, 8, 9) か経験的に DES 中の長さスケールを急激に変化させる(10) ことで手法の改善を試みている．
一方，壁面摩擦を直接モデル化する手法 2 に関して，
Kawai & Larsson(11) は従来の壁面モデルを用いた LES
では壁からの数点は必ず数値的なエラーと subgrid-scale

Copyright c 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C09-2
Laminar boundary layer
development

Laminar separation bubble,
transition and turbulent reattchment

Wall-modeled LES meshturbulent separation

Inner-layer wall-model meshFig. 1: Airfoil geometry and resultant flow features on the suction
side: laminar separation, turbulent transition, turbulent reattachment, and turbulent separation.

Fig. 2: Wall-modeled LES and inner-layer wall-model meshes.

モデルのエラーを含んでいる事を圧縮性の平衡境界層方程式を解くことによって示し，当手法においても問題であった “log-layer mismatch” の主な原因である事を明らかにした．また，彼らは壁面境界付近での乱流の渦スケールの振る舞いを考慮する事によってこの “log-layer mismatch”
を除去し，正確な壁面剪断応力を予測できる手法を提案した．さらにこの手法を基に壁面モデル内で非定常な RANS
方程式を解く非平衡なダイナミック壁面モデルを提案し，高レイノルズ数における超音速乱流平板境界層や衝撃波/
乱流境界層干渉で生じる剥離流れ(12, 13) を精度良く予測できる事を示した．
本研究は最近の Kawai & Larsson (11, 12, 13) の研究成果をもとに彼らのモデルを遷移や剥離を伴う流れへ適用
範囲を拡張する．Kawai & Larsson の平衡及び非平衡モ
デルの両モデル(11, 12, 13) で新たに乱流遷移の取り扱いを提案し，高レイノルズ数（Rec = 2.1 × 106 ），亜音速の A-Airfoil 周りの流れに適用する．対象とする流れ場
は Fig. 1 に示す様に負圧側で層流剥離，乱流遷移，乱流
再付着，乱流剥離といった複雑な現象を含んでいる．得
られた結果は実験(14) 及び壁面近傍を解像する LES(15) のデータと比較し，壁面モデルの LES への適用可能性につて議論をおこなう．

2.

致しており，その位置での LES の瞬間値を壁面モデルの境界条件として与える．ここで y = hwm は壁面モデル格子の上境界と LES 格子が一致する位置（マッチングポイント）であり，hwm = y1（壁面から 1 点目）である必要はないが，過去の研究では例外無く LES の壁面から 1 点目の物理量を壁面モデルのインプットとして壁面モデル計
算を行なっていた．Kawai & Larsson(11) はこの hwm に関して，LES の壁面近傍の格子点では数値エラーが必然的に大きく，そのエラーの大きい格子点情報を使って壁面モ
デルを駆動しても正確な壁面摩擦や壁面熱流束を見積もることが出来ない事を示した．またその解決法として過去の全ての研究の慣例，壁面から 1 点目の情報を壁面モデル
へのインプットするのではなく，LES で正確に解像できている位置での物理情報を壁面モデルへのインプットとすることで，正確な壁面摩擦を LES にフィードバックできることを示した．本研究でも Kawai & Larsson(11) の研究に従い，LES の壁面から 4 点目 (hwm = y4 , h+
wm = 84)の物理量を壁面モデルの上境界インプットとして壁面モデルに与え，数値エラーの大きい領域をバイパスし，正確な LES データを用いて物理的に正しい壁面モデルを駆動させる．壁面モデルから LES にフィードバックされる
情報は壁面剪断応力と熱流束のみであり，一方，LES から壁面モデルに与えられる情報は上境界の物量のみである事を注記しておく．

Wall-modeled LES framework

本研究で用いる壁面モデルは壁面摩擦を直接モデル化する手法 (手法 2) をベースとしている．本手法のベースとなるアイディアは，レイノルズ数依存のほとんどない
境界層の 90%以上を占める外層域の乱流構造は LES として直接格子で解像するが，LES の格子で解像しないレ
イノルズ数依存の大きい内層域は壁面モデルでモデル化
し，非定常に変動する壁面摩擦応力 τw や壁面熱流束 qwを境界条件として LES にフィードバックする手法である．
従って，LES で用いる格子は Fig. 2 に示すように，外層域の乱流構造を解像する格子，すなわち境界層厚さでスケーリングされる格子を用い，粘性低層は解像しない(∆xi ≈ 0.05δ)．内層域を解く壁面モデルは Fig. 2 の右に示される格子の様に LES 格子の壁面近傍部分を取り出し，
壁面垂直方向に細かく格子を再分布させたものを別に用いる．壁面モデルの上境界は y = hwm で LES の格子と一2.1 Outer-layer LES equations
LES 計算では以下に示される空間的にフィルタリング
された圧縮性ナヴィ・エストークス方程式を外層域のみを解像するように用意された LES 用格子を用いて解く：

2

Copyright c 2012 by JSFM∂ρ∂+(ρuj ) = 0,
∂t
∂xj
∂
∂
∂p
∂τij
(ρui ) +(ρui uj ) +
=
,
∂t
∂xj
∂xi
∂xj
∂E
∂
∂
∂qj
+
[(E + p)uj ] =(τij ui ) −
,
∂t
∂xj
∂xj
∂xj
p
1
E=
+ ρu2 , p = ρRT,γ−1 2 k

(1)
(2)
(3)
(4)

各物理量は空間的にフィルタリングされており，ρ は密
度，u は速度，p は静圧，E は全エネルギー，T は温度，

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

γ(= 1.4) は比熱比，R は気体定数である．応力テンソル
τij 及び熱流束ベクトル qj は以下の様に表される．
2
τij = 2(µ + µt )Sij + [β − (µ + µt )]Skk ,
3

1 ∂ui
∂uj
Sij =
,
+
2 ∂xj
∂xi

 2
1
µt
µ
∂cs
qj =
+
,
γ − 1 P r P rt ∂xj

壁面モデル計算の壁面境界条件は滑りなしの断熱壁とし，この時の上境界の境界条件は (y = hwm = y4 ) での LESの基本変数（ρ, ui , p）の値をインプットとして設定する．
壁面モデルで得られた瞬間の壁面剪断応力と熱流束の値
は LES の壁面流束境界条件としてフィードバックされる．

(5)

(6)

ここで µ はサザーランドの式から計算される粘性係数，βは体積粘性率（本研究では β = 0），P r(= 0.72) はプラ
ントル数，Cs は音速である．subgrid-scale モデルとしてダイナミック Smagorinsky モデル(16) に Lilly の修正(17)を用い，乱流渦粘性 µt および乱流プラントル数 P rt を算出する．
LES の方程式を解く際の壁面境界条件として，壁面モ
デルからフィードバックされる流束による境界条件を用いる．すなわち壁面での壁面垂直方向の対流項流束や粘性による流束 (τij ui ) を 0 とし，壁面摩擦 τw や壁面熱流
束 qw は，壁面モデル計算で得られる値を使う．この流束境界条件の意味するところは，LES は内層域を格子で解像しないが，壁面上で滑り無し ui = 0 という事実は変わらず，ただ内層域での物理量の勾配を正確に計算できな
いということを意味している．言い換えると，滑り無し
壁条件による壁面上での物理量を用いて，LES 壁面付近の格子点における物理量の勾配を正確に見積もる事が困難であるということを意味している．壁面モデルを用いた LES が内層域を格子で解像しないという事実は，計算に用いる数値スキームにもいくつかの変更を加える必要が生ずる．簡潔に言うと，LES の格子がない内層域よりも上の格子点と下の格子点間 (すなわち壁から 1 点目の格子 j = 1 と壁 j = 0 との間) の差分は定義できず，そのまま j = 1 と j = 0 の間で計算する差分は不正確になる．本研究では，粘性流束の計算に用いる壁面垂直方向の差分として，j = 1 の点では完全な 2 次精度片側差分，
j = 2 では 2 次精度中心差分，その他の点では 3 重対角の6 次精度コンパクト差分法を用いた．以上の差分の取り扱いにより，LES の方程式は壁面上での速度や密度，温度を陽に定義する事無く計算を進める事が出来る．

2.2

2.2.1 Equilibrium wall model full-RANS 方程式に平衡流近似を用いた平衡境界層方程式は壁面モデル領域（y = 0 から y = hwm ）で以下の 2 つのカップルした
常微分方程式に単純化される(11) ．


dUk
d
(µ + µt,wm )
= 0,
dy
dy

dUk
d
(µ + µt,wm ) Uk
(7)
dy
dy



µt,wm dT
µ
+cp
+
= 0,
P r P rt,wm dy
ここで Uk は壁面に沿った速度，T は温度である（従って壁面垂直方向にのみ方程式を解けばよい）．これらの方
程式は流れ方向の運動量保存の式及びエネルギー保存の式に一般的な平衡平板境界層の近似を用いることで導かれる(18) ．この単純なモデルでは圧力は壁面垂直方向に一定の値が想定され，上境界 y = hwm と同じ，すなわち外層の LES 計算から受け取った値と同じ値が設定されることになる．また，本計算では壁面剪断応力と熱流束の値が収束するまで上記の方程式を数値的に解いている．熱輸送のない非圧縮の極限ではこの単純な壁面モデルは対数領域において広く知られている log-law を示す．
乱流渦粘性 µt,wm は mixing-length 渦粘性モデルに vanDriest ダンピングを用いて評価する．
rτw
D,µt,wm = κmod ρy
ρ
(8)

2
D = 1 − exp(−y + /A+ ) ,p
ここで τw /ρ は境界層内の局所的な瞬間密度と壁面剪
断応力で見積られる速度スケール，y + = ρwall uτ y/µwall
は壁からの距離の viscous ユニットで，A+ = 17 である．
モデリングパラメータとして κmod は通常の von Kármán定数と同じ κmod = 0.41 を用い，乱流プラントル数は定数として P rt,wm = 0.9 を用いる．

Inner-layer wall-model equations

本研究では RANS の概念を基にした 2 つの異なる壁面モデルを比較する．1 つは平衡流近似を行った平衡(EQBL) モデル(11) ，もう 1 つは近似なしの full-RANS を解く非平衡 (NonEQBL) モデル(11, 12, 13) である．EQBL,
NonEQBL のいずれのモデルを用いた LES も高レイノル
ズ数下，圧力勾配なしの付着平板乱流境界層において壁面応力及び乱流統計量を精度よく予測可能である事が確認されている(11, 12, 13) ．本研究ではこれら 2 つのモデル
を乱流遷移や剥離の伴う高レイノルズ数翼流れに適用し，壁面モデルの LES への適用可能性について議論を行う．

2.2.2 Non-equilibrium wall model 非平衡壁面モ
デルでは粘性拡散項に加えて，対流項，圧力項の効果を含む full-RANS 方程式を別途用意した RANS 用格子を用いて内層域のみで計算する．壁面モデル計算領域の壁面
境界条件は滑り無し断熱壁条件を用い，y = hwm = y4 における上境界条件には，時々刻々と変化する LES で得られる瞬間の物理量が用いられる．壁面モデルの計算では，
壁面垂直方向に格子をストレッチさせているため，RANSという観点から乱流境界層を解像する．従って，各時間
ステップ毎に，時間精度を持った非平衡壁面モデルの計

3

Copyright c 2012 by JSFM

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

算結果から壁面摩擦や壁面熱流束を算出でき，LES の壁
面境界条件としてフィードバックすることが出来る．
RANS の乱流渦粘性 µt,wm 及び乱流プラントル数
P rt,wm を見積る方法として Kawai & Larsson によって提案された方法(12, 13) を用いる．RANS の乱流渦粘性 µt,wm
は (8) 式で示される mixing-length 渦粘性モデルを基に求められるが，壁面モデルの計算内で µt,wm , P rt,wm の値を
ダイナミックに修正する．ここで用いるダイナミック非平
衡壁面モデルのキーは 1) マッチングポイント（y = hwm ）において全剪断応力を LES と壁面モデルとで一致させる
事と 2) 格子で解像している乱流スケールと解像できていないスケールが壁面垂直方向に変化しおり，その影響を近似的にモデルに取り込む事である．
1 つめのマッチング条件に関して，LES における対数領域（非粘性領域 30 . y + . 0.1δ + ）の全剪断応力バラ
ンスは運動量保存則から簡単に求める事ができ，以下の様に表される．
′′ v ′′ = τ ,µLES ∂y ue − ρu]w(9)

この式は平均速度勾配 ∂y u
e は全剪断応力 τw と解像される剪断応力，平均の subgrid 渦粘性 µLES のバランスで決まる事を示している．境界層内層の RANS 計算において壁面モデルの上境界の速度及び温度は LES と一致しているため，それらの値から計算される解像される応力と熱流束の値も LES と一致している．従って，マッチング
ポイントにおいて全応力と全熱流束を一致させるには解像されていない応力及び熱流束を一致させる必要がある．
Wang & Moin(19) が提案するマッチングのアイディアと同様に，本モデルでも位置 y = hwm において全せん断応力を近似的に一致させるために，乱流渦粘性と乱流熱伝導率（Cp µt /P rt ）を LES と RANS の計算で一致させ
cr t とる．マッチング条件にあう κmod と P rt を κb 及び P
すると
hµt,LES i
q
κ
b=
,
(10)
hρy τρw Diの様にダイナミックに κb を決定する事で乱流渦粘性をマッ
チングの位置において一致させる事ができる．hi はスパ
ン方向の平均を表し，従って κb は流れ方向及び時間と共
cr t,RANS = hP rt,LES i とすに変化する．また，熱流束は Pる事で RANS と LES の計算の間で一致させる．
Wang & Moin(19) もこの様にマッチング条件に合うκmod と P rt をダイナミックに決定し，マッチングの位置で求めた κ
b を境界層の内層域内全てで用いているが， 2
つめのキーとなる項目として述べた様に格子で解像しているものとしていないものの寄与分を考慮する必要がある．そこで本計算で用いるモデルでは κmod の値を壁面モデル計算領域内壁面垂直方向で変化させる．まず一般に知られているように，log-layer における支配的な乱流スケールは壁面からの距離 y に比例する(18) ．そこで壁面に平行な方向の log-layer における支配的な乱流スケー
ルを Lk と定義すると，Lk = Ck y と書くことが出来る(ここで Ck は内層域の物理で決まる定数)．添字の k は
壁面に平行な方向の長さスケールを表しており，例えば
流れ方向は L1 ，スパン方向は L3 となる．また代表的な壁面平行方向の格子幅 ∆k と乱流スケール Lk の比を取ると，Lk /∆k = Ck y/∆k となる．ここで代表的格子幅は∆k = max(∆x, ∆z) と定義され，∆k は格子がサポートできる最も小さなスケールと考えることが出来る．すなわ
ち Lk /∆k がある値 α (Lk に対し格子が何点あるかを表す
定数．壁面モデルで用いる数値スキームに依存する) より′′ v ′′
も小さければ，格子で解像できているせん断応力 −ρu]は無視できるほど小さいと考えられる．すなわちここでモデル定数は通常の RANS 解析で用いられる κmod = 0.41，
P rt = 0.9 を用いるべきである．逆に Lk /∆k > α の領域
では，格子で解像している Reynolds 応力が存在するためモデル化による応力 (κmod ) を位置 y = hwm のマッチ
ングの値 κ
b に向けて小さくしていくべきである．ここで，壁面モデルで用いる計算格子は RANS という観点から乱+流境界層を解像しているため (ywall< 1)，∆y は L2 を十分に解像できるほど細かく，∆y ≪ ∆k である．本研究では線形的なダンピング関数 K を用いて κmod と P rt を
以下のように定義する．
κmod= 0.41K + κb(1 − K) ,c
= 0.9K + P r t,wm (1 − K) ,

P rt
ここで

K = min
Lk,top /∆k − Lk /∆k, 1 ,Lk,top /∆k − α
Lk,top = Ck ytop = Ck hwm

(11)
(12)

(13)

もしくは，少し式変形をして
ytop − y
K = min, 1 ,ytop − ycrit
αycrit =
∆k = α′ ∆k
Ck

(14)

と表せる．後者の定式化 (式 14) からも明らかなように，
K は流れ場依存の関数ではなく，事前に決定することが出来る壁からの距離 y とパラメータ α′ = α/Ck のみの関数である．またここで α′ は物理から決まる定数 Ck と壁面モデルで用いる計算スキームで決まる定数 α の関数である．図は，関数 K が異なる α′ でどのような振る舞いを示すか表した図である．過去の研究から(12, 13) 本研究
では α′ = 0.48 とする．

2.2.3 Treatment of transition in the wall models 本研究で対象とする流れ場は翼前縁付近で乱流遷移を含んでいる．通常の LES では乱流遷移を直接解像するが，壁面モデルを用いた LES では境界層内層域の物理現象は解像しない．そのために乱流遷移を特別に扱う必要がある．そこで本研究では乱流への遷移点より上流の領
4

Copyright c 2012 by JSFM

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

ルの LES への適用可能性について議論を行う．また非平衡 (NonEQBL) モデル(11) と平衡 (EQBL) モデル(12, 13)を比較し，逆圧力勾配下で剥離を生じる流れ場の予測において非平衡効果がいかに重要性であるかを示す．

1

0.8

K

0.6

0.4

3.1

A-Airfoil 周りの計算格子は C-type のトポロジーを採用
し，コード長方向 (x) 及び垂直方向（y ）に 20c（20 chord
length），スパン方向（z ）に 1.7%c の計算領域をとる．
スパン方向の相関長さに関して参考となる実験データはないが，スパン方向に 1.2%c の計算領域を用いた内層域を解像する LES の先行研究(14) において x/c ≈ 0.8 以降
で短いスパン領域を用いた影響が見られると述べられているもののそれ以外の領域では実験と良い一致を示している．そこで本研究では Mary & Sagaut らの研究(14) よりもスパン方向にわずかに広い計算領域を用いる．

0.2

0

0

0.01

0.02

0.03

0.04

0.05

0.06

y/δ
Fig. 3: Linear blending function K with α′ = 0.48.

域で壁面モデル計算の乱流渦粘性を µt,wm = 0 とする事
で乱流遷移を取り扱う．Mary & Sagaut(14) の境界層内層域を解像する LES によると x/c ≈ 0.14 で乱流遷移を起こすと言及されている．そこで，本研究では遷移点を
x/c = 0.14 としている．

2.3

LES に用いる格子解像度は外層域を解像する様に設定し，従って格子幅は局所的な境界層厚さ δ によってスケーリングされる．本研究では Fig. 4 に示す様な 2813 × 65 ×
149（負圧側に 2263 点，正圧側に 350 点，後流に 200 点）の格子を用いた．壁面モデルを用いた LES では格子解像
度は局所的な境界層厚さ δ によって決められるので，流れ方向の格子は前縁付近に集められその後境界層の発達に従って滑らかにストレッチングしている．本研究と同じ内層域のモデル化と空間離散化を用いて Kawai & Larsson
が行った平板乱流境界層の格子収束に関する研究(11) に基づき，発達する境界層に従って局所的な境界層厚さ δ を流れ方向，壁面垂直方向に 15 − 20 点以上で解像するように計算格子を作成した（δ/∆x & 15 − 20, δ/∆y & 15 − 20）．
またスパン方向には局所的な境界層厚さ δ に対して翼前縁近傍で少なくとも 7 点，x/c ≈ 0.3 以降では 15 − 20 点程度の格子点が入るように解像度を設定した．

Numerical schemes

空間離散化は保存形で，内部点には 6 次精度コンパク
ト差分法(20) を用い，境界の 1 点目には 2 次精度片側差分，2 点目には 2 次精度中心差分を用いている．エイリ
アシングエラー等を許容するため，8 次精度の low-pass
フィルター(20, 21) を用いる（αf = 0.45）．
壁面モデルを用いる本 LES は外層スケールのみを直接格子で解像し，内層域下部や粘性底層は格子で解像しな
いため，壁から 1 点目の格子点位置が本研究では y1+ ≈ 21となる．壁面近傍でこのような荒い格子を用いることから，壁面モデルを用いた LES は通常の LES と比べ 1 桁程
度大きい時間刻み幅を取ることができるというメリットがある．時間積分法として，LES の計算では 3 次精度 Total
variation diminishing Runge-Kutta 法(22) を，壁面モデ
ル計算には壁面近傍での CFL 条件を緩和するため 2 次
精度陰解法(23, 24) に内部反復法を組み合わせたものを用∗
い，双方で時間刻み幅は同じにした (∆t = ∆t∗ U∞/c∗ =4
3.3 × 10 ，ここで上付きの ∗ は有次元量を表す)．

3.

Computational grid and time-step size

壁面モデルを用いた LES では境界層内層域の物理現象を解像しないために格子解像度の評価指標で重要なのは局所的な境界層厚さ δ でスケーリングされる格子解像度であり，viscous wall unit(内層域の渦スケール) は意
味をなさない事に注意されたい．しかしながら viscouswall unit を基準とした格子解像度は通常の境界層内層を解像する LES との比較において有用である．文献(14)の最大壁面摩擦応力から計算した viscous wall unit を基準とした場合の本研究の LES 用格子の解像度は ∆x+ =+21 − 400, ∆z + = 56, ∆yw= 21, h+wm ≈ 84 である．

Results

本研究では亜音速 (M∞ = 0.15)，高レイノルズ数流れ（コード長基準のレイノルズ数 Rec = 2.1 × 106 ）において，失速迎角付近（α = 13.3）の Aerospatiale AAirfoil 周り流れの LES に壁面モデルを適用した．今回対象にする流れ場は幾つかの複雑な物理現象を含んでいる．例えば負圧側において翼前縁付近で層流剥離，乱流
遷移，乱流再付着し，翼後縁付近で再度乱流剥離する．本
研究の気流条件は過去に ONERA によって行われた豊富な実験結果(15) を元にしており，同様の条件で従来の乱流
境界層内層まで解像する LES が Brite-Euram LESFOIL
project(14, 25) の下で行われている．本研究で行った壁面
モデルを用いた LES の結果とこれらを比較し，壁面モデ

本研究では LES と壁面モデルの両計算で同じ時間刻み∗
幅 ∆t = ∆t∗ U∞/c∗ = 3.3 × 104 （上付きの ∗ は有次元数を表す）を用いており，LES における最大クーラン数
は 0.6 である．壁面モデルを用いた LES では従来の LES
よりも荒い格子を用いるために大きな時間刻み幅を採用することが可能となり，1 桁程度従来の境界層内層域を解像する LES よりも計算コストを削減できる事を強調しておく．よって壁面モデルを用いた LES では２つの理由で計算コストを削減する事ができる（格子点数の削減及

5

Copyright c 2012 by JSFM

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

-4

Cp

-3

-2

-1

0

Fig. 4: Computational grid with close-up views near transition
and trailing-edge regions overlaid with instantaneous streamwise(x) velocity contours for the LES mesh 2813 × 65 × 149. Everytenth grid point is shown.

1
0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

x/c
0.024

び時間刻み幅の制限緩和）．また，Chapman の解析(1) を参考にするとレイノルズ数が増せば増すほど劇的に計算コストを削減できる事になる．

Cf

0.016

3.2

0.008

Comparisons of statistics

Figure 5 に NonEQBL 及び EQBL モデルを用いた LES
2の翼面上の平均圧力係数 Cp = 2(p−p∞ )/ρU∞
分布，平均2
壁面応力係数 Cf = 2τw /ρU∞ 分布を示す．いずれのモデ
ルも翼後縁付近で壁面応力を大きく見積もっており，実験値に見られる様な平らな圧力分布は予測できていない．しかしながら，再付着点付近での急激な壁面応力の上昇の様な大域的な積分量は精度よく予測できている．NonEQBLモデルと EQBL モデルの違いは主に乱流再付着より上流（x/c . 0.15）の領域に見られ，特に NonEQBL を用いた
LES では周期的に層流剥離剪断層から放出される渦を捉える事ができている（詳しくは 3.3 節にて議論する）．また NonEQBL モデルでは x/c ≈ 0.1 付近で剥離している
事を示す平らな圧力分布や負の壁面応力分布を予測する事ができているが，EQBL モデルではこの様な分布は見られず，流れが付着したままであるのがわかる．EQBL
モデルでは粘性拡散項の効果のみ考慮しているのに対して，NonEQBL モデルでは圧力項及び対流項の効果も取り扱っているために逆圧力勾配によって生じる層流剥離を NoEQBL モデルだけが捉えられたという事実は理に適っていると言える．
壁面モデルを用いた LES と実験(15) 及び境界層内層を解像する LES(14) のコード長方向速度プロファイルを壁
面垂直方向に沿って x/c = 0.1 − 0.99 の位置で比較したものを Fig. 6 に示す．NonEQBL モデルを用いた LES では翼前縁の x/c ≈ 0.1 付近に生じる層流剥離泡から乱流遷移，乱流再付着，発達する乱流境界層領域の x/c ≈ 0.5付近まで精度よく速度プロファイルを予測することができている．一方 EQBL を用いた LES では層流剥離は存在せず，乱流遷移，乱流再付着，乱流境界層発達の初期段
階で速度プロファイルを大きく見積もっている．この主な原因は EQBL モデルが翼前縁付近の逆圧力勾配の効果を取り扱え無いために，適切な壁面応力が LES にフィードバックされないためである．乱流再付着してから十分

0-0.008
0

0.2

0.4

x/c
Fig. 5: Mean pressure and skin friction coefficients along the
airfoil compared to the corresponding experiments (15) (circles),and wall-resolved LES (14) (crosses). NonEQBL (solid line, blue);
EQBL (dashed line, red).

に乱流境界層が発達した x/c ≈ 0.5 以降の領域では壁面モデルの種類に対する平均速度分布の感度は低くなって
おり，NonEQBL モデル，EQBL モデルを用いた LES 共によく似た乱流境界層の発達を示している．x/c ≈ 0.7 より下流では壁面モデルを用いた LES のコード長方向平均
速度プロファイルは徐々に実験値から外れていき値を大きく見積もってしまっている．この様に x/c ≈ 0.7 以降
で予測精度が落ちるのは本研究では比較的短いスパン長（1.7 % c）を採用した事が一因として考えられる．文献(14)
中でも不十分なスパン長が与える影響について言及さ
れている．
√
√
次 に レ イ ノ ル ズ 応 力 分 布（ u′ u′ /U∞ ,v ′ v ′ /U∞ ,′ ′2−u v /U∞ ）を Figs. 7-9 に示す．NonEQBL モデルを用いた LES では速度プロファイルと同様に翼前縁の x/c ≈ 0.1付近に生じる層流剥離泡から乱流遷移，乱流再付着，発達する乱流境界層領域の x/c ≈ 0.7 付近まで乱流統計量を精度よく予測することができている．しかしながら，非平衡効果を考慮していない EQBL モデルを用いた LES では層流剥離及び乱流再付着付近においてレイノルズ応力を小さく見積もってしまっている．平均速度分布と同様に乱流再付着し，十分に乱流境界層が発達した x/c ≈ 0.5より下流の領域では，壁面モデルの種類による影響はさほど見られなくなる．x/c ≈ 0.835 以降では両モデル共

6

Copyright c 2012 by JSFM

0.012

0.012

0.01

0.01

(y − yw )/c(y − yw )/c

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

0.008
0.006
0.004
0.002

0.008
0.006
0.004
0.002

0

0
0

1.4

2.8

4.2

5.6

7

0

0.3

0.09

0.09

0.08

0.08

0.07

0.07

(y − yw )/c(y − yw )/c

U/U∞

0.06
0.05
0.04
0.03

1.2

1.5

0.6

0.9

1.2

1.5

0.05
0.04
0.03
0.02

0.01

0.01
0
0

1.4

2.8

4.2

5.6

7

0

U/U∞
Fig. 6: Mean streamwise velocity (U/U∞ ) profiles as a functionof wall-normal distance at x/c = 0.1, 0.15, 0.2, 0.3, 0.5 on top figure and x/c = 0.7, 0.825, 0.87, 0.93, 0.99 on bottom. Each plot isseparated by a horizontal offset of 1.4. Symbols and lines as inFig. 5.

0.3√
u′ u′ /U∞

Fig.
7:
Resolved streamwise Reynolds normal stress√( u′ u′ /U∞ ) as a function of wall-normal distance at x/c =
0.1, 0.15, 0.2, 0.3, 0.5 on left and x/c = 0.7, 0.825, 0.87, 0.93, 0.99
on right. Each plot is separated by a horizontal offset of 0.3.
Symbols and lines as in Fig. 5.

に予測精度は落ち，境界層上方でレイノルズ応力を小さく見積もっている．この予測精度低下の原因としてスパン長が短いために境界層上方の大きな乱流構造が非物理的な 2 次元的な構造に拘束された事が考えられる．

3.3

0.9

0.06

0.02

0

0.6

√
u′ u′ /U∞

Instantaneous flowfield

Figure 10 に NoneEQBL, EQBL 両ケースのコード長
方向 (x) 速度の瞬間場を示す．上に流れ場全体の側面図，
下に層流剥離後，乱流遷移し再付着する付近の拡大図（上：側面図，下：上面図）を示している．上面図は y = hwm =y4 (h+
wm ≈ 84) の位置での値を示している．NonEQBL モ
デルを用いた LES では青色で表される剥離領域からスパン方向の 2 次元的な渦が間欠的に放出された後，3 次元的な細かな渦構造を持った乱流へ遷移し，再付着する一連の物理現象を捉えることができている．一方で EQBL モデルを用いた LES では層流剥離領域が見られず，NonEQBL
モデルのケースで見られた様な 2 次元的な渦が崩壊して乱流へと遷移する過程なしに，突然乱流へと遷移している．
翼前縁付近の層流剥離から生じた 2 次元的な渦が 3 次
元的な渦へと崩壊する様子を速度勾配テンソルの第 2 不
変量 Q で Fig. 11 に示す．Q 値の等値面は渦構造の高さがわかるようにコード長方向速度で色付けしている．
NonEQBL モデルのケースでは剪断層の不安定性から放7出されたスパン方向の 2 次元的な渦がヘアピン状の渦へと崩壊していく様子が見てとれる．ヘアピン渦の崩壊後，下流ではより細かくランダムな渦構造を持つ乱流へと発達し，流れの再付着を促進している．この一連の現象（層流剥離剪断層から放出された 2 次元的な渦がヘアピン渦をへて崩壊していく現象）は境界層の外層域で起こっているために, 壁面モデルを用いた LES で解像可能であったと考えられる．この様に NonEQBL モデルでは非平衡効果と遷移を適切に取り扱う事で，層流剥離から乱流再付着までの現象を自動的に再現すると共に，遷移後の適切な壁面応力を予測できている．一方，EQBL モデルを用いた LES は Fig. 5 に示す様に翼前縁付近で壁面応力を大きく見積もってしまっているために流れは付着したままとなっている. また x/c = 0.14 の位置から乱流粘性µt,wm を付加し始め，壁面応力 τw (LES に流束境界条件
として与えられる) が急に増加するが，それによって生成される始める渦の位置は非常に壁面に近い位置となっている. このように x/c & 0.5 の位置での統計量はほぼ同じであるのに NonEQBL モデルと EQBL モデルとでは乱流への発展過程が全く異なっている事がわかる．

Copyright c 2012 by JSFM

0.012

0.012

0.01

0.01

(y − yw )/c(y − yw )/c

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

0.008
0.006
0.004
0.002

0.008
0.006
0.004
0.002

0

0
0.3

0.6

0.9

√
v ′ v ′ /U∞

1.2

1.5

0

0.09

0.09

0.08

0.08

0.07

0.07

(y − yw )/c(y − yw )/c

0

0.06
0.05
0.04
0.03

0.028

0.042

0.056

0.07

0.056

0.07

2
−u′ v ′ /U∞

0.06
0.05
0.04
0.03

0.02

0.02

0.01

0.01

0

0.014

0
0

0.3

0.6

0.9

√
v ′ v ′ /U∞

1.2

1.5

0

0.014

0.028

0.042

′ ′

2
−u v /U∞2 ) as a function
Fig. 9: Resolved Reynolds shear stress (−u′ v′ /U∞of wall-normal distance at x/c = 0.1, 0.15, 0.2, 0.3, 0.5 on left andx/c = 0.7, 0.825, 0.87, 0.93, 0.99 on right. Each plot is separated
by a horizontal offset of 0.014. Symbols and lines as in Fig. 5.

Fig.
8:
Resolved wall-normal Reynolds normal stress√( v′ v′ /U∞ ) as a function of wall-normal distance at x/c =
0.1, 0.15, 0.2, 0.3, 0.5 on left and x/c = 0.7, 0.825, 0.87, 0.93, 0.99
on right. Each plot is separated by a horizontal offset of 0.3.Symbols and lines as in Fig. 5.

Acknowledgement
4.

Conclusions

高レイノルズ数付着乱流境界層を正確に予測可能な壁面モデルを高レイノルズ数下で剥離，遷移，再付着を伴う A-airfoil 周り流れの LES に適用し，複雑な流れへの適用可能性について議論を行なった．本研究では最近の研究成果に基づいた平衡モデル及び非平衡モデルを剥離
及び遷移を考慮したものへと拡張した．
対流項及び圧力項の効果を考慮した非平衡モデルを用いた LES では x/c ≈ 0.7 より上流の層流剥離，乱流遷移，
乱流再付着，発達する乱流境界層の領域で平均量，乱流統
計量共に精度よく予測することが可能であった．さらに，
非平衡モデルは層流剥離から 2 次元的なスパン方向の渦が放出され，逆圧力勾配による不安定性で崩壊し微細な乱流渦構造へと遷移しながら再付着するという一連の物
理現象を捉えることに成功している．一方で，非平衡の効果を考慮していない平衡モデルではこれら一連の特徴的な物理現象を捉えることはできず，非平衡モデルの重要性を示唆している．また，本研究の Rec = 2.1 × 106 における翼流れの LES において，壁面モデルを用いた LES
は境界層の内層領域まで解像する LES に比べて（格子点数の削減及び最小格子幅に依存する時間刻みの制限緩和
により）実に 100 倍近いオーダーで計算コストを削減できることを強調しておく．

8

本研究は JAXA International Top Young Fellow-ship
及び科学研究費（若手 (B)No. 24760570），東京大学
Global COE Program のによる支援を受けた．本研究の計算は JAXA の Japan Aerospace Exploration Agency
Supercomputer System (JSS) を使用した．また本論文
著者・浅田は Visiting scholar としてスタンフォード大学の Lele 教授に受け入れて頂き，有意義なコメントを多く頂いた．ここに謝意を表す．

References(1) D. R. Chapman. Computational aerodynamics development and outlook. AIAA Journal, Vol. 17,
No. 12, pp. 1293–1313, December 1979.
(2) J. W. Deardorff. A numerical study of threedimensional turbulent channel flow at large
Reynolds number. Journal of Fluid Mechanics,
Vol. 41, pp. 453–480, 1970.
(3) U. Schumann. Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. Journal of Computational Physics,
Vol. 18, No. 4, pp. 376–404, 1975.
Copyright c 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C09-2(4) P. R. Spalart. Detached-eddy simulation. Annual
Review of Fluid Mechanics, Vol. 41, pp. 181–202,
2009.
(5) U. Piomelli and E. Balaras. Wall-layer models forlarge-eddy simulations. Annual Review of Fluid
Mechanics, Vol. 34, pp. 349–374, 2002.

0.12

0.14

0.16

0.18

0.20

0.22

(6) N. V. Nikitin, F. Nicoud, B. Wasistho, K. D.
Squires, and P. R. Spalart. An approach to wallmodeling in large-eddy simulations. Physics of Fluids, Vol. 12, No. 7, pp. 1629–1632, July 2000.

0.24

x/c

(7) U. Piomelli, E. Balaras, H. Pasinato, K. D. Squires,and P. R. Spalart. The inner-outer layer interface in large-eddy simulations with wall-layer models. International Journal of Heat and Fluid Flow,
Vol. 24, No. 4, pp. 538–550, August 2003.

(a) LES with NonEQBL wall model(8) J. E. Brasseur and T. Wei. Designing large-eddy
simulation of the turbulent boundary layer to capture law-of-the-wall scaling. Physics of Fluids,
Vol. 22, No. 2, p. 021303, February 2010.

0.12

0.14

0.16

0.18

0.20

0.22

(9) J. Larsson, F. S. Lien, and E. Yee. Feedbackcontrolled forcing in hybrid LES/RANS. International Journal of Computational Fluid Dynamics,
Vol. 20, No. 10, pp. 687–699, 2006.

0.24

x/c

(b) LES with EQBL wall model(10) M. L. Shur, P. R. Spalart, M. Kh. Strelets, andA. K. Travin. A hybrid RANS-LES approachwith delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow,
Vol. 29, No. 6, pp. 1638–1649, December 2008.

Fig. 10: Instantaneous streamwise velocity in a side-view plane(top), and in a close-up side-view and top-view plane near transition region (bottom) at y = hwm = y4 (h+wm = 84) obtained by
LES with NonEQBL and EQBL wall models. Contours range of
0.0 (blue) ≤ u/cs ≤ 0.4 (red).

(11) S. Kawai and J. Larsson.
Wall-modelingin large eddy simulation: Length scales, grid
resolution, and accuracy.
Physics of Fluids, Vol. 24, No. 1, p. 015105, January 2012.
http://dx.doi.org/10.1063/1.3678331.
(12) S. Kawai and J. Larsson. A non-equilibrium wallmodel for les of shock/boundary layer interaction
at high reynolds numbers. AIAA Paper 2012–2226,
AIAA, June 2012.

(a) LES with NonEQBL wall model(13) S. Kawai and J. Larsson. Dynamic non-equilibrium
wall-modeling for large-eddy simulation at highreynolds numbers. Physics of Fluids, 2012. underreview.
(14) I. Mary and P. Sagaut. Large eddy simulation offlow around an airfoil near stall. AIAA Journal,
Vol. 40, No. 6, pp. 1139–1145, July 2002.

(b) LES with EQBL wall modelFig. 11: Instantaneous snapshots of vortex structures in thelaminar-turbulent transition region near leading edge obtained by
LES with NonEQBL and EQBL wall models. Isosurfaces of thesecond invariant of velocity gradient tensor Q colored by streamwise velocity. Every fifth grid point on the airfoil surface is shown.

(15) C. Gleyzes. Opíon décrochage–résultats des essaisá la souferie f2. Technical Rept. RT-OA 19/5025,
ONERA, June 1978.
9

Copyright c 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C09-2(16) P. Moin, K. Squires, W. Cabot, and S. Lee. A
dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A,
Vol. 3, No. 11, pp. 2746–2757, November 1991.
(17) D. K. Lilly. A proposed modification of theGermano subgrid-scale closure method. Physics of
Fluids A, Vol. 4, No. 3, pp. 633–635, 1992.
(18) S. B. Pope. Turbulent Flows. Cambridge University Press, 2000.
(19) M. Wang and P. Moin. Dynamic wall modeling forlarge-eddy simulation of complex turbulent flows.
Physics of Fluids, Vol. 14, No. 7, pp. 2043–2051,
July 2002.
(20) S. K. Lele. Compact finite difference schemes with
spectral-like resolution. Journal of Computational
Physics, Vol. 103, No. 1, pp. 16–42, 1992.
(21) D. V. Gaitonde and M. R. Visbal. Padé-type
higher-order boundary filters for the Navier-Stokes
equations. AIAA Journal, Vol. 38, No. 11, pp.
2103–2112, November 2000.
(22) S. Gottlieb and C. W. Shu. Total variation diminishing Runge-Kutta scheme. Mathematics ofComputation, Vol. 67, No. 221, pp. 73–85, January
1998.
(23) S. Obayashi, K. Fujii, and S. Gavali. Navier-Stokes
simulation of wind-tunnel flow using LU-ADI factorization algorithm. NASA-TM-100042, NASA,
February 1988.
(24) N. Iizuka. Study of Mach Number Effect on the Dynamic Stability of a Blunt Re-entry Capsule. PhD
thesis, University of Tokyo, February 2006.
(25) C. P. Mellen, J. Frohlich, and W. Rodi. Lessons
from lesfoil project on large-eddy simulation of flow
around an airfoil. AIAA Journal, Vol. 41, No. 4,
pp. 573–581, April 2003.

10

Copyright c 2012 by JSFM