JSME TED Newsletter, No.87, 2019

TED Plaza[連載企画] 乱流熱伝達モデル入門（第２回）
服部博文

名古屋工業大学 課長
技術部情報解析技術課
hattori@nitech.ac.jp
４．渦拡散（渦粘性）係数のモデリング
レイノルズ応力と乱流熱流束は，勾配拡散型モデルを用いて次式のようにモデル化されるが，このモデリングで定義される運動量もしくは熱の渦拡散係数は場の変数であるため，さらにモデル化が必要であることは本企画第１回で述べた．
̅𝑗
̅𝑖 𝜕𝑈
𝜕𝑈
1
+
) − 𝑢𝑘 𝑢𝑘 𝛿𝑖𝑗𝜕𝑥𝑗 𝜕𝑥𝑖3𝜕𝑇̅−𝑢𝑖 𝑢𝑗 = 𝜈𝑡 (−𝑢𝑖 𝑡 = 𝛼𝑡𝜕𝑥𝑖

ここで，これら渦拡散係数は，本企画第１回の３・２節で述べたように，乱れの速度𝑉×乱れの空
間スケール𝐿（もしくは，乱れの時間スケールτを用いて，𝐿 = 𝑉τとして空間スケールを表現する）
で表現される．
𝜈𝑡 = 𝐶𝑉 ∙ 𝐿
𝛼𝑡 = 𝐶𝑉 ∙ 𝑉τ
勾配拡散型モデルを用いた乱流モデルの分類は，これら乱れのスケールに対して，新たに解かれる方程式の数で決められるが，ここでは，その渦拡散係数に対する各種方程式モデルのモデリ
ング方法について具体的に説明するとともに，そのモデルを用いた単純な壁面乱流熱伝達場の予測値を示しながら，モデルの特徴を述べる．

４・１ 速度場モデル
４・１・１ ０方程式モデル
０方程式モデル（Zero equation model）はプラントルの混合距離理論（Prandtl’s Mixing Length
Theory （Prandtl, 1925））として知られているが，混合距離とは，図 3 に示すような平均速度分布を持つ乱流境界層（Turbulent boundary layer）において，微小距離ℓだけ離れた２つの流体塊を考えたとき，これら流体塊が自分の特性を有したまま動ける距離のことをいう．つまり，１つの流体塊がℓだけ動いて合体してしまうとその性質を失うため，その合体するまでの距離である．この時，これら流体塊の速度差𝛿𝑈を乱れの代表速度𝑉とする．ℓは十分に小さいので，𝛿𝑈はテイラー展
開を用い，高次項を無視することで以下のように近似できる．

- 15 -

JSME TED Newsletter, No.87, 2019

Fig. 3 Mixing length in turbulent boundary layer

̅
̅
̅
𝜕𝑈
1 𝜕2𝑈
𝜕𝑈
2
ℓ+
ℓ
+
⋯
−
𝑈(𝑦)
≈
ℓ
(30)
𝜕𝑦
2 𝜕𝑦 2𝜕𝑦
ここで，𝑦は壁からの垂直距離である．
̅ ⁄𝜕𝑦)と混合距離ℓを用いて，運動量の渦拡散係数を表式(30)で導かれた乱れの代表速度𝑉 = ℓ(𝜕𝑈す．
̅ (𝑦 + ℓ) − 𝑈
̅ (𝑦) = 𝑈̅ (𝑦) +
𝛿𝑈 = 𝑈𝜈𝑡 = 𝐶𝑉 ∙ 𝐿 = 𝐶 (ℓ

̅
̅
𝜕𝑈
𝜕𝑈
) ∙ ℓ = 𝐶ℓ2 ( )
𝜕𝑦
𝜕𝑦

(31)

ここで，モデル定数𝐶は１である．運動量の渦拡散係数は正値のため，速度勾配を絶対値として
以下のように書く．
𝜈𝑡 = ℓ2 |

̅
𝜕𝑈
|
𝜕𝑦

(32)

混合距離ℓは経験則に立脚してℓ = 𝜅𝑦として与えられ，𝜅はカルマン定数（von Kármán constant）と呼ばれる定数（𝜅 ≈ 0.4）であるため，渦拡散係数について新たな方程式を解く必要はない．これにより，レイノルズせん断応力は次式でモデル化される．
−𝑢𝑣 = 𝜈𝑡

̅
̅ 𝜕𝑈
̅
̅ 𝜕𝑈
̅
𝜕𝑈
𝜕𝑈
𝜕𝑈
= ℓ2 | |= (𝜅𝑦)2 | |
𝜕𝑦
𝜕𝑦 𝜕𝑦
𝜕𝑦 𝜕𝑦

(33)

この０方程式モデルは，乱流境界層に対して重要な関係を導く．その関係をここで導いてみる．
乱流境界層において，壁面せん断力𝜏は次式のように与えられる．
𝜏𝑤 = 𝜏𝑙 + 𝜏𝑡 = 𝜇

̅
𝜕𝑈
− 𝜌𝑢𝑣
𝜕𝑦

(34)

- 16 -

JSME TED Newsletter, No.87, 2019

̅ ⁄𝑈
𝑈

1

0.5

= 2442

DNS (Abe & Kawamura)
Zero-equation

0
0

0.5

𝑦⁄ 𝛿

1

(a) Normarized by channel width and bulk velocity

DNS (Abe & Kawamura)
Zero-equation

̅+
𝑈

20

̅ + = 1⁄𝜅𝑈

𝑦+ + 510= 2442

0 0
10

101

102

𝑦+

103

(b) Normarized by friction velocity
Fig. 4 Prediction of turbulent channel flow using zero-equation model

ここで𝜏𝑙 は粘性によるせん断力，𝜏𝑡 は乱流によるせん断力である．乱流境界層は，内層（Inner layer）
における壁面極近くに粘性底層（Viscous sublayer）が存在し，その上の遷移層（Buffer layer）に続き対数領域（Logarithmic region）が存在することが知られている．この対数領域ではせん断力は乱流によるせん断に支配されている．そのため，式(34)は𝜏𝑤 = −𝜌𝑢𝑣となる．これに０方程式モデルの式(33)を代入する．
𝜏𝑤 = −𝜌𝑢𝑣 = 𝜌(𝜅𝑦)2 |

̅ 𝜕𝑈
̅ 𝜏𝑤
̅ 𝜕𝑈
̅
𝜕𝑈
𝜕𝑈
|
⇒
= (𝜅𝑦)2 | |
𝜕𝑦 𝜕𝑦
𝜌
𝜕𝑦 𝜕𝑦

(35)

̅+ = 𝑈̅⁄𝑢𝜏
ここで，摩擦速度（Friction velocity）𝑢𝜏 = √𝜏𝑤 ⁄𝜌，無次元距離𝑦 + = 𝑦𝑢𝜏 ⁄𝜈，無次元速度𝑈を用いて式(35)を変形すると次式を得る．
̅ 𝜕𝑈
̅
̅
̅+
𝜏𝑤
𝜕𝑈
𝜕𝑈
𝜕𝑈
= 𝑢𝜏 2 = (𝜅𝑦)2 | |⇒ 𝑢𝜏 = 𝜅𝑦⇒ 1 = 𝜅𝑦 + +
𝜌
𝜕𝑦 𝜕𝑦
𝜕𝑦
𝜕𝑦
式(36)を積分すると次式を得る．

- 17 -(36)

JSME TED Newsletter, No.87, 2019

30
ℓ+

ℓ+
𝜅𝑦+

15

−𝑢𝑣 +
ℓ =
̅ +⁄𝜕𝑦+
𝜕𝑈
+

0
0

30

𝑦+

60

Fig. 5 Distribution of mixing length in boundary layer

̅+ = ∫
∫ 𝜕𝑈

1
1
̅+ =
𝜕𝑦 + ⇒ 𝑈
𝜅𝑦 +𝜅

𝑦+ + 𝐶(37)

ここで𝐶は積分定数である．この関係式は，対数速度分布（Logarithmic velocity profile）と呼ばれる乱流境界層を解析する際の重要な関係式であるが，定数とされている𝜅と𝐶が乱流境界層の状態によって変化するため，普遍的な関係式とならないのが，乱流モデルを使用するにあたっての難点となっていく．
さて，ここで，０方程式モデルによる完全発達した平行平板間乱流の予測結果を図 4 に示す．
計算において，境界条件は壁面で滑りなし条件を用いている．
図 4(a)はバルク速度で無次元化した結果を示し，図 4(b)は摩擦速度で無次元化した結果を示している．また，比較したデータは直接数値シミュレーション（Direct numerical simulation: DNS）の結果である（Abe and Kawamura, 2003）．図 4(b)には，式(37)の対数速度分布に対し，平行平板間
乱流の典型的な値である𝜅 = 0.4と𝐶 = 5とした分布を入れた．DNS による速度分布は，この対数
速度分布によく沿っていることが分かるが，０方程式モデルの予測値は全く沿っていない．左図で見ると DNS の結果とよく合っているように見えるが，摩擦速度で無次元化した結果では，なぜこれほどの違いが出るのであろうか？まず考えられるのが，予測された摩擦速度が大きく違っていることがある．それは，モデルの適用に何らかの問題があることが考えられ，モデルの使用が適切でないため，摩擦速度の予測に影響していることとなる．その問題を探るために，この混合̅ ⁄𝜕𝑦)として，DNS データ（Hattori et al., 2007）から逆算した結果を図 5 に示
距離をℓ = √−𝑢𝑣⁄(𝜕𝑈
す．明らかに，混合距離は壁面近傍で短くなり，ℓ = 𝜅𝑦の関係が成立していないことが分かる．
このことから，式(35)～(37)の関係から導出された対数速度分布は壁極近傍では成立しないことも分かる（図 5 からは，約𝑦 + > 50で混合距離が短くなっていることが分かるが，図 4(b)の DNS の速度分布からは約𝑦 + > 30で対数領域となっているように見える）．そのため，この不一致を補正するために，次の補正関数をつけた混合距離が提案された（Van Driest, 1956）．
ℓ = 𝜅𝑦 [1 − exp (−

𝑦+
)]
𝐴

(3 )

ここで𝐴はモデル定数で，𝐴 = 26である．この関数形は，以後，０方程式モデルのみならず，
多くの乱流モデルの壁面近傍挙動の補正に Van Driest 型補正関数（Van Driest’s function），もしく
は壁面補正関数（Wall dumping function）として多用されるようになる．そして逆にこのことから，
０方程式モデルから導出される対数速度分布は，対数領域で速度分布を与えることが出来るため，- 18 -

JSME TED Newsletter, No.87, 2019̅+𝑈

DNS (Abe & Kawamura)
Zero-equation
20
Zero-equation with WDF+
̅ = 1 ⁄𝜅𝑈
𝑦+ + 510= 2442

0 0
10

101

102 𝑦 +

103(a) with Van Driest’s function̅+𝑈

DNS (Abe & Kawamura)
Zero-equation with WDF
20
Zero-equation with OLM
Zero-equation+̅𝑈 = 1⁄ 𝜅𝑦+ + 510= 2442

0 0
10

101

102

𝑦+

103

(b) with outer layer modelFig. 6 Prediction of turbulent channel flow using zero-equation model with modifications

境界層の予測において壁面で厳密な境界条件を用いずに，対数速度分布を境界条件として予測計算を行う，いわゆる高レイノルズ数型乱流モデル（High-Reynolds-Number Turbulence Model）に応用されることとなる．
さて，この Van Driest 型補正関数を導入した０方程式モデルによる完全発達した平行平板間乱流の予測結果を図 6(a)に示す．明らかに補正関数を導入したモデルが，壁面から対数領域にかけて，補正関数を用いない場合より良好な予測値を与えていることが分かる．なお，𝑦 +が大きい場
所で予測値が合わないのは，そこは対数領域を超えた外層領域（Outer layer，もしくは後流領域；
Wake region）であり，ℓ = 𝜅𝑦の関係が成立しないからであり，さらなる補正が必要となる．そのため，外層に対して，以下のモデル（Cebeci, 1973）を入れて計算した結果を図 6(b)に示す．
∞

̅𝑒 − 𝑈̅ )𝑑𝑦| ; 𝑦𝑖 ≤ 𝑦 ≤ 𝛿𝜈𝑡 = α |∫ (𝑈

(39)

0

ここで，α = 0.0168 であり（Cebeci, 1973），𝑦𝑖 は境界層厚さ𝛿99 の 30%とした．図 6(b)より，外層領域における速度増加を適切に予測していることが分かる．なお，外層の混合距離に対するモデ- 19 -

JSME TED Newsletter, No.87, 2019
リングとして，実験データに基づいた以下の式も提案されている（Kays and Crawford, 1993）．
ℓ = λ𝛿99 ; 𝑦 ≥ 0.7𝛿99(40)

ここで，𝜆 = 0.0 5である．このモデルを用いて行った計算においても，図 6(b)と同様の結果が得られる．
以上示したように，０方程式モデルは，補正関数を導入することで壁面から対数領域までは，ℓ = 𝜅𝑦の関係式が使用できるが，この関係が成立しない外層領域では，別のモデリングが必要となる．カルマン定数𝜅も，様々な乱流場では変化し，また外層のモデルに対してモデル定数が普遍的ではないことや，壁面噴流など境界層厚さが規定し難い乱流場を考慮すると，実用的な乱流予測に用いるには難しい．
４・１・２ １方程式モデル
４・１・２（a） 𝑘 − ℓ１方程式モデル
０方程式モデルにおける運動量の渦拡散係数のモデリングは，乱流境界層に対して主に対数領域でのみ成立し，他の領域については別途補正等が必要であった．これを解消するためには，運動量の渦拡散係数を構成するスケールの１つ以上は輸送方程式を解くことが望ましい．そこで，乱れの長さスケールは混合距離を用い，乱れの速度スケールについては，乱流エネルギー𝑘（=𝑢𝑖 𝑢𝑖 ⁄2）̅̅̅̅̅
（Turbulent kinematic energy）から定まる√𝑘をそのスケールとし，次式のように運動量の渦拡散係数を与える（Kolmogorov, 1942, Prandtl, 1945）．
𝜈𝑡 = 𝐶ℓ √𝑘 ∙ ℓ(41)

ここで，𝐶ℓ はモデル定数であり，その値については後述する．そして，このモデルでは，速度スケール√𝑘の輸送方程式を解くことから，𝑘 − ℓ１方程式モデルと呼ぶ．
さて，乱流エネルギーの輸送方程式は，本企画第１回の２・１節で述べたように，Navier-Stokes
方程式より導出することが出来る．まず，式(9)で与えらえた Navier-Stokes 方程式の速度，圧力
̅𝑖 + 𝑢𝑖 ，𝑝 = 𝑃̅ + 𝑝のようにレイノルズ分解を施すと次式を得る．
に，𝑢𝑖 = 𝑈̅𝑖 𝑈̅𝑗 𝜕𝑢𝑖 𝑢𝑗 𝜕𝑈̅𝑖 𝑢𝑗 𝜕𝑢𝑖 𝑈
̅𝑗 𝜕𝑢𝑖 𝑢𝑗 1 𝜕𝑃̅ 1 𝜕𝑝
̅𝑖 𝜕𝑢𝑖 𝜕𝑈
̅𝑖𝜕𝑈
𝜕2𝑈
𝜕 2 𝑢𝑖
+
+
+
+
+
+
=
+
+𝜈
+𝜈
𝜕𝑡
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜌 𝜕𝑥𝑗 𝜌 𝜕𝑥𝑗
𝜕𝑥𝑘 𝜕𝑥𝑘𝜕𝑥𝑘 𝜕𝑥𝑘(42)

ここで，式(9)の左辺第２項は，連続の式（𝜕𝑢𝑗 ⁄𝜕𝑥𝑗 = 0）を用いて，𝑢𝑗 𝜕𝑢𝑖 ⁄𝜕𝑥𝑗 = 𝜕𝑢𝑖 𝑢𝑗 ⁄𝜕𝑥𝑗 と保存系に変形してからレイノルズ分解を施した．さらに式(40)の時間平均をとると式(13)と同様の次式が得られる．
̅𝑖
̅𝑖
𝜕𝑈
𝜕
1 𝜕𝑃̅
𝜕2𝑈
̅𝑖 𝑈
̅𝑗 + ̅̅̅̅̅
(43)
+
(𝑈
𝑢𝑖 𝑢𝑗 ) = −+𝜈𝜕𝑡 𝜕𝑥𝑗𝜌 𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗
次に式(40)から式(41)を引くと，変動速度成分𝑢𝑖 に関する輸送方程式を得る．
𝜕𝑢𝑖
𝜕
1 𝜕𝑝
𝜕 2 𝑢𝑖̅𝑖 𝑢𝑗 + 𝑢𝑖 𝑈
̅𝑗 + 𝑢𝑖 𝑢𝑗 − ̅̅̅̅̅)+(𝑈𝑢𝑖 𝑢𝑗 =
+𝜈𝜕𝑡 𝜕𝑥𝑗𝜌 𝜕𝑥𝑗
𝜕𝑥𝑘 𝜕𝑥𝑘(44)

̅𝑗 ⁄𝜕𝑥𝑗 = 0）を用いて，左辺第２項を変形し，次式を得る．
ここで，連続の式（𝜕𝑢𝑗 ⁄𝜕𝑥𝑗 = 0，𝜕𝑈
̅𝑖
𝜕𝑢𝑖
𝜕𝑈
𝜕𝑢
𝜕
1 𝜕𝑝
𝜕 2 𝑢𝑖
̅𝑗 𝑖 ++ 𝑢𝑗+𝑈(𝑢𝑖 𝑢𝑗 − ̅̅̅̅̅)𝑢𝑖 𝑢𝑗 =
+𝜈𝜕𝑡𝜕𝑥𝑗
𝜕𝑥𝑗 𝜕𝑥𝑗𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘
さらに，式(43)における添字𝑗を𝑘に置換すると，式(18)を得る．

- 20 -(45)

JSME TED Newsletter, No.87, 2019
̅𝑖
𝜕𝑢𝑖
𝜕𝑈
𝜕𝑢
𝜕
1 𝜕𝑝
𝜕 2 𝑢𝑖
̅𝑘 𝑖 +(𝑢𝑖 𝑢𝑘 − 𝑢𝑖 𝑢𝑘 ) =
+ 𝑢𝑘
+𝑈
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘(46)

ここで，変動速度成分𝑢𝑗 に関する輸送方程式を得るために，添字𝑖を𝑗に置換すると，式(19)を得
る．
̅𝑗
𝜕𝑢𝑗
𝜕𝑈
𝜕𝑢
𝜕 2 𝑢𝑗𝜕
1 𝜕𝑝
̅𝑘 𝑗 ++ 𝑢𝑘+𝑈(𝑢𝑗 𝑢𝑘 − 𝑢𝑗 𝑢𝑘 ) =
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘(47)

式(46) × 𝑢𝑗 + 式(47) × 𝑢𝑖 とすることで，第 1 回の式(20)で示したレイノルズ応力𝑢𝑖 𝑢𝑗 の輸送方
程式を得る．
̅̅𝜕𝑢𝑖 𝑢𝑗𝜕𝑢 𝑢
𝜕𝑈𝜕𝑢 𝑢 𝑢𝜕𝑈𝑝 𝜕𝑢 𝜕𝑢̅𝑘 𝑖 𝑗 = − (𝑢𝑖 𝑢𝑘 𝑗 + 𝑢𝑗 𝑢𝑘 𝑖 ) − 𝑖 𝑗 𝑘 + ( 𝑖 + 𝑗 )
+𝑈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜌 𝜕𝑥𝑗 𝜕𝑥𝑖
𝜕 2 𝑢𝑖 𝑢𝑗
1 𝜕𝑝𝑢𝑖 𝜕𝑝𝑢𝑗𝜕𝑢𝑖 𝜕𝑢𝑗
− (
+
)+𝜈
− 2𝜈
𝜌 𝜕𝑥𝑗
𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘(4 )乱流エネルギー𝑘に対する輸送方程式は，式(46)の添字𝑗を𝑖と置けば得られる．
̅̅
𝜕𝑢𝑖 𝑢𝑖𝜕𝑢 𝑢
𝜕𝑈𝜕𝑈𝜕𝑢 𝑢 𝑢𝑝 𝜕𝑢 𝜕𝑢̅𝑘 𝑖 𝑖 = − (𝑢𝑖 𝑢𝑘 𝑖 + 𝑢𝑖 𝑢𝑘 𝑖 ) − 𝑖 𝑖 𝑘 + ( 𝑖 + 𝑖 )
+𝑈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜌 𝜕𝑥𝑖 𝜕𝑥𝑖
1 𝜕𝑝𝑢𝑖 𝜕𝑝𝑢𝑖
𝜕 2 𝑢𝑖 𝑢𝑖
𝜕𝑢𝑖 𝜕𝑢𝑖
− (
+
)+𝜈
− 2𝜈
𝜌 𝜕𝑥𝑖
𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘(49)

ここで，𝑘 = ̅̅̅̅̅
𝑢𝑖 𝑢𝑖 ⁄2であるので，式(47)は次式のようになる．
̅𝑖 𝜕𝑢𝑗 𝑘 ′ 1 𝜕𝑢𝑗 𝑝
𝜕𝑘
𝜕𝑘
𝜕𝑈
𝜕2𝑘
𝜕𝑢𝑖 𝜕𝑢𝑖
̅𝑗+𝑈= −𝑢𝑖 𝑢𝑗
−
−
+𝜈
−𝜈
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜌 𝜕𝑥𝑗
𝜕𝑥𝑗 𝜕𝑥𝑗𝜕𝑥𝑗 𝜕𝑥𝑗(50)

ここで，𝑘 ′ = 𝑢𝑖 𝑢𝑖 ⁄2，微分の添字は𝑗と置換した．この輸送方程式において，両辺の各項はそれぞれ，左辺第１項から時間項，対流項（Convection term），右辺第１項から生成項（Productionterm），乱流拡散項（Turbulent diffusion term），圧力拡散項（Pressure diffusion term），粘性
拡散項（molecular diffusion term），散逸項（Dissipation term）である．この輸送方程式を解け
ば乱流エネルギー𝑘の分布が得られることになるが，右辺第２項の乱流拡散項，同３項の圧力拡散項，同５項の散逸項は，このままでは解くことが出来ないため，モデル化が必要である．
乱流拡散項と圧力拡散項は，まとめて勾配拡散型モデルを用いて次式のようにモデル化される．
𝜕𝑘′̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
− (𝑢𝑗 𝑘 ′ + 𝑢𝑗 𝑝⁄𝜌) = −𝑢𝑗 (𝑘 + 𝑝⁄𝜌) = (𝜈𝑡 )𝑘
𝜕𝑥𝑗

(51)

ここで，(𝜈𝑡 )𝑘 は乱流エネルギーの渦拡散係数である．(𝜈𝑡 )𝑘 も場の変数であるため，まだモデル化の必要がある．ここでは，プラントル数（Prandtl number）（𝑃𝑟 = 𝜈⁄𝛼）のように，場の拡散係数の比で表す量を，乱流量に拡大定義した乱流プラントル数（Turbulent Prandtl number）の概念を用いて，以下のようにモデル化する（勾配拡散型モデルの一般化，乱流プラントル数の詳細については後述する）．
𝜎𝑘 =𝜈𝑡(𝜈𝑡 )𝑘(52)

上式は，乱流エネルギーの拡散に対する乱流プラントル数𝜎𝑘 を，運動量と乱流エネルギーの渦- 21 -

JSME TED Newsletter, No.87, 2019

DNS (Abe & Kawamura)
One-equation
Zero-equation

̅+
𝑈

20

̅ + = 1⁄𝜅𝑈

𝑦+ + 510= 2442

0 0
10

101

102 𝑦 +

103(a) without outer layer model for mixing length̅+𝑈

DNS (Abe & Kawamura)
One-equation20
One-equation with OLM
̅ + = 1⁄𝜅𝑈
𝑦+ + 510= 2442

0 0
10

101

102

𝑦+

103

(b) with outer layer model for mixing length
Fig. 7 Prediction of turbulent channel flow using one-equation model

拡散係数の比で定まる量と定義してモデル化している．厳密には，この𝜎𝑘 も場の関数であり，このモデル化によって(𝜈𝑡 )𝑘 が与えられることにはならないが，通常𝜎𝑘 を定数として扱い，(𝜈𝑡 )𝑘 を与える．そのため式(46)は，次式となる．
′
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅−𝑢𝑗 (𝑘 + 𝑝⁄𝜌) =𝜈𝑡 𝜕𝑘𝜎𝑘 𝜕𝑥𝑗(53)

ここで，𝜎𝑘 は典型的な値は𝜎𝑘 = 1であるが，1⁄𝜎𝑘 = 0.3 も使用される．次に，散逸項のモデル化であるが，次元解析的に混合距離を用いて次式で与える．
𝜈𝜕𝑢𝑖 𝜕𝑢𝑖𝑘 3⁄2= 𝐶𝐷𝜕𝑥𝑗 𝜕𝑥𝑗ℓ(54)

ここで，𝐶𝐷 はモデル定数である．以上により，乱流エネルギーのモデル化された輸送方程式は- 22 -

JSME TED Newsletter, No.87, 2019
次式のようになる．
̅𝑖
𝜕𝑘
𝜕𝑘
𝜕
𝜈𝑡 𝜕𝑘
𝜕𝑈
𝑘 3⁄2
̅𝑗
+𝑈
=
[(𝜈 + )
] − 𝑢𝑖 𝑢𝑗− 𝐶𝐷𝜕𝑡𝜕𝑥𝑗 𝜕𝑥𝑗
𝜎𝑘 𝜕𝑥𝑗𝜕𝑥𝑗ℓ(55)

モデル定数𝐶ℓと𝐶𝐷 については，理論的な考察により決めることが出来る．まず，対数速度分布を求めた式(35)より，次式を得る．
−𝑢𝑣 ≅

̅ 𝜕𝑈
̅
̅ 𝑢𝜏
𝜏𝑤
𝜕𝑈
𝜕𝑈
= 𝑢𝜏2 = ℓ2 | |
⇒
=
𝜌
𝜕𝑦 𝜕𝑦
𝜕𝑦
ℓ

(56)

ここで，新たに−𝑢𝑣⁄𝑘 = √𝐶μ という乱流における構造パラメーター（Structural parameter）を導入し，式(35)に代入すると次式を得る．
−𝑢𝑣 = 𝑢𝜏2 ⇒

𝑢𝜏2
𝑢𝜏 4
= √𝐶μ ⇒ 𝐶μ = ( )
𝑘
√𝑘

(57)

この時，構造パラメーター𝐶μ と，モデル定数𝐶ℓ ，𝐶𝐷 の関係は𝐶μ = 𝐶𝐷 𝐶ℓ である（この関係式は，後述する𝑘 − 𝜀２方程式モデルにおける渦拡散係数のモデリングと整合性を取っているためである）．また，乱流エネルギーの輸送方程式は，対数領域で局所平衡（local equilibrium）が成立すると考え，そこでは生成＝散逸であるため，この関係に式(56)と式(57)を代入すると，次式を得る．
−𝑢𝑣

̅
𝜕𝑈
𝑘 3⁄2
𝑢𝜏
𝑘 3⁄2
𝑢𝜏 3
3⁄4
= 𝐶𝐷
⇒ 𝑢𝜏2 ∙ = 𝐶𝐷
⇒ 𝐶𝐷 = ( ) = 𝐶μ
𝜕𝑦
ℓ
ℓ
ℓ
√𝑘

(5 )

実験や DNS データより，この構造パラメーター√𝐶μ は，約 0.3（𝐶ℓ = 0.09）となることが知られている．そのため，この値を用いると𝐶𝐷 =0.164 を得る．また，𝐶μ = 𝐶𝐷 𝐶ℓの関係より，𝐶ℓ = 0.55となるが，研究の進展により𝐶ℓ = 0.56，𝐶𝐷 =0.18 が使われる．
図 7 に 1 方程式モデルで予測した完全発達した平行平板間乱流の平均速度分布を示す．計算においては，壁面から解くため，混合距離は壁面補正された式(38)を用い，モデル定数𝐶𝐷 は𝐶𝐷 = 1とした．比較のため，DNS の結果と壁面補正を入れた方程式モデルの予測結果を図に入れてある．
1 方程式モデルの予測結果を見ると，外層における後流領域の予測が不十分で，０方程式モデルの結果と全く変わらない．この原因としては，長さスケールに混合距離を用いる以上，この外層領域の予測は難しく，そのため，混合距離を用いた散逸項のモデル化にも影響が出たと考えられる．そのため，式(40)で与えられた外層における混合距離モデルを入れて行った計算結果を図 7(b)に示す．混合距離のモデルを適切にすることによって，外層における速度分布を適切に予測できることが分かる．１方程式モデルにおいては，混合距離のモデルを１つにして，その不足分を補えるだけの速度スケールが算出できることを期待したが，単純な乱流場においても，速度スケー
ルとした乱流エネルギーの輸送方程式を解くだけでは予測性能の向上が見出せず，今日においても𝑘 − ℓ１方程式モデルは，RANS であまり利用例を見ない．
４・１・２（b） 渦拡散係数輸送方程式モデル
１方程式モデルには，前述した運動量の渦拡散係数を，乱れの速度スケールと乱れの長さスケールでモデル化する試みの他，直接，渦拡散係数の輸送方程式で解くというモデルも存在する（Neeand Kovasznay, 1969, Baldwin and Barth, 1990, Spalart and Allmaras, 1994, Menter, 1997, Nagano, Pei
and Hattori, 2000）．熱流体の支配方程式の他，渦拡散係数の輸送方程式を１つ解くだけなので，これも１方程式モデルで分類される．そのモデル化の方法は，後述の𝑘 − 𝜀２方程式モデルを基本としたモデルと，様々な考察を経て構築されたモデルが最近の主流である．前者は，BB モデル（Baldwin and Barth, 1990），Menter モデル（Menter, 1997），NPH モデル（Nagano, Pei and Hattori,
2000），後者は SA モデル（Spalart and Allmaras, 1994）がある．ここで，乱流モデル名は，その開発者名の頭文字を繋げて呼ばれることが多いため，ここでもそのような表記とした．

- 23 -

JSME TED Newsletter, No.87, 2019
まず，𝑘 − 𝜀２方程式モデルを基本としたモデルについて記す．このモデルは，その導出過程が明確なため，比較的理解しやすいモデルである．
運動量の渦拡散係数の輸送方程式の導出については，その渦拡散係数を，乱れの速度スケールには𝑘 − ℓ１方程式モデルと同じく乱流エネルギーから定まる√𝑘を，乱れの長さスケールには，乱流エネルギーの散逸率εを用いて定義した乱れの長さスケール𝑘 3⁄2⁄εを用いて，𝜈̃𝑡 = 𝐶μ √𝑘 ∙(𝑘 3⁄2⁄ε) = 𝐶μ 𝑘 2⁄εと置くことから始まる（この乱れの長さスケールの決定方法については，２方程式モデルの節で述べる）．ただし，運動方程式における渦粘性係数は，𝜈𝑡 = 𝐷2 𝜈̃𝑡 で与える．この関係から，渦拡散係数の実質微分を取ると，次式のようになる．
𝐷𝜈̃𝑡
𝐷
𝑘2
2𝑘 𝐷𝑘 𝑘 2 𝐷ε
2𝜈̃𝑡 𝐷𝑘 𝜈̃𝑡 𝐷ε=(𝐶μ ) = 𝐶μ (− 2 ) = 𝐶μ (
−
)
𝐷𝑡
𝐷𝑡
ε
ε 𝐷𝑡 ε 𝐷𝑡
𝑘 𝐷𝑡 ε 𝐷𝑡(59)

̅𝑗 𝜕⁄𝜕𝑥𝑗 である．右辺の乱流エネルギーの実質微分𝐷𝑘⁄𝐷𝑡（𝐷𝑘⁄𝐷𝑡 =式
ここで，𝐷⁄𝐷𝑡 = 𝜕⁄𝜕𝑡 + 𝑈(50)右辺とする．ここで散逸項𝐶𝐷 (𝑘 3⁄2 ⁄ℓ)は，𝐶𝐷 (𝑘 3⁄2 ⁄ℓ) = εと置き換えるが，導出に用いる混合̅⁄𝜕𝑦⁄(𝜕 2𝑈̅⁄𝜕𝑦 2 )|である
距離は von Kármán が提案した混合距離ℓ𝐸 = κ|𝜕𝑈（von Kármán, 1930））と，⁄
乱流エネルギー散逸率の実質微分𝐷ε 𝐷𝑡（散逸率ε の輸送方程式の導出は２方程式モデルの節で述べる）にそれぞれの輸送方程式を代入し，さらに𝑘とεを，𝜈̃𝑡 を用いて表現すると，𝜈̃𝑡 の輸送方程式を得る（Menter, 1997）．
𝐷𝜈̃𝑡𝜕𝜈̃𝑡 𝜕𝜈̃𝑡𝐸𝑘𝜀=[(𝜈 + )
] + 𝐶1 𝐷1𝜈̃𝑡 Ω − 𝐶2 𝐶3𝐸𝐵𝐵 ta h ()
𝐷𝑡𝜕𝑥𝑗𝜎𝜈 𝜕𝑥𝑗
𝐶3 𝐸𝐵𝐵(60)

̅𝑖 ⁄𝜕𝑥𝑗 − 𝜕𝑈̅𝑗 ⁄𝜕𝑥𝑖 )⁄2
ここで，𝜎𝜈 ，𝐶1，𝐶2，𝐶3はモデル定数であり，Ω = √Ω𝑖𝑗 Ω𝑖𝑗 は渦度，Ω𝑖𝑗 = (𝜕𝑈
は渦度テンソルである．また，モデル関数𝐷1は次式で与えられる．
𝐷1 =𝜈𝑡 + 𝜈
,𝜈̃𝑡 + 𝜈𝜈𝑡 = 𝐷2 𝜈̃𝑡 ,𝐷2 = 1 − exp (−

𝜈̃𝑡 2
)
𝐴+ 𝜅𝜈

(61)

さらに，右辺最終項における𝐸𝐵𝐵 ，𝐸𝑘𝜀 は次式で与えられる（Menter, 1997）
．

𝐸𝐵𝐵 =

𝜕𝜈̃𝑡 𝜕𝜈̃𝑡
,𝜕𝑥𝑗 𝜕𝑥𝑗̅𝑖
𝜕2 𝑈

𝐸𝑘𝜀 = 𝜈̃𝑡 2 (̅𝑖
𝜕2𝑈𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑘 𝜕𝑥𝑘
)
̅𝑠 𝜕𝑈
̅𝑠
𝜕𝑈

(62)

𝜕𝑥𝑙 𝜕𝑥𝑙

次に，上述のモデルのように，𝑘 − 𝜀２方程式モデルの輸送方程式から数学的に導出することな
く，様々な考察を経て構築された渦拡散係数輸送方程式モデルについて述べる．このモデルは SA
モデルと呼ばれ，遷移問題も扱うことが出来るため，上述のように航空機の翼周りなどの計算によく用いられるようである．そのモデルの輸送方程式は次式で与えられる（Spalart and Allmaras,
1994）
．
2

𝐷𝜈̃𝑡
1 𝜕
𝜕𝜈̃𝑡
𝐶𝑏2 𝜕𝜈̃𝑡𝜈̃𝑡 2=[(𝜈 + 𝜈̃𝑡 )] + 𝐶𝑏1 𝜈̃𝑡 𝑆̃ −( ) − 𝐶𝑤1𝑓𝑤 ( )𝐷𝑡𝜎𝜈 𝜕𝑥𝑗𝜕𝑥𝑗
𝜎𝜈 𝜕𝑥𝑗𝑑(63)

ここで，運動方程式における渦粘性係数は𝜈𝑡 = 𝑓𝑣1 𝜈̃𝑡 で与えられ，𝐶𝑏1，𝐶𝑏2，𝜎𝜈 はモデル定数であ3 )る ． モ デ ル 関 数 は そ れ ぞ れ 𝑓𝑣1 = 𝜒 3⁄(𝜒 3 + 𝐶𝑣1， 𝜒 = κ𝑦 + ， 𝑆̃ = 𝑆 + (𝜈̃𝑡 ⁄𝜅 2 𝑑2 )𝑓𝑣2 ， 𝑓𝑣2 = 1 −𝜒⁄(1 + 𝜒𝑓𝑣1 )であり，関数中のκはカルマン定数，𝑑は壁への距離，𝑦 + は壁面垂直方向無次元距離
̅𝑖 ⁄𝜕𝑥𝑗 + 𝜕𝑈̅𝑗 ⁄𝜕𝑥𝑖 )⁄2は歪テンソルである．
である．また，𝑆 = √2𝑆𝑖𝑗 𝑆𝑖𝑗 は歪パラメーター，𝑆𝑖𝑗 = (𝜕𝑈
モデルの輸送方程式(63)と，前述の𝑘 − 𝜀２方程式モデルをベースとした渦拡散係数輸送方程式(60)を比較すると，右辺第４項の消滅項に違いがみられるのみであり，結果としてほぼ同様の形となっているのは，２方程式モデルを含めた長いモデル研究の蓄積からである．

- 24 -

JSME TED Newsletter, No.87, 2019

DNS (Abe & Kawamura)
One-equation with OLM

̅+
𝑈

20

̅ + = 1⁄ 𝜅𝑈𝑦+ + 5Menter model
SA model10= 24420 0
10102 𝑦 +

101

103

Fig. 8 Prediction of turbulent channel flow using eddy viscosity transport models

さて，この SA モデルで，翼面上などの層流から乱流へ遷移する場合の計算を行うには，式(63)を次式のようにして扱う（Spalart and Allmaras, 1994）
．
2

𝐷𝜈̃𝑡
1 𝜕
𝜕𝜈̃𝑡
𝐶𝑏2 𝜕𝜈̃𝑡
𝐶𝑏1
𝜈̃𝑡 2
̃
=
[(𝜈 + 𝜈̃𝑡 )] + 𝐶𝑏1 (1 − 𝑓𝑡2 )𝜈̃𝑡 𝑆 −( ) − (𝐶𝑤1 𝑓𝑤 − 2 𝑓𝑡2 ) ( ) + 𝑓𝑡1 Δ𝑈 2
𝐷𝑡𝜎𝜈 𝜕𝑥𝑗𝜕𝑥𝑗
𝜎𝜈 𝜕𝑥𝑗𝜅𝑑(64)

ここで，Δ𝑈は遷移を考える点と壁面との速度差のノルムであり，モデル関数は𝑓𝑡1 ，𝑓𝑡2 次式で与えられる．
𝑓𝑡1 = 𝐶𝑡1 𝑔𝑡 exp [−𝐶𝑡2𝜔𝑡2(𝑑 2 + 𝑔𝑡2 𝑑𝑡2 )] ,Δ𝑈 2

𝑓𝑡2 = 𝐶𝑡3 exp(−𝐶𝑡4𝜒 2 )(65)

ここで，𝐶𝑡1，𝐶𝑡2，𝐶𝑡3は，モデル定数，モデル関数𝑔𝑡 は，𝑔𝑡 = mi (0.1, Δ𝑈⁄𝜔𝑡 Δ𝑥𝑡 )で与えられ，
Δ𝑥𝑡 は遷移点における壁面上の格子幅，𝑑𝑡 は空間から遷移点への距離，𝜔𝑡 は遷移点における壁面上の渦度である．
図 8 にこれらの運動量の渦拡散係数輸送方程式モデルによって予測した発達した平行平板間乱流の平均速度分布を示す．Menter モデルが遷移域で若干の過少予測を示しているのは，壁面近傍におけるレイノルズ応力を過剰予測しているため，対数領域への遷移が早く予測されたためである．また，SA モデルは外層における速度増加を若干過少予測し，Menter モデルは最大速度である
中心速度はより良く予測しているように見えるが，外層モデルを導入した 1 方程式モデルのような，対数速度分布より上昇した明確な後流領域は分かり難い結果を与えていることが分かる．しかしながら，これら渦拡散係数輸送方程式モデルは，熱流体の支配方程式の他，1 つの輸送方程式を解くだけでほぼ良好な計算結果が得られることより，高レイノルズ数流れなど，計算負荷をより低く抑えたい熱流体の解析に用いられることが多い．
４・１・３ ２方程式モデル
現在の乱流熱流体解析で，その予測精度と計算負荷のバランスがとれているため，最も使用されているのが２方程式モデルである．１方程式モデルでは，運動量の渦拡散係数を表現する乱れの速度スケール√𝑘はその輸送方程式から解き，乱れの長さスケールを実験的な経験則である混合距離理論で定めていたが，乱れの長さスケールについて，乱流の特性周波数𝑓 = √𝑘⁄𝐿を微分方程式で解き長さスケールを与えたモデル（Kolmogorov, 1942），もしくは乱れの長さスケールそのものを輸送方程式で解き与えたモデル（Rotta, 1951, Hallow and Nakayania, 1967）が，運動量の渦拡- 25 -

JSME TED Newsletter, No.87, 2019
散係数を与えるために２つ方程式を解いたため，２方程式モデル（Two-equation model）の始まりとなっている．その後，長さスケールの概念に，エネルギー保有渦（Energy containing eddy）の持
つ乱流エネルギーの減衰率𝜀とで定まる長さスケール（𝐿 = 𝑘 3⁄2⁄𝜀 ），言い換えると，エネルギー
保有渦が持つ乱流エネルギーが消滅するまでの時間スケール（τ = 𝐿⁄𝑉 = 𝑘⁄𝜀）の概念を導入した
𝑘 − 𝜀２方程式モデルが提案され（Launder and Spalding, 1974），現在の主流モデルとなっている．
その他，乱れの時間スケールを乱流エネルギーと長さスケールと組み合わせた𝑘𝐿で表現し，その輸送方程式を解いた𝑘 − 𝑘𝐿２方程式モデル（Rodi and Spalding, 1970），渦度を用いた𝑘 − 𝑊２方程式モデル
（Spalding, 1982）が開発された．
さらに，上述した Kolmogrov のモデル（Kolmogorov, 1942）
を起源として，後に Saffman（1970）によって提唱され，その後応用工学問題用に改良された𝑘 − ω
２方程式モデル（Wilcox and Alber, 1972, Wilcox, 1988, Wilcox, 1998）が最近はよく使われている，
このように様々なバージョンが存在する２方程式モデルは，多くの研究者により改良が行われ，ますます乱流熱流体解析の主流モデルとなった．さらに，本企画第１回の３・１節で述べたように，渦拡散係数の概念を導入した勾配拡散型モデルの原理的な限界もあるが，その限界を克服す
るモデリング研究も進んでいる．ここでは，𝑘 − 𝜀２方程式モデルを中心に，モデルの改良例も示しながら，説明する．

４・１・３（a） 𝑘 − 𝜀２方程式モデル
𝑘 − 𝜀 ２方程式モデルでは，運動量の渦拡散係数を乱れの速度スケール𝑉を乱流エネルギーから定まる√𝑘と，エネルギー保有渦（Energy containing eddy）における乱流エネルギーの減衰率（散逸率）𝜀（Dissipation rate of turbulent kinematic energy）から定まる乱れの長さスケール𝐿=𝑘 3⁄2⁄𝜀を用い，次式のように表現する．
𝜈𝑡 = 𝐶𝜇 √𝑘 ∙

𝑘 3⁄2
𝑘2
= 𝐶𝜇
𝜀
𝜀

(66)

ここで，𝐶𝜇 はモデル定数である．この渦粘性係数を表現する変数である，乱流エネルギー𝑘と，その散逸率𝜀の２つの輸送方程式を解くことから２方程式モデルと呼ばれ，特に𝑘 − 𝜀 ２方程式モデル，𝑘 − 𝜀モデルと呼ばれる．
このエネルギー保有渦における乱流エネルギーの減衰率（散逸率）𝜀は，次式で定義される．
𝑑𝑘
= −𝜀
𝑑𝑡

(67)

エネルギー保有渦𝑡 = 𝑡0𝑡=𝑡

粘性により消散（熱に変わる）流れ

𝑘=0
𝑘 = 𝑘0𝑡
Fig. 9 Cascade process of energy containing eddies

これは，図 9 に示されるように，格子によって生成され，初期状態（𝑡 = 𝑡0）で，ある乱流エネ
ルギー（𝑘 = 𝑘0）を持つ乱流渦を考え，その乱流渦が流れに乗りながら時間の経過とともに減衰
し，時間𝑡 = 𝑡で，その乱流エネルギーが粘性によって消散し，熱に変わるまでの時間変化率を𝜀で- 26 -

JSME TED Newsletter, No.87, 2019
表現した式である．これより，乱流渦が減衰し，消散するまでの時間は，式(67)を積分することによって求めることが出来る．
𝑡𝑡𝑑𝑘𝑑𝜏 = − ∫ 𝜀 𝑑𝜏 ⇒ 𝑘(𝑡) − 𝑘(𝑡0 ) = −[𝜀(𝑡) − 𝜀(𝑡0 )]𝜏
𝑡0 𝑑𝑡
𝑡0

∫

(6 )

これより，𝑘(𝑡) = 0，𝜀(𝑡0 ) = 0なので，𝜏 = 𝑘⁄𝜀を得る．この𝜏が，エネルギー保有渦が減衰して熱エネルギーに変わるまでの時間である．よって，長さスケールは，時間スケールと速度スケール
から，𝐿 = 𝑉 ∙ 𝜏 = √𝑘 ∙ (𝑘 2 ⁄𝜀) = 𝑘 3⁄2⁄𝜀と求めることが出来るが，この概念から，わざわざ長さス
ケールを用いずに，時間スケールを用いて，渦粘性係数を𝜈𝑡 ~ 𝑉 ∙ 𝑉𝜏 ~ 𝑉 2𝜏と，乱れのエネルギー×乱れの時間スケールで定義すればよいことになる．この概念は，後述の温度場モデルで有用となる．
さて，乱流エネルギーの散逸率であるが，具体的には次式で定義される．

𝜀=𝜈

𝜕𝑢𝑖 𝜕𝑢𝑖 𝜕𝑢𝑗𝜕𝑢𝑖 𝜕𝑢𝑖(+)=𝜈𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗(69)

これは，式(50)で示した乱流エネルギーの輸送方程式，本企画第１回の２・１節の式(3)で示したエネルギー方程式に現れる．そのため，エネルギー保有渦が消散し，熱へと変化した乱流エネ
ルギーは，エネルギー方程式に反映されていることとなるが，非圧縮性を仮定した流体では無視できることが多い．
𝑘 − 𝜀 ２方程式モデルでは，この散逸率を輸送方程式で解くこととなるため，その輸送方程式を導出する．導出方法は，式(42)～(50)で示した乱流エネルギー（レイノルズ応力）の輸送方程式の導出と同様である．まず，レイノルズ応力の輸送方程式(49)の右辺最終項にある非等方散逸率̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜀𝑖𝑗 = 𝜈(𝜕𝑢𝑖 ⁄𝜕𝑥ℓ )(𝜕𝑢𝑗 ⁄𝜕𝑥ℓ )の式を作るため，変動速度𝑢𝑖 の式(46)と，変動速度𝑢𝑗 の式(47)に対し，
以下の操作を施す（式(46)と(47)は再掲）．
𝜈̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑗 𝜕
𝜕𝑢𝑖 𝜕
(equatio (46)) + 𝜈(equatio (47))
𝜕𝑥ℓ 𝜕𝑥ℓ
𝜕𝑥ℓ 𝜕𝑥ℓ

̅𝑖
𝜕𝑢𝑖
𝜕𝑈
𝜕𝑢
𝜕
1 𝜕𝑝
𝜕 2 𝑢𝑖
̅𝑘 𝑖 +(𝑢𝑖 𝑢𝑘 − 𝑢𝑖 𝑢𝑘 ) =
+ 𝑢𝑘
+𝑈
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘
̅
𝜕𝑢𝑗
𝜕𝑈𝑗
𝜕𝑢
𝜕 2 𝑢𝑗𝜕
1 𝜕𝑝
̅𝑘 𝑗 ++ 𝑢𝑘+𝑈(𝑢𝑗 𝑢𝑘 − 𝑢𝑗 𝑢𝑘 ) =
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘

(70)

(46)
(47)

ここで，式(70)の第１項は次式となる．
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅𝑖𝜕𝑢𝑗 𝜕 2 𝑢𝑖
𝜕𝑢𝑗 𝜕
𝜕𝑢𝑗 𝜕
𝜕𝑢
𝜕𝑈
𝜕𝑢
𝜕2
̅𝑘 𝑖 ) + 𝜈 𝑗(𝑢 𝑢 − 𝑢𝑖 𝑢𝑘 )
+𝜈
(𝑢𝑘
)+𝜈
(𝑈
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑡
𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘 𝑖 𝑘̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑗 𝜕 2 𝑝𝜕𝑢𝑗𝜈 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 3 𝑢𝑖=+ 𝜈2
𝜌 𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑗
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥𝑘𝜈(71)

同様に第２項は次式となる
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅𝑗𝜕𝑈𝜕𝑢𝜕𝑢𝑖 𝜕 2 𝑢𝑗
𝜕𝑢𝑖 𝜕
𝜕𝑢𝑖 𝜕
𝜕𝑢
𝜕2
̅𝑘 𝑗 ) + 𝜈 𝑖
+𝜈
(𝑢𝑘
)+𝜈
(𝑈
(𝑢 𝑢 − 𝑢𝑗 𝑢𝑘 )𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑡
𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘 𝑗 𝑘̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 3 𝑢𝑗𝜈 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕 2 𝑝
𝜕𝑢𝑖
=
+ 𝜈2
𝜌 𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑗
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥𝑘𝜈よって，非等方散逸テンソルに対する輸送方程式は次式で与えられる．

- 27 -(72)

JSME TED Newsletter, No.87, 2019̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑗 𝜕𝑢𝑖𝜕𝑢𝑗 𝜕𝑢𝑖
𝜕
𝜕
̅𝑘
(𝜈
)+𝑈
(𝜈
)
𝜕𝑡 𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅𝑗̅𝑖̅𝑘𝜕𝑢𝑗 𝜕𝑢𝑘 𝜕𝑈𝜕𝑢𝑗 𝜕𝑢𝑖 ̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑈𝜕𝑢𝑖 𝜕𝑢𝑗 𝜕𝑈
= −𝜈
−𝜈
−𝜈(
+
)
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ22̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅𝑗̅𝑖𝜕𝑢𝑗 𝜕 𝑈𝜕𝑢𝑗 𝜕𝑢𝑘 𝜕𝑢𝑖𝜕𝑢𝑖 𝜕 𝑈𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑗
−𝜈𝑢𝑘
− 𝜈𝑢𝑘
−𝜈
−𝜈
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 𝜕𝑢𝑗 𝜕 𝜕𝑢𝑖
𝜕 𝑝
𝜕 𝜕𝑢𝑗
𝜕 𝜕𝑢𝑖
−2𝜈 2
(
)
(
)+𝜈
( )[
( )+
( )]
𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ𝜕𝑥ℓ 𝜌 𝜕𝑥ℓ 𝜕𝑥𝑖
𝜕𝑥ℓ 𝜕𝑥𝑗̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑗
𝜕𝑢𝑗 𝜕𝑢𝑖𝜕 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑗
𝜕 𝑝
𝜕𝑢𝑖
𝜕2
−
[𝑢𝑘 𝜈 ()+( ) (𝛿𝑖𝑘 𝜈+ 𝛿𝑗𝑘 𝜈)] + 𝜈(𝜈)𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ
𝜕𝑥ℓ 𝜌
𝜕𝑥ℓ
𝜕𝑥ℓ
𝜕𝑥𝑘 𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥ℓ(73)̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
乱流エネルギーの散逸率𝜀は，𝜀 = 𝜀𝑖𝑖 = 𝜈(𝜕𝑢𝑖 ⁄𝜕𝑥ℓ )(𝜕𝑢𝑖 ⁄𝜕𝑥ℓ )であるので，その輸送方程式は次式となる．
̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅
̅𝑖
̅𝑘
̅𝑖
𝜕𝜀
𝜕𝜀
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑈𝜕𝑢𝑖 𝜕𝑢𝑖 𝜕𝑈𝜕𝑢𝑖 𝜕 2 𝑈
̅𝑘
+𝑈
= −2𝜈
− 2𝜈
− 2𝜈𝑢𝑘𝜕𝑡𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
′̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑖 𝜕𝑢2𝜈 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 𝜕𝑢𝑖 𝜕𝑝𝜕2𝜀
𝜕 2 𝑢𝑖𝜕 2 𝑢𝑖
𝑘𝜀
−2𝜈
−
−
(
)+𝜈
− 2𝜈 2 ()()
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘𝜕𝑥𝑘𝜌 𝜕𝑥𝑖 𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑘 𝜕𝑥𝑘𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ(74)

ここで，𝜀 ′ = 𝜈(𝜕𝑢𝑖 ⁄𝜕𝑥ℓ )(𝜕𝑢𝑖 ⁄𝜕𝑥ℓ )である．この輸送方程式において，両辺の各項はそれぞれ，
左辺第１項から時間項，対流項，右辺第１項から第４項までは生成項，第５項は乱流拡散項，第６項は圧力拡散項，第７項は粘性拡散項，そして最終項は消滅項（Destruction term）である．この輸送方程式を解けば乱流エネルギーの散逸率𝜀の分布が得られることになるが，粘性拡散項以外
すべての項のモデル化が必要である．この散逸方程式に対しては，各項が DNS から算出されたデータを再現するモデル化もされているが（Rodi and Mansour, 1993, Nagano and Shimada, 1995），
ここでは様々な手法によるモデル化を紹介する．
，
まず，右辺第１項の生成項（混合生成項（Mixed production）と呼ぶ（Rodi and Mansour, 1993））
に対しては，非等方テンソルモデル（Anisotropic tensor model）を用いる．このモデルは，次式で示すように，非等方散逸テンソル𝜀𝑖𝑗 と等方散逸率𝜀の比を，レイノルズ応力の非等方成分と等方成分で近似するものである．
𝜈̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ̅̅̅̅̅̅̅̅ =𝜕𝑢 𝜕𝑢𝜈 𝜕𝑥 𝑖 𝜕𝑥 𝑖ℓ̅̅̅̅̅̅̅̅̅̅
𝜀𝑖𝑘 𝑢𝑖 𝑢𝑘
𝜕𝑢𝑖 𝜕𝑢𝑘
𝑢𝑖 𝑢𝑘
~
⇒ 2𝜈
=𝜀
𝜀
𝑢𝑖 𝑢𝑖
𝜕𝑥ℓ 𝜕𝑥ℓ𝑘(75)ℓ

ここで，𝑢𝑖 𝑢𝑖 = 2𝑘である．右辺第２項（平均速度勾配生成項（Production by mean velocity gradient））
については，このモデルが適用できないため，第１項のモデルに含めてモデル化する（第２項を第１項を含めるには疑問が残ったが， DNS から算出されたそれぞれの項のオーダーを比較すると，ほぼ同じオーダーであったため，第２項を含めることには問題がない）．よって，これらの項は以下のようにモデル化される．
−2𝜈

̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅𝑖̅𝑘̅𝑖
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑈𝜕𝑢𝑖 𝜕𝑢𝑖 𝜕𝑈
𝜀
𝜕𝑈
− 2𝜈
= −𝐶𝜀1𝑓𝜀1 𝑢𝑖 𝑢𝑘𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ
𝑘
𝜕𝑥𝑘

(76)

ここで，𝐶𝜀1 はモデル定数，𝑓𝜀1 はモデル関数である．次に右辺第３項（勾配生成項（Gradient
Production））については，勾配拡散型モデルを用いてモデル化する．ここで，式(51)で示した乱流エネルギー輸送方程式における乱流拡散項のモデル化に用いた勾配拡散型モデルについて，Prandtl の混合距離理論を用いて一般化することを説明する．
𝑥𝑗 方向に勾配を持つ変数𝛷について，- 28 -

JSME TED Newsletter, No.87, 2019
̅を𝛷，
その変動成分を𝜙，混合距離をℓ𝑗 とすると，混合距離理論により，次式が成り立つ（図 3 の𝑈𝑦を𝑥𝑗 ，ℓをℓ𝑗 とされたい）．
𝛿𝛷 = 𝛷(𝑥𝑗 + ℓ𝑗 ) − 𝛷(𝑥𝑗 ) ≃

𝜕𝛷
ℓ
𝜕𝑥𝑗 𝑗

(77)

この差𝛿Φが変動成分と考えられるため，𝛿𝛷 = 𝜙 ≃ (𝜕𝛷⁄𝜕𝑥𝑗 )ℓ𝑗 となる．よって，変動速度𝑢𝑖 と変動成分の相関値は以下のようにモデル化される．
̅̅̅̅̅
̅̅̅̅̅
−𝑢
𝑖 𝜙 ∝ 𝑢𝑖 ℓ𝑗

𝜕𝛷
𝜕𝑥𝑗

(7 )

ここで，乱れの時間スケールτを導入すると，ℓ𝑗 = 𝑢𝑗 τより，次式を得る．これが一般化された
勾配拡散型モデル（Generalized Gradient Diffusion Hypothesis, GGDH（一般形勾配拡散仮説））
である．
̅̅̅̅̅
̅̅̅̅̅
−𝑢
𝑖 𝜙 = −𝐶𝑢𝑖 ℓ𝑗

𝜕𝛷
𝜕𝛷
= −𝐶𝑢
̅̅̅̅̅τ
𝑖 𝑢𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗

(79)

また，𝑢
̅̅̅̅̅を非等方速度スケールとみなし，非等方渦拡散係数を(𝜈̅̅̅̅̅τと定義すると，
𝑖 𝑢𝑗𝜙 ) = 𝐶𝑢𝑖 𝑢𝑗𝑖𝑗

渦拡散係数を用いた一般化された勾配拡散型モデルを得る．
̅̅̅̅̅−𝑢𝑖 𝜙 = −(𝜈𝜙 )

𝜕𝛷
𝑖𝑗 𝜕𝑥
𝑗

( 0)

よって，右辺第３項は次式のようにモデル化される．
̅̅̅̅̅̅̅̅
̅𝑖
̅𝑖
̅𝑖
𝜕𝑢𝑖 𝜕 2 𝑈
𝜕 𝜕𝑈
𝜕2𝑈
−2𝜈𝑢𝑘
= 𝐶𝜀3 𝜈(𝜈𝜙 )
(
)(
)
𝑘 𝜕𝑥
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘( 1)

ここで，𝐶𝜀3はモデル関数であり，𝑘 − 𝜀モデルでは(𝜈𝜙 )

𝑘

(𝜈𝜙 )

の添字が𝑘 = 𝑚，τ = 𝑘⁄𝜀であるため，= 𝐶𝑢̅̅̅̅̅̅̅̅τ𝑢 = 𝐶 𝑘 2⁄𝜀 となる．また，式(52)で導入された乱流プラントル数の概念を用いると，(𝜈𝜙 )

= 𝜈𝑡 ⁄𝜎𝜙 となり，モデル定数𝜎𝜙 を𝐶𝜀3に含めれば，式(81)は次式で与えられる．

̅̅̅̅̅̅̅̅
̅𝑖
̅𝑖
̅𝑖
𝜕𝑢𝑖 𝜕 2 𝑈
𝜕2𝑈
𝜕2𝑈
−2𝜈𝑢𝑘
= 𝐶𝜀3 𝜈𝜈𝑡 ()()
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝜕𝑥ℓ 𝜕𝑥𝜕𝑥ℓ 𝜕𝑥( 2)

もう一つの生成項である，右辺第４項は乱流生成項（Turbulent Production）と呼ばれ，非等方テ
ンソルモデルを用いて次式のようにモデル化される．
−2𝜈

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑖
𝜀 ̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖
𝜀 ̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖2𝜀 𝜕 ̅̅̅̅̅̅̅2= −𝐶𝜀4 𝑢𝑖 𝑢𝑘= −𝐶𝜀4 𝑢𝑘= −𝐶𝜀4𝑢 𝑢
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘𝑘𝜕𝑥𝑘
2𝑘 𝜕𝑥𝑘2𝑘 𝜕𝑥𝑘 𝑘 𝑖( 3)2̅̅̅̅̅̅̅
ここで，𝐶𝜀4はモデル関数であり，𝑢𝑘 𝑢𝑖 に勾配拡散型モデルを適用すると，次式となる．

̅̅̅𝑖2
𝜕𝑢
𝜕𝑘
̅̅̅̅̅̅̅
𝑢𝑘 𝑢𝑖2 = −𝐶τ𝑢̅̅̅̅̅̅
= −2𝐶τ𝑢̅̅̅̅̅̅
𝑘 𝑢𝑗
𝑘 𝑢𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗

( 4)

よって乱流生成項は次式のようにモデル化される．

- 29 -

JSME TED Newsletter, No.87, 2019−2𝜈̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑖𝜀 𝜕𝜕𝑘= −𝐶𝜀4(τ𝑢̅̅̅̅̅̅
)𝑘 𝑢𝑗
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝑘 𝜕𝑥𝑘
𝜕𝑥𝑗

( 5)

式(82)と同様に，式(85)は𝑘 − 𝜀モデルでは次式となる．
−2𝜈

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑖
𝜀 𝜕 𝜈𝑡 𝜕𝑘= −𝐶𝜀4()
𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘
𝑘 𝜕𝑥𝑘 𝜎𝑘 𝜕𝑥𝑘( 6)

乱流拡散項については，勾配拡散型モデルを用いて次式のようにモデル化される．
−′
̅̅̅̅̅̅𝜕𝑢𝜕 𝜈𝑡 𝜕𝜀
𝑘𝜀
=
(
)
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜎𝜀 𝜕𝑥𝑘( 7)

圧力拡散項は，まず非等方テンソルモデルを用いて次式のようにモデル化される．
−

2𝜈 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 𝜕𝑢𝑖 𝜕𝑝𝜕 𝜀 ̅̅̅̅̅
𝑝() = 𝐶𝜀5( 𝑢𝑖 )
𝜌 𝜕𝑥𝑖 𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑖 𝑘 𝜌()

ここで，𝐶𝜀5はモデル関数であり，圧力の変動成分𝑝に対し，そのポアソン方程式から速度成分に置き換えることが出来るため，勾配拡散型モデルと併用して𝑘 − 𝜀モデルでは次式となる．
𝜕𝑢̅̅̅̅̅
2𝜈 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 𝜕𝑢𝑖 𝜕𝑝𝜕 𝜀 ̅̅̅̅̅
𝑝
𝜕 𝜀
1 𝜕𝑃̅
𝜕
𝑖 𝑢𝑗
(
) = 𝐶𝜀5( 𝑢𝑖 ) = 𝐶𝜀5( ̅̅̅̅̅τ𝑢𝑖 𝑢𝑗) = 𝐶𝜀5(𝑘)𝜌 𝜕𝑥𝑖 𝜕𝑥ℓ 𝜕𝑥ℓ𝜕𝑥𝑖 𝑘 𝜌
𝜕𝑥𝑖 𝑘
𝜌 𝜕𝑥𝑗
𝜕𝑥𝑖
𝜕𝑥𝑗
̅𝑖
𝜕
𝜕
𝜕𝑈
= −𝐶𝜀5[𝑘(𝜈𝑡)]𝜕𝑥𝑖 𝜕𝑥𝑗
𝜕𝑥𝑗

−

( 9)

最後に，消滅項のモデル化であるが，この項のような変動の２階微分に対するモデル化技術が確立していないため，次元解析によるモデル化を行う．
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 2 𝑢𝑖𝜕 2 𝑢𝑖𝑉 2 𝑉4−2𝜈 2 ()() ≃ (𝑉𝐿 ∙ 2 ) = 4
𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ
𝐿
𝐿

(90)

これより，渦拡散係数のモデル化で用いた乱れの速度スケール𝑉 = √𝑘，乱れの長さスケール𝐿 =𝑘 3⁄2⁄𝜀を代入すると，次式を得る．
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕 2 𝑢𝑖𝜕 2 𝑢𝑖𝜀2−2𝜈 2 ()(
) = −𝐶𝜀2𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ𝑘(91)

ここで，𝐶𝜀2はモデル関数である．以上，モデル化された乱流エネルギーの散逸率𝜀の輸送方程
式は次式で与えられる．
̅𝑖
̅𝑖 2
𝜕𝜀
𝜕𝜀
𝜕
𝜈𝑡 𝜕𝜀
𝜀
𝜕𝑈
𝜀2
𝜕2𝑈
̅𝑗
+𝑈
=
[(𝜈 + )
] − 𝐶𝜀1𝑓𝜀1 𝑢𝑖 𝑢𝑗− 𝐶𝜀2𝑓𝜀2 + 𝐶𝜀3 𝜈𝜈𝑡 ()
𝜕𝑡𝜕𝑥𝑗 𝜕𝑥𝑗
𝜎𝜀 𝜕𝑥𝑗𝑘𝜕𝑥𝑗𝑘
𝜕𝑥ℓ 𝜕𝑥(92)

ここで，生成項と消滅項には補正モデル関数𝑓𝜀1 と𝑓𝜀2 が付加されている．また，式(86)でモデル化された乱流生成項，式(89)でモデル化された圧力拡散項は，通常無視される（𝐶𝜀4 = 𝐶𝜀5 = 0と考えても良い）．
以上のようにモデル化技術を駆使して乱流エネルギー散逸率の輸送方程式をモデル化したが，最近の研究（Rodi and Mansour, 1993, Nagano and Shimada, 1995）では，生成項は右辺第１項と第２項から導出されるのではなく，上記の説明では無視された第３項の乱流生成項から導出される．
ただし，その場合には消滅項も同時にモデル化し，以下のように表現される．

- 30 -

JSME TED Newsletter, No.87, 2019−2𝜈̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅𝑖𝜕𝑢𝑖 𝜕𝑢𝑘 𝜕𝑢𝑖
𝜕 2 𝑢𝑖𝜕 2 𝑢𝑖𝜀𝜕𝑈𝜀2− 2𝜈 2 ()() = −𝐶𝜀1 𝑓𝜀1 𝑢𝑖 𝑢𝑗− 𝐶𝜀2 𝑓𝜀2𝜕𝑥ℓ 𝜕𝑥ℓ 𝜕𝑥𝑘𝜕𝑥𝑘 𝜕𝑥ℓ 𝜕𝑥𝑘 𝜕𝑥ℓ
𝑘
𝜕𝑥𝑗
𝑘

(93)

前述のモデリングでは，直感的に平均速度勾配がある項が主要生成項であると考えていたが，DNS の結果により，これらの項は寄与が小さいことが分かり，式(93)で与えられるモデリングがより妥当である結果となっている（Rodi and Mansour, 1993, Nagano and Shimada, 1995）．
ところで，非常に難しいモデル化を施して乱流エネルギー散逸率の輸送方程式をモデル化したが，実のところ，もっと簡単にモデル化する方法がある．𝑘 − 𝜀モデルの渦拡散係数の表現に用いた乱れの時間スケール𝜏 = 𝑘⁄𝜀の実質微分をとると次式となる．
𝐷𝜏 𝐷 𝑘
1 𝐷𝑘 𝑘 𝐷ε=( )=−𝐷𝑡 𝐷𝑡 εε 𝐷𝑡 ε2 𝐷𝑡(94)

これより，時間スケールの時間微分は𝐷𝜏⁄𝐷𝑡 = 0であることから，散逸率に対する実質微分を得る．
𝐷ε ε 𝐷𝑘=
𝐷𝑡 𝑘 𝐷𝑡(95)

つまり，次式で与えられる𝑘 − 𝜀モデルで使用される乱流エネルギーの輸送方程式（式(55)の散逸項をεと置く）に𝜀 ⁄𝑘を掛けた式であることは，式(92)と比較すれば自明である（式(92)にある，勾配生成項はこのモデリング方法では導出されないが，この項が必要であれば付加すればよい）．
̅𝑖
𝐷𝑘
𝜕
𝜈𝑡 𝜕𝑘
𝜕𝑈
=
[(𝜈 + )
] − 𝑢𝑖 𝑢𝑗−𝜀𝐷𝑡 𝜕𝑥𝑗𝜎𝑘 𝜕𝑥𝑗𝜕𝑥𝑗(96)

以上により，𝑘 − 𝜀モデルにおける乱流エネルギーとその散逸率に対する，モデル化された輸送方程式を導出したが，様々なモデル定数や，モデル関数が出現した．それらを決定するにも何らかの手法が必要であるため，以下ではその説明を行う．
まず，運動量の渦拡散係数にあるモデル定数𝐶𝜇 である．𝐶𝜇 の定義は，上述にもあるが，構造パラメーターとして次式で与えられる．
−𝑢𝑣
̅̅̅̅
= √𝐶𝜇
𝑘

(97)

0.5

−𝑢𝑣⁄𝑘

0.4
𝐶𝜇 ~0.3

0.3
0.2

DNS
BL:
BL:
CH:
CH:

0.1
0
0

0.4
Fig. 10 Stractual parameters- 31 -= 1000
= 300
2,
= 395
= 590

2,

𝑦⁄ 𝛿99 , 𝑦⁄𝛿 0.8

JSME TED Newsletter, No.87, 2019̅ ⁄𝑑𝑦) = 𝐶𝜇 (𝑘 2⁄𝜀 )(𝜕𝑈
̅ ⁄𝑑𝑦)
この関係式は，レイノルズ応力と渦拡散係数のモデル式−𝑢𝑣̅̅̅̅ = 𝜈𝑡 (𝜕𝑈からも成立する．乱流境界層の壁面近くに定応力層（Constant stress layer）が存在し，そこでは，
̅ ⁄𝑑𝑦) = 𝜀）．この関係を，レイノルズ応力と渦拡散係数のモデ
局所平衡が成立している（−𝑢𝑣̅̅̅̅(𝜕𝑈
ル式に代入し，𝐶𝜇 について求めると，式(97)と同じ関係式を得ることから，𝐶𝜇 を渦拡散係数のモデル定数とすることには整合性が取れていることが理解できる．
̅ ⁄𝑑𝑦) −𝑢𝑣
𝜀 −𝑢𝑣
̅̅̅̅
−𝑢𝑣
̅̅̅̅(𝜕𝑈̅̅̅̅
−𝑢𝑣̅̅̅̅ 2==()̅ ⁄𝑑𝑦)
̅ ⁄𝑑𝑦)(𝜕𝑈
𝑘 2 (𝜕𝑈
𝑘2
𝑘

𝐶𝜇 =

(9 )

この𝐶𝜇 の分布について，乱流境界層の DNS（Hattori et al., 2007, 2016）と，完全発達した平行平
板間乱流の DNS（Moser, et al., 1997）から算出した値を図 10 に示す．壁面極近傍では値が小さくなっているが，どのデータもほぼ 0.3 付近に分布していることが分かる．そのため，この実験的事実から，𝐶𝜇 は 0.09 で与えられる（提案された時代には DNS ではなく，装置実験による測定データから同定されたが，DNS も同様の結果を与えていることから，妥当な値である）．
次に，乱流エネルギー散逸率のの輸送方程式に関するモデル定数の決定法について記す．式(92)において，定応力層では次式が成立する．
0=̅𝜕 𝜈𝑡 𝜕𝜀𝜀𝜕𝑈𝜀2() − 𝐶𝜀1 𝑢𝑣− 𝐶𝜀2
𝜕𝑦 𝜎𝜀 𝜕𝑦
𝑘
𝜕𝑦
𝑘

(99)

̅ ⁄𝑑𝑦)）
ここで，式(57)の関係式−𝑢𝑣 = 𝑢𝜏2 ，０方程式モデルから導出される関係式(36)（𝑢𝜏 = 𝜅𝑦(𝜕𝑈，
上述の乱流エネルギー輸送方程式における生成項と散逸項の釣合いを用いると，散逸率に対する
次の関係式を得る．
−𝑢𝑣

̅
̅
𝜕𝑈
𝜕𝑈
𝑢𝜏3
= 𝜀 ⇒ 𝜀 = 𝑢𝜏2
⇒ 𝜀=
𝜕𝑦
𝜕𝑦
𝜅𝑦

(100)

構造パラメーターの式(97)より，𝑘 = 𝑢𝜏2 ⁄√𝐶𝜇 であるから，これらの関係を式(99)に代入すると，次式を得る．
2

𝐶𝜇 𝜕 𝜅𝑦 4 𝜕 𝑢𝜏3𝑢𝜏3 √𝐶𝜇 𝑢𝜏3𝑢𝜏3 √𝐶𝜇[( 3 ) 𝑢𝜏( )] + 𝐶𝜀1
−
𝐶
(
)
𝜀2
𝜎𝜀 𝜕𝑦 𝑢𝜏𝜕𝑦 𝜅𝑦𝜅𝑦 𝑢𝜏2 𝜅𝑦
𝜅𝑦
𝑢𝜏2
𝐶𝜇 1
√𝐶𝜇
√𝐶𝜇
0=
+
𝐶
−
𝐶
𝜀1
𝜀2
(𝜅𝑦)2
(𝜅𝑦)2
𝜎𝜀 𝑦 20=(101)

式を整理すると，モデル定数の関係式を得る．
𝐶𝜀1 = 𝐶𝜀2 −

𝜅2

(102)

𝜎𝜀 √𝐶𝜇次に，図 9 で示したような減衰乱流を考えると，乱流エネルギーとその散逸率の輸送方程式は次の式となる．
𝑑𝑘
= −𝜀,
𝑑𝑡

𝑑𝜀
𝜀2
𝑑2𝑘
1 𝑑𝑘 2= −𝐶𝜀2 𝑓𝜀2 ⇒= 𝐶𝜀2 𝑓𝜀2 ( )
𝑑𝑡
𝑘
𝑑𝑡 2
𝑘 𝑑𝑡

(103)

減衰乱流では，実験的に減衰則としてに𝑘 = 𝐶𝑡 −𝑛 の関係となることが知られている（Bachelor
and Tawnsend, 1948）ため，これを式(103)に代入すれば，𝐶𝜀2 𝑓𝜀2 = (𝑛 + 1)⁄𝑛を得る．減衰初期では，𝑛 = 1 ∼ 1.25，𝑓𝜀2 = 1のため，𝑛 = 1とすると，𝐶𝜀2 = 2を得る．𝜎𝜀 を 1 とすると，式(102)より，
𝐶𝜀1 = 1.47を得る（𝜅 = 0.4とした）．標準𝑘 − 𝜀モデルでは，𝐶𝜀1 = 1.44，𝐶𝜀2 = 1.92，𝜎𝜀 = 1.3の値- 32 -

JSME TED Newsletter, No.87, 2019が取られている．また，乱流エネルギーの輸送方程式におけるモデル定数𝜎𝑘 は，経験的に𝜎𝑘 = 1とされるが，応力方程式モデルの研究成果により（Hanjalić and Launder, 1976），𝜎𝑘 ⁄𝜎𝜀 = 1.09が提案されたため，𝜎𝑘 = 1.4と置くモデルもある（Nagano and Tagawa, 1990）．モデル定数は，数々提案された𝑘 − 𝜀モデルで，すべて同じ値ではなく，計算対象とした乱流場に対して数値最適化されている場合が多いが，上記に示した関係式を大きく逸脱していない．これは，ある複雑乱流場の計算において，予測結果が上記の関係式を逸脱したモデル定数を用いた時に良好な結果を示したとしても，その結果をもって乱流モデル定数の最適値というべきでないことを，多くの研究成果が語っている．
次にモデル関数について述べる．０方程式モデルの節で述べたが，壁面近くの乱流挙動の変化を取り入れてモデルを補正しないと，モデルは良好な予測値を与えない．一方，標準𝑘 − 𝜀モデル
は，上記で導出した対数速度分布といった定応力層における関係式を，壁面から少し離れた場所に人工的な壁面境界条件として用いるため，壁面近くの補正は不要であった（これを高レイノル
ズ数形乱流モデル（High-Reynolds-number turbulence model）と呼ぶ）．これに対して，壁面で物理的な境界条件（流れ場においては粘着条件）を用いて，壁面から解くモデルを低レイノルズ形乱流モデル（Low- Reynolds-number turbulence model）と呼ぶ．高レイノルズ数形乱流モデルの計算方法については，第３回のニュースレターで詳しく述べるが，ここでは，低レイノルズ数形𝑘 − 𝜀
モデルのモデル関数について述べる．
図 10 で示した構造パラメーターが，壁面近くで急激に減少していることが分かるが，このことから𝑘 − 𝜀モデルにおいて，渦拡散係数のモデル形𝜈𝑡 = 𝐶𝜇 (𝑘 2⁄𝜀)が，𝐶𝜇 を一定値としているため，𝑘 2⁄𝜀がこの変化をカバーすればモデル関数を付加する必要はない．しかしながら，そこまで上手
く変化しないため，０方程式モデルで導入されたように壁面影響関数𝑓𝜇 を次式のように取り入れた（Jones and Launder, 1972）𝜈𝑡 = 𝐶𝜇 𝑓𝜇

𝑘2
𝜀

(104)

ここで，壁面影響関数𝑓𝜇 の関数形は，最初に提案された低レイノルズ形𝑘 − 𝜀モデル（Jones andLaunder, 1972，
JL モデル）では，
壁面近くで乱流レイノルズ数 𝑡 = 𝑘 2⁄𝜈𝜀（Turbulent Reynolds number）
が小さくなることと，乱流から層流へ遷移した時に 𝑡 = 0となっても発散しないように，次の関数形を用いた．
𝑓𝜇 = exp (−2.5)
1 + 𝑡 ⁄50

(105)

このことから，低レイノルズ数形の低レイノルズ数というのは，内部流れで用いられる流れのバルクレイノルズ数や，境界層などで用いられる運動量厚さによるレイノルズ数が低いという意
味ではなく，乱流レイノルズ数が低いことということが分かる（ただし，乱流レイノルズ数は，境界層の外縁においても低い値を取るため，筆者の解釈は，無次元距離𝑦 + = 𝑦𝑢𝜏 ⁄𝜈が小さい場所に適用できるモデルと考えている．それは，無次元距離𝑦 +の関数形を見ると，実はレイノルズ数になっていることからも考えられる）．
さらに JL モデルでは，乱流エネルギーの散逸率のモデル輸送方程式(92)におけるモデル関数を，
𝑓𝜀1 は 1 と置くが，𝑓𝜀2 については，上記の減衰乱流の末期を再現するために，以下の関数形で与えている．
𝑓𝜀2 = 1 − 0.3exp(− 𝑡 )(106)

さらにさらに，乱流エネルギーの散逸率の壁面境界条件は，厳密には𝜀𝑤 = ν(𝜕 2 𝑘⁄𝜕𝑦 2 )|𝑤 であるが，これを与えるには壁面上で困難であるため，乱流の壁面漸近論（Wall limiting behavior of
turbulence）を用いて，𝜀𝑤 = 0と出来るようモデル化されている．そのため，乱流エネルギーとその散逸率の輸送方程式は次式で与えられる．

- 33 -

JSME TED Newsletter, No.87, 2019

DNS (Abe & Kawamura)
JL modelSA model̅+𝑈20

̅ + = 1⁄ 𝜅𝑈𝑦+ + 510= 24420 0
10102 𝑦 +

101103(a) Mean velocity

6𝑘+

DNS (Abe & Kawamura)
JL model3= 2442

0 0
10

101

102 𝑦 +

103(b) Turbulence kinematic energy
Fig. 11 Prediction of turbulent channel flow using LRN 𝑘 − 𝜀 model

̅𝑖
𝐷𝑘
𝜕
𝜈𝑡 𝜕𝑘
𝜕𝑈
=
[(𝜈 + )
] − 𝑢𝑖 𝑢𝑗−𝜀𝐷𝑡 𝜕𝑥𝑗𝜎𝑘 𝜕𝑥𝑗
𝜕𝑥𝑗

(107)

̅𝑖
̅𝑖 2
𝜕𝜀̃
𝜕𝜀̃
𝜕
𝜈𝑡 𝜕𝜀̃
𝜀̃
𝜕𝑈
𝜀̃ 2
𝜕2𝑈
̅𝑗
+𝑈
=
[(𝜈 + )
] − 𝐶𝜀1𝑓𝜀1 𝑢𝑖 𝑢𝑗− 𝐶𝜀2𝑓𝜀2 + 2𝜈𝜈𝑡 ()
𝜕𝑡𝜕𝑥𝑗 𝜕𝑥𝑗
𝜎𝜀 𝜕𝑥𝑗
𝑘
𝜕𝑥𝑗
𝑘
𝜕𝑦𝜕𝑦

(10 )

ここで，𝜀̃はその壁面境界条件を 0 とするための，疑似散逸率である．乱流エネルギーの輸送方
程式(107)に対する散逸率𝜀については，実質散逸率を与えるため，補正項を用いて次式で与える．
2𝜕√𝑘𝜀 = 𝜀̃ + 2𝜈 (
)
𝜕𝑦

(109)

なお，この補正項は壁面近くの𝜕√𝑘⁄𝜕𝑦 > 0となる箇所まで付加し，その後は付加しない．また，
渦拡散係数と乱流レイノルズ数の算出には𝜀̃を用いているため，この付加項と合わせて計算には注- 34 -

JSME TED Newsletter, No.87, 2019
意が必要である．また，JL モデルにおけるモデル定数は，𝐶𝜀1 = 1.55，𝐶𝜀2 = 2，𝜎𝑘 = 1，𝜎𝜀 = 1.3である．
図 11 に，JL モデルによる予測した，完全発達した平行平板間乱流の平均速度と乱流エネルギ
ーの分布を示す．比較のため，１方程式モデルである SA モデルによる平均速度の予測結果を，図に入れてある．JL モデルの予測結果は，遷移域で若干低くなっているが，全体を通して良好な結果を与えている．乱流エネルギーの予測値は，DNS の結果と比較して最大値付近で過少予測されていることが分かる．しかしながら，DNS と比較して悪いモデルであるというのではなく，
DNS が発達した現代のように，壁面近傍の乱流データが豊富ではない時代に開発されたことを考えると，それでもこれだけの予測値を与えるモデルであることに敬服するしかない．また，JL モ
デルが，その時代の最高の乱流モデルであったことは，この JL モデルをベースとして，多くの改良低レイノルズ数形モデルが提案されていったことからも言えるであろう．
第 3 回までの予定で，乱流熱伝達モデル入門と題しながら，この原稿を執筆しているが，第２回目は速度場モデルの話題に終始してしまった．最終回の第 3 回は，タイトル通りに，熱伝達（温度場）モデルについて述べるが，速度場モデルについても，実はまだ書き切れていない話題が沢山ある．Web 掲載のため，ページ数は好きなだけ書いてよいと言われているので，次回，最終回
は，基本事項は抑えつつ，執筆時間が許すだけ筆者が経験してきた乱流モデル開発の話題にも触れる予定である．

文献
Abe, H. and Kawamura, H., DNS of turbulent heat transfer in a channel flow up to𝜏 = 640 (Withrespect to Reynolds- and Prandtl-number effects), Transactions of the JSME (B), Vol.69, pp.2291-2298,
2003 (in Japanese).
Bachelor, G. K. and Tawnsend, A. A., Decay of turbulence in the final period, Proc. Roy. Soc. London, Vol.
A194 (1948), pp. 527-543.
Baldwin, B.S. and Barth, T.J., A one-equation turbulence transport model for high Reynolds number
wall-bounded flows, NASA TM 102847 (1990).
Cebeci, T., Variation of the Van Driest damping parameter with mass transfer, AIAA Journal, Vol. 11 (1973),
p. 237.
Hanjalić, K. and Launder, B. E, Contribution towards a Reynolds-stress closure for low-Reynolds-number
turbulence, Journal of Fluid Mechanics, Vol. 74 (1976), pp. 593-610.
Harlow, F. H. and Nakayania, P. I., Turbulence transport equations, Physics of Fluids, Vol. 10, No. 11(1967), p. 2323.
Hattori, H., Houra, T. and Nagano, Y., Direct numerical simulation of stable and unstable turbulent thermalboundary layers, International Journal of Heat and Fluid Flow, Vol. 28 (2007), pp. 1262-1271.
Hattori, H., Kono, A. and Houra, T., Direct numerical simulation of thermally-stratified turbulent boundary
layer subjected to adverse pressure gradient, International Journal of Heat and Fluid Flow, 61A (2016), pp.
213-228.
Hishida, M., Nagano, Y. and Tagawa, M., Transport processes of heat and momentum in the wall region of
turbulent pipe flow, Proceedings of 8th International Heat Transfer Conference, San Francisco, Vol. 3 (1986),
pp. 925-930, Hemisphere.
Jones, W. P. and Launder, B. E., The prediction of laminarization with a two-equation model of turbulence,
International Journal of Heat and Mass Transfer, Vol. 15 (1972), pp. 301-314.
Kays W. M. and Crawford, M. E., Convective heat and mass transfer, International editions 1993,
McGraw-Hill.
Kolmogorov, A. N., Equations of turbulent motion of an incompressible fluid. Izv. Akad. Nauk SSSR Ser.
Fiz. 6 (1942), pp. 56-58
Launder B. E. and Spalding, D. E., The numerical computation of turbulent flows, Comp. Meths. Appl.
Mech. Engrng, Vol. 3 (1974), pp. 269-289.
Menter, F. R., Eddy viscosity transport equation and their relation to the 𝑘 − 𝜀 model, Journal of Fluids
Engineering, Vol. 119 (1997), pp. 876-884.
Moser, R. D., Kim, J. and Mansour, N. N., Direct numerical simulation of turbulent channel flow up to𝜏 = 590, Physics of Fluids, Vol. 11, No.4 (1999), pp. 943-945.
Nagano, Y. and Hattori, H., Direct numerical simulation and modelling of spanwise rotating channel flow

- 35 -

JSME TED Newsletter, No.87, 2019with heat transfer, Journal of Turbulence, Institute of Physics, 4-010 (2003).
Nagano Y. and Shimada, M., Rigorous modeling of dissipation-rate equation using direct simulations,
JSME International Journal, Series B, Vol. 38, No. 1 (1995), pp. 51-59.
Nagano Y. and Tagawa, M., An Improved k-ε Model for Boundary Layer Flows, ASME Journal of Fluids
Engineering, Vol. 112, No. 1 (1990), pp. 33-39.
Nee, V. W. and Kovasznay, L. S. G., Simple phenomenological theory of turbulent shear flows, Physics of
Fluid, Vol. 12 (1969), pp. 473-484.
Prandtl, L., Z. angew, Math. Mech. Vo. 5, No. 1 (1925), pp. 136-139.
Prandtl, L.: Über ein neues Formelsystem für die ausgebildete Turbulenz, Nachr. Akad. Wiss. Göttingen,Math.-Phys. Klasse, (1945), pp. 6–19.
Rodi, W. and Mansour, N. N., Low Reynolds number 𝑘 − 𝜀 modelling with the aid of direct simulation
data, Journal of Fluid Mechanics, Vol. 250 (1993), pp. 509-529.
Rodi, W. and Spalding, D. B., A two-parameter model of turbulence and its application to free jets,
Wärme-und Stoffübertragung Vol. 3, No. 2 (1970), pp. 85-95.
Rotta, J. C., Statistische theorie nichthomogener turbulenz, Zeitsch, für Physik Vol. 129 (1951), p. 547 andVol. 131 (1951), p. 51.
Saffman, P. G., A model for inhomogeneous turbulent flow, Proceedings of the Royal Society of London,
Series A: Mathematical and Physical Sciences, Vol. 317, No. 1530 (1970), pp. 417–433.
Spalart, P. R. and Allmaras, S. R., A one-equation turbulence model for aerodynamics flows, La Recherche
Aerospatiale, No. 1 (1994), pp. 5-21.
Spalding, D. B., The vorticity-fluctuations (KW) model of turbulence, early papers’ Imperial College
CFDU Rep. CFD/82/17 (1982).
Van Driest, E. R., On turbulent flow near a wall, Journal of the Aeronautical Sciences, Vol. 23, No.11(1956), pp. 1007-1011.
Vogel, J. C. and Eaton, J.K., Combined heat transfer and fluid dynamic measurements downstream of a
backward-facing step, Transaction of ASME, Journal of Heat Transfer, Vol. 107 (1985), pp. 922-929.
von Kármán, Th, Mechanische Ähnlichkeit und Turbulenz, Nachr. Ges. Wiss Göttingen Math. Phys. Kasse,
Vol.58 (1930) and Proceedings of 3rd International Congress Appl. Mech., Stockholm Prat I, 85, NACA
TM 611 (1931).
Wilcox, D. C., and Alber, I. E., A turbulence model for high speed flows, Proceedings of the 1972 Heat
Transfer and Fluid Mechanics, Institute, Stanford Univ. Press, Stanford, CA (1972), pp. 231–252.
Wilcox, D. C., Reassessment of the scale determining equation for advanced turbulence models, AIAA
Journal, Vol. 26, No. 11 (1988), pp. 1299–1310.
Wilcox, D. C., Turbulence modeling for CFD, 2nd ed., DCW Industries, Inc., La Cañada, CA (1998).

- 36 -