JSME TED Newsletter, No.86, 2018

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

名古屋工業大学 課長
技術部情報解析技術課
hattori@nitech.ac.jp
１．はじめに熱や物質の輸送現象（Transport Phenomena of Heat and Mass Transfer）を物理的や工学的に正確かつ詳細に把握することは，伝熱問題，物質混合問題，気象問題，環境問題を解決するために必須事項である．その応用問題として，熱交換機器の熱交換効率問題や，機器や密閉空間からの排熱・換気問題，また異常熱による機器破損などの熱問題，化学薬品等の物質の混合効率問題，気象の風況・気温予測（大きく言えば天気予報）問題，都市部におけるヒートアイランド問題，これらをまとめたエネルギー問題等々があり，熱や物質輸送現象は多岐にわたる問題に絡んでいる．
さて，熱や物質が輸送される現象の一つとして対流現象があり，工業製品や気象，環境など様々な場面で出くわす現象である．
対流の形態には，浮力が影響しない強制対流場（Forced Convection），
浮力が支配的な自然対流場（Natural Convection），それらが混合している共存対流場（複合対流場：Combined Convection）がある．それゆえ，熱に関する「温度場」とともに「流れ場」の解析が重要となってくる．その流れ場であるが，層流と乱流状態があり，現実の現象としては乱流となっていることが多い．乱流はその定義として，「流れ場の無次元パラメーターであるレイノル
ズ数（Reynolds Number）が高い時の流れの運動方程式の解」，「3 次元的な渦度の変動がある」，「大きな混合と拡散（Mixing and Diffusion）を伴う」，そして，「不規則で偶然的な時間，空間における変動を伴う流れ」であるとされている．この定義の最後に「不規則で偶然的な時間，空間における変動を伴う流れ」とあるが，流れの運動方程式（Navier-Stokes Equation もしくは Equationof Momuntum）と，連続の式（Equation of Continuty）を連立させて正確に解くことが出来れば，
不規則で偶然的である乱流現象を捉えることは可能である（ということは，結果として不規則で偶然的ではないことになる）．このことは，バースト現象（Bursting Phenomena）やストリーク構
造（Streak Structure）に代表される乱流の組織構造や準秩序構造（Coherent Structure）の発見により，乱流とは「必然と偶然の中間的な現象である」とされた．乱流の組織構造は，実験的研究で解明されたが，支配方程式（微分方程式）を乱流の条件（レイノルズ数が高い時）で解析的に解くことは困難である（Navier-Stokes 方程式の解の存在と滑らかさについては，数学のミレニアム
問題（Millennium Prize Problems）の一つとなっており，数学的には一般解を求めることが要求されていると思われるが，初期条件が違えば層流解と乱流解が存在する時点で不可能ではないかと思われる？？）．そのため，計算機で数値的に解くこととなるが，時々刻々と変化し，様々なス
ケールを持つ 3 次元的な渦度の変動を数値的に解くには，計算機に多大な性能が要求される．
上述のような乱流の組織構造を発見した実験的手法では，乱流中における速度や圧力の空間分
布を計測することや，可視化することが行わている（有名な Osborne Reynolds の実験（1883）では，染料による円管内流れの可視化によって臨界レイノルズ数（Critical Reynolds Number）の存在を証明した）．乱流空間における速度成分の計測器には，ピトー管や熱線流速計（佐々木，玉野，
1941）などの接触式計測器，非接触式計測器として，レーザーを用いた流速計（LDV: Laser Dopplar- 14 -

JSME TED Newsletter, No.86, 2018
Velocimetory; Durst, et al., 1972 など）や，粒子画像流速計（PIV: Particle Image Velocimetory; Willert
and Gharib, 1990 など），超音波ドップラー計測器などがある．前者のうち，熱線流速計では流速の時間変化の計測が可能でかつ計測精度が高いが，基本的にポイント計測となるため空間計測には多くの手間を生じることや，乱流場で重要となる３次元的な渦度の変動を計測する場合や逆流を伴う場合にはその適用が難しい．また後者の PIV は，逆流を伴う流れ場を含め，様々な流れ場の空間的な計測が可能であるため，広く普及してきている．さて，この乱流場に熱や物質移動が伴うことは，現実的な場では普通である．そのため，流れ場に加え，熱や物質の空間分布を計測することが重要となる．熱に関しては，その温度分布を計測することとなるが，乱流場において
空間的に時間変動する温度分布を計測する方法は，熱電対や抵抗線を用いる方法が一般的であるが（例えば，田川他，1998），レーザーによる温度計測法もある（Nakajima, et al., 1990 など）．
乱流熱・物質輸送現象を解析するためには，速度変動の直交する成分の相関であるレイノルズせん断応力（Reynolds Shear Stress）や，速度変動と温度もしくは濃度などのスカラー物質変動の相関値である乱流熱流束もしくは乱流スカラー流束（Turbulent Heat Flux もしくは Turbulent ScalarFlux）を正確に計測する必要がある．特に乱流熱流束の計測では，速度場と温度場の瞬時場を同時計測する必要があるが，熱線流速計と抵抗線の組合せによる速度温度同時計測(菱田他，1977)，もしくは LDV による速度温度同時計測(Nakajima, et al., 1990)，LDV と熱電対の組合せによる速度
温度同時計測(田川他，1999)などの計測例がある．
本稿で扱う乱流熱・物質輸送現象は，境界条件の一つに「壁面」を持ち，そのため乱流境界層（Turbulent Boundary Layer）を形成する場における現象を扱う．また，より簡単に解説するため，作動流体の空間における密度と物性値を一定と考え，つまり非圧縮性流体と仮定して扱う．加熱された壁面上を流れる乱流境界層のように，壁面で熱伝達（Heat Transfer）を伴う場を解析するには，薄い乱流境界層内の更に薄い粘性底層（Viscous Sublayer）や伝導底層（Conductive Sublayer）
も含んだ内層領域（Inner Region）における乱流熱流束を計測しなければならないが，上記に挙げた計測方法でも困難であることが多い．そのため，数値的な解析手法で乱流境界層内まで解くことが出来れば，実験的手法で計測が困難な領域における乱流熱流束分布が取得されるため，実験による計測値と合わせて乱流境界層における乱流熱・物質輸送現象をより正確かつ詳細に把握することが出来るであろう．しかしながら，乱流境界層を数値的に正確に解くためには，膨大な計算機資源を必要とする．上記のような乱流現象を正確に解くためには，支配方程式による連立微分方程式を解くためのより精度の高い計算スキームと，適切な計算格子が必要である．より精度の高い計算スキームは，スペクトル法による方法(Kim, et al., 1987)があるが，ある方向の境界条件に周期境界条件を用いる必要があるため，適用に限界がある．また，より汎用性が高い差分法による高精度スキームが考案され(Morinishi, et al., 1998)，こちらの方法によると様々な壁面形状を持つ乱流境界層にも適用できる．これらは，熱・物質輸送現象を支配する方程式を，何ら物理的もしくは数学的な付加項を用いずに直接解くため，直接数値シミュレーション（ DNS: Direct
Numerical Simulation）と呼ばれ，速度場に限れば，乱流スケールの最小スケールであるコロモゴ
ロフスケール（Kolmogorov Scale）まで捉えるが，そのための必要な 3 次元空間における格子点数
𝑁は，代表速度𝑈と代表長さ𝐿，動粘性係数𝜈で定義される流れ場のパラメーターであるレイノルズ数(Re = 𝑈𝐿⁄𝜈)に対し，概ね𝑁~Re9⁄4に比例する．そのため，例えば人が歩いている時の周りの流れの DNS を想定すると，代表速度𝑈 = 1 [m/s]（時速約 4 km/h），代表長さ𝐿 = 0.25 [m]（胴囲を0.8m），空気の動粘性係数𝜈 = 1.512 × 10−5 [m2 ⁄s]（気温が 293.15K の時）であるので，レイノルズ数は約 17,000 となり，必要格子点数は約 33 億点である．この格子点数に，４変数（3 速度と圧力）を 2 時刻で配置する必要があり，倍精度計算とすると，その記憶容量は約 100 ギガバイト（Giga=109）必要である．これは現在市販されているワークステーションで計算可能な範囲であるが，計算速度を考えると，初期条件から時々刻々と変化する速度と圧力を計算しながら，レイ
ノルズ応力や乱流熱流束を算出するための十分な統計量（Tubulence Statistics）を取得するためには数カ月の時間が必要となるであろう．ここで統計量という用語を出したが，上述のレイノルズ
応力や乱流熱流束は，次章に示す「レイノルズ平均（Reynolds Average）」(Reynolds, 1895)という統計処理を行った結果生じる乱流統計量である．
上述の事柄は，乱流熱・物質輸送現象を物理的や工学的に正確かつ詳細に把握するための手段として記述したが，これらの事を行うには膨大な時間と費用が掛かり，またそれ相当の知識も必要となるため，非常に敷居が高い．例えば，加熱された円管内を流れ，下流で発達する乱流熱伝- 15 -

JSME TED Newsletter, No.86, 2018
達場（Turbulent Heat Transfer of Fully-Developed Pipe Flow）の乱流統計量を得るとしよう．実験で統計量を得るためには，装置の開発から，測定装置の準備が必要である（場合によっては，測定装置の開発も必要となる）．これらの準備が終わったとしても，有意な計測結果を出すためには数日を有することであろう．では，DNS ではどうか？DNS の経験があれば，計算格子と境界条件の設定，初期値の準備が終われば計算に取り掛かることが出来るが，経験がなければ初期値の準備方法からして大きな問題となる．さらに適切な条件を設定しないと，計算も収束するとは限らない．また，両者に共通する事項として，狙ったレイノルズ数における実験装置の大きさ（格子の大きさ），計測装置の空間分解能（壁面近くの格子幅）等，装置や DNS を設計する上で考えなければならないことが幾つもある．そのため，測定対象の乱流熱伝達場からどんな乱流統計量が得られるかが「予測（Prediction）」出来れば，その計画段階で適切な設計が行える．
本稿では，乱流熱・物質輸送現象に対する予測手法で用いる「乱流熱伝達モデル（Turbulent Heat
Transfer Model）」について，3 回に分けて解説する．そこでは，乱流熱伝達モデルに関する基礎的な知識を始めとして，乱流熱伝達モデルの開発秘話？等も出来る限り記述していくが，拙い表
現がある場合にはご容赦いただければ幸いである．
２. 基礎方程式
2・1 支配方程式
デカルト座標系において，流体とエネルギー輸送の支配方程式は，次式で与えらえる（例えば，Hinze, 1975，エネルギー方程式については，田川, 2004 に解説がある）．
𝜕𝑢𝑖
𝜕𝑢𝑖
𝜕𝑝
𝜕
𝜕𝑢𝑖 𝜕𝑢𝑗 2 𝜕𝑢𝑘
𝜌(
+ 𝑢𝑗
)=−
+
[𝜇 (
+
− 𝛿
)] + 𝑓𝑖𝜕𝜏𝜕𝑥𝑗𝜕𝑥𝑖 𝜕𝑥𝑗𝜕𝑥𝑗 𝜕𝑥𝑖 3 𝑖𝑗 𝜕𝑥𝑘𝜕𝜌 𝜕𝜌𝑢𝑖
+
=0
𝜕𝑡
𝜕𝑥𝑖
𝜕𝑡
𝜕𝑡
𝐷𝑝
𝜕
𝜕𝑡
𝜕𝑢𝑖 𝜕𝑢𝑗 𝜕𝑢𝑖 2 𝜕𝑢𝑘 2
𝜌𝑐𝑝 ( + 𝑢𝑗
)=𝑞+
+
(𝜆
) + 𝜇 [(
+
)
− (
) ]
𝜕𝜏
𝜕𝑥𝑗
𝐷𝜏 𝜕𝑥𝑗 𝜕𝑥𝑗𝜕𝑥𝑗 𝜕𝑥𝑖 𝜕𝑥𝑗 3 𝜕𝑥𝑘

(1)
(2)
(3)

ここで，式(1)は，Navier-Stokes（運動）方程式，式(2)は連続の式，式(3)はエネルギー方程式（Equaitonof Energy）であり，𝐷⁄𝐷𝜏 = 𝜕⁄𝜕𝜏 + 𝑢𝑗 (𝜕⁄𝜕𝑥𝑗 )は実質微分である．式中の記号は，𝑢𝑖 は𝑥𝑖 方向速度，
𝑥𝑖 は𝑖方向，𝑝は圧力，𝜌は密度，𝜇は粘性係数，𝑓𝑖 は外力，𝜏は時間，𝑡は温度，𝑐𝑝 は定圧比熱，𝜆は熱伝導率，𝑞は外部から与えられる熱量である．これらの式のように，本稿では式の記述にテンソルを用いるため，以下にテンソルで用いられる規約や式等を示す．まず，式(1)含まれる𝛿𝑖𝑗 はクロネッカーのデルタ（Kronecker delta）であり，以下で定義される．
𝛿𝑖𝑗 = 1: 𝑖 = 𝑗}
𝛿𝑖𝑗 = 0: 𝑖 ≠ 𝑗(4)

例えば，このクロネッカーのデルタを用いると，微分の添字を以下のように変えることが出来る．
𝛿𝑖𝑗

𝜕
𝜕
=
𝜕𝑥𝑗 𝜕𝑥𝑖(5)

このため，式（1）の右辺第 1 項は，𝜕𝑝⁄𝜕𝑥𝑖 = 𝛿𝑖𝑗 (𝜕𝑝⁄𝜕𝑥𝑗 )となり，右辺第 2 項の中に，圧力を含めた応力としてまとめることが出来る．また，テンソルにおける重要な規約として，同じ項に同じ添字がある場合には，以下のように添字を１～３（空間の方向数）まで加算する総和規約がある．
(6)
𝛿𝑖𝑖 = 𝛿11 + 𝛿22 + 𝛿33 = 1 + 1 + 1 = 3
𝜕
𝜕
𝜕
𝜕
(7)
𝛿𝑖𝑗
= 𝛿𝑖1
+ 𝛿𝑖2
+ 𝛿𝑖3
𝜕𝑥𝑗
𝜕𝑥1
𝜕𝑥2
𝜕𝑥3
他に知っておくと便利なテンソル式として，次式で定義されるエディントンのイプシロン（Levi-Civita symbol）がある．

- 16 -

JSME TED Newsletter, No.86, 2018
𝜖𝑖𝑗𝑘 = 1: 𝑖𝑗𝑘 = 123, 231, 312
𝜖𝑖𝑗𝑘 = −1: 𝑖𝑗𝑘 = 132, 321, 213 }𝜖𝑖𝑗𝑘 = 0: 𝑖 = 𝑗, 𝑖 = 𝑘, 𝑗 = 𝑘(8)

本稿では作動流体の空間における密度変化を一定（𝜌 =一定）および物性値も一定と考え，つまり非圧縮性流体と仮定して扱うため，支配方程式は以下のように簡単になる．
𝜕𝑢𝑖
𝜕𝑢𝑖
1 𝜕𝑝
𝜕 2 𝑢𝑖
+ 𝑢𝑗
=−
+𝜈
+ 𝑓𝑖
𝜕𝜏
𝜕𝑥𝑗
𝜌 𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗
𝜕𝑢𝑖
=0
𝜕𝑥𝑖
𝜕𝑡
𝜕𝑡
𝜕2𝑡
+ 𝑢𝑗
=𝛼
+𝑞
𝜕𝜏
𝜕𝑥𝑗
𝜕𝑥𝑗 𝜕𝑥𝑗

(9)
(10)
(11)

ここで，𝜈 (= 𝜇⁄𝜌)は動粘性係数，𝛼 (= 𝜆⁄𝑐𝑝 𝜌)は温度伝導率である．
DNS では，これらの微分方程式を連立させて解いているが，ここで連立方程式として解けるかの確認をしておく．式の数を数えると，式（9）は𝑢1 , 𝑢2 , 𝑢3 の 3 方向速度について 3 本の式となり，式（10）の連続の式，式（11）のエネルギー方程式と合わせて 5 本の式がある．一方，変数の数は，𝑢1 , 𝑢2 , 𝑢3 の速度成分と，圧力𝑝，温度𝑡の 5 つである．つまり，式の数＝変数の数となり，方程式は閉じているため解くことが出来る（𝑓𝑖 や𝑞については，速度や温度の関数であるとして，変数の数は増えないとしている）．
2・2 レイノルズ平均
さて，乱流現象は時々刻々と変化するため，瞬時の情報だけを見てもその全体の情報を把握することは困難である．そのため，現象を時間平均で考えるため，前章で述べた次式で与えられる
レイノルズ平均（Reynolds Average，もしくはレイノルズ分解（Reynolds Decomposition））が考えられた（図１参照）．
̂𝑖 (𝑡, 𝑥𝑖 ) = 𝑈
̅𝑖 (𝑥𝑖 ) + 𝑢𝑖 (𝑡, 𝑥𝑖 )𝑈

(12)

ここで，左辺は瞬時速度（Instaneous Velocity），右辺第 1 項は時間平均値（Time-Mean Velocity），
̅𝑖 + 𝑢𝑖 ），圧力
第 2 項は乱れ成分（変動値，Velocity Fluctuation）である．この関係を速度（𝑢𝑖 = 𝑈̅̅（𝑝 = 𝑃 + 𝑝）
，温度（𝑡 = 𝑇 + 𝑡）に適用し，式（9）～（11）に代入して，時間平均を取ると次式
̂𝑖 (𝑡, 𝑥𝑖 ) = 𝑢𝑖を得る（なお，式（9）～（11）では，瞬時速度の記号を小文字として標記しているため，𝑈であるが，右辺の乱れ成分とは区別する）．

𝑢𝑖
̂𝑖 (𝑡 )
𝑈

̅𝑖
𝑈

𝑡
Fig. 1 Reynolds Decomposition- 17 -

JSME TED Newsletter, No.86, 2018
̅𝑖
̅𝑖
𝜕𝑢
̅̅̅̅̅
𝜕𝑈
1 𝜕𝑃̅
𝜕2𝑈
𝑖 𝑢𝑗
=−
+𝜈
−
+ 𝐹̅𝑖
𝜕𝑥𝑗
𝜌 𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗
𝜕𝑥𝑗
̅𝑖
𝜕𝑈
=0
𝜕𝑥𝑖
̅̅̅̅
𝜕𝑇̅
𝜕 2 𝑇̅
𝜕𝑢
𝑖𝑡
̅𝑗
𝑈
=𝛼
−
+ 𝑄̅
𝜕𝑥𝑗
𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗
̅𝑗
𝑈

(13)
(14)
(15)

上式を導出するにあたり，平均には次式の関係があることに注意されたい．
̅ +𝜙 =Φ
̅ + 𝜙̅ = Φ
̅ + 𝜙̅ = Φ̅𝜙̅ = ̅̅̅̅̅̅̅̅
Φ̅̅̅̅̅̅̅̅
̅ 𝐴Φ
̅𝐵 = Φ
̅𝐴 + Φ
̅𝐵 = Φ
̅𝐴 + Φ
̅𝐵
}
Φ
̅̅̅̅̅
̅
̅
̅̅̅̅̅̅̅̅̅̅
Φ𝜙 = Φ𝜙 = Φ𝜙 = 0, 𝜙𝐴 𝜙𝐵 ≠ 0(16)式(13)と式(15)を見ると，右辺に新たな項が出現している．これらは，左辺の慣性項からレイ
ノルズ平均をすることによって現れた乱流の本質を支配する項であり，それぞれレイノルズ応力，
乱流熱流束である．レイノルズ応力は，その成分を行列で示すと以下のようになる．
̅̅̅𝑢12
𝑢𝑖 𝑢𝑗 = (̅̅̅̅̅̅
̅̅̅̅̅
𝑢2 𝑢1
𝑢3 𝑢1
̅̅̅̅̅̅

̅̅̅̅̅̅
𝑢
1 𝑢2
̅̅̅
𝑢2

𝑢1 𝑢3
̅̅̅̅̅̅
(17)̅̅̅̅̅̅
𝑢2 𝑢3 )
̅̅̅𝑢3 𝑢2̅̅̅̅̅̅
𝑢32
2行列の対角成分は，レイノルズ垂直応力（Reynolds Normal Stress）と呼ばれ，また非対角成分はレイノルズせん断応力（Reynolds Shear Stress）と呼ばれるが，対称性（例えば，𝑢̅̅̅̅̅̅
𝑢2 𝑢1
1 𝑢2 = ̅̅̅̅̅̅）̅̅̅̅によりレイノルズ応力の全成分は 6 成分となる．また，乱流熱流束𝑢𝑡は，𝑥~𝑥のそれぞれの方𝑖
1
3
向の 3 成分である．
これらの項は定数ではなく，場の未知数であるため，方程式の数が増えてないにもかかわらず，変数が合計 9 つ増えたことになる．そのため，方程式の数≠変数の数となり，連立方程式が閉じていないことになる．つまり，なんらかの方法でレイノルズ応力，乱流熱流束を与えなければ連立方程式は解けないこととなる．まず，直観的にレイノルズ応力の輸送方程式を導出し，解けばよいのではと？考える．レイノルズ応力の輸送方程式は以下の手順で作成できる．
最初に，Navier-Stokes 方程式より，変動速度成分𝑢𝑖 と𝑢𝑗 の輸送方程式を導出する．
̅𝑖
𝜕𝑢𝑖
𝜕𝑈
𝜕𝑢
𝜕
1 𝜕𝑝
𝜕 2 𝑢𝑖
̅𝑘 𝑖 +(𝑢𝑖 𝑢𝑘 − 𝑢𝑖 𝑢𝑘 ) =
+ 𝑢𝑘
+𝑈
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘
̅𝑗
𝜕𝑢𝑗
𝜕𝑈
𝜕𝑢𝑗
𝜕 2 𝑢𝑗
𝜕
1 𝜕𝑝
̅𝑘
+ 𝑢𝑘
+𝑈
+
(𝑢𝑗 𝑢𝑘 − 𝑢𝑗 𝑢𝑘 ) =
+𝜈
𝜕𝑡
𝜕𝑥𝑘
𝜕𝑥𝑘 𝜕𝑥𝑘𝜌 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘(18)(19)(18) × 𝑢𝑗 + (19) × 𝑢𝑖 とすることで，次のレイノルズ応力𝑢𝑖 𝑢𝑗 の輸送方程式を得る．
̅𝑘
𝑈

̅𝑗
̅𝑖
𝜕𝑢𝑖 𝑢𝑗𝜕𝑈𝜕𝑢𝑖 𝑢𝑗 𝑢𝑘 𝑝 𝜕𝑢𝑖 𝜕𝑢𝑗𝜕𝑈= − (𝑢𝑖 𝑢𝑘+ 𝑢𝑗 𝑢𝑘
)−
+ (
+
)
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜕𝑥𝑘
𝜌 𝜕𝑥𝑗 𝜕𝑥𝑖

𝜕 2 𝑢𝑖 𝑢𝑗
1 𝜕𝑝𝑢𝑖 𝜕𝑝𝑢𝑗𝜕𝑢𝑖 𝜕𝑢𝑗
(20)
− (
+
)+𝜈
− 2𝜈
𝜌 𝜕𝑥𝑗
𝜕𝑥𝑖
𝜕𝑥𝑗 𝜕𝑥𝑗𝜕𝑥𝑘 𝜕𝑥𝑘
ここで，注目すべきは，右辺第２項に 3 次モーメント（Third-order Moment）である𝑢𝑖 𝑢𝑗 𝑢𝑘 を含む項が出現していることである．つまり，2 次モーメント（Second-order Moment）であるレイノ
ルズ応力𝑢𝑖 𝑢𝑗 の輸送方程式を作成しても，その 1 次高いモーメントが出現し，これも未知数であるため，この方法では永遠に方程式が閉じることがない．そのため，人為的，強制的に方程式を閉じさせる必要がある．このため，この問題のことをクロージャー問題（Closure Problem）という．そのクロージャー問題で行うことは，レイノルズ応力や乱流熱流束をなんらかの形で与える- 18 -

JSME TED Newsletter, No.86, 2018
ことに他ならない．乱流の本質であるレイノルズ応力や乱流熱流束を，なんらかの形で与えるため，このことを「乱流のモデリング（Turbulence Modeling）」という．
３. 乱流のモデリング
3・1 勾配拡散型モデル
レイノルズ平均を，乱流熱伝達現象の支配方程式に施すことにより出現したレイノルズ応力と乱流熱流束のモデリングに対し，一番基礎的なモデリングは勾配拡散型モデル（Gradient Diffusion
Type Model）と呼ばれる乱流モデルである．これは，Boussinesq が最初に提唱した渦粘性（Eddy
Viscosity）の概念(Boussinesq, 1877 を用い，次式で与えらえる．
̅𝑗
̅𝑖 𝜕𝑈
𝜕𝑈
1
−𝑢𝑖 𝑢𝑗 = 𝜈𝑡 (+) − 𝑢𝑘 𝑢𝑘 𝛿𝑖𝑗𝜕𝑥𝑗 𝜕𝑥𝑖3̅𝜕𝑇−𝑢𝑖 𝑡 = 𝛼𝑡
𝜕𝑥𝑖

(21)
(22)

ここで𝜈𝑡 は，渦動粘性係数（Eddy Viscosity）もしくは運動量の渦拡散係数（Eddy Diffusivity for
Momentum），𝛼𝑡 は熱の渦拡散係数（Eddy Diffusivity for Heat）である．また，式(21)はテンソル表
記のため，添字が𝑖 = 𝑗となった時に右辺＝左辺とするために右辺に項が付加されている．
このモデルの利点は，これらの関係式を支配方程式(13)と(15)に代入すると理解できる．
̅𝑖
̅𝑖
𝜕𝑈
1 𝜕𝑃̅
𝜕
𝜕𝑈
=−
+
[(𝜈 + 𝜈𝑡 )] + 𝐹̅𝑖𝜕𝑥𝑗
𝜌 𝜕𝑥𝑖 𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑇̅
𝜕
𝜕𝑇̅
̅𝑗
𝑈
=
[(𝛼 + 𝛼𝑡 )] + 𝑄̅
𝜕𝑥𝑗 𝜕𝑥𝑗
𝜕𝑥𝑗

̅𝑗
𝑈

(23)
(24)

これらの式を見ると，粘性係数（拡散係数）と渦拡散係数が，実効拡散係数（𝜈𝑒 = 𝜈 + 𝜈𝑡 ，𝛼𝑒 =𝛼 + 𝛼𝑡 ）としてまとめられることが分かる．層流の場合には𝜈𝑡 = 0，𝛼𝑡 = 0となり乱流の作用が消え，乱流の場合にはこれらが値を持つことから実効拡散係数が増加し，運動量や熱の拡散が活発となることが説明できる（第 1 章の乱流の定義で記述した「大きな混合と拡散」は，これで説明が出来る）．
また，このモデルの欠点としては，レイノルズ応力の垂直応力成分（乱れの各成分）が正確に算出できないことにある．例えば，平行平板間で十分に発達したスパン方向に一様な乱流を考え̅ (=る．流れ方向を𝑥 (= 𝑥1 )，壁面垂直方向を𝑦 (= 𝑥2 )，スパン方向を𝑧 (= 𝑥3 )流れ方向平均速度を𝑈
̅
̅
̅
̅
̅
̅
̅
𝑈1 )，壁面垂直方向速度を𝑉 (= 𝑈2 )，スパン方向速度を𝑊 (= 𝑈3 )とすると，条件により𝑉 = 𝑊 = 0，𝜕⁄𝜕𝑥 = 0，𝜕⁄𝜕𝑧 = 0となる．そのため，式(21)より，レイノルズ応力の垂直成分はそれそれ次式で与えらえる．
̅ 1
𝜕𝑈1+ 𝑞= 𝑞
𝜕𝑥 3
3
𝜕𝑉̅ 1
1
2
𝑣 = −2𝜈𝑡+ 𝑞= 𝑞
𝜕𝑦 3
3
̅ 1
𝜕𝑊
1
𝑤 2 = −2𝜈𝑡
+ 𝑞= 𝑞
𝜕𝑧 3
3

(25)

𝑢2 = −2𝜈𝑡

(26)
(27)

ここで，𝑞 = 𝑢𝑘 𝑢𝑘 である．上式から， 𝑢2 = 𝑣 2 = 𝑤 2 (= 𝑞⁄3)となることが分かるが，発達した平衡平板間乱流における実際の各方向乱れ分布は．実験や DNS の結果（図 2 参照（Nagano and Hattori,
2003））からも明らかのように，𝑢2 > 𝑤 2 > 𝑣 2と非等方（Anisotoropy）的な分布である．このように，関連する平均勾配が存在しない場合には，レイノルズ応力を等方（Isotoropy）的に与えてしまう．

- 19 -

JSME TED Newsletter, No.86, 20188𝑅𝑒𝜏 =150

𝑢2
𝑣2
𝑤2

4

0
-1

0

𝑦 ⁄𝛿

1

Fig. 2 Distributions of Reynolds normal stresses in filly-developed channel flow

さて，これら渦拡散係数は，作動流体の物質で決まる物性値ではなく，場の状態によって決まる変数である．そのため，勾配拡散型モデルでレイノルズ応力，乱流熱流束をモデル化しても，方程式は閉じていることにならない．よって，さらに渦拡散係数をモデル化する必要がある．
3・2 渦拡散係数を用いた勾配拡散型モデルの分類
渦拡散係数𝜈𝑡 ，𝛼𝑡 のモデル化は，分子粘性係数𝜈 [m2 ⁄s] が，分子の平均速度[m⁄s]×平均自由行
程[m]で算出されることを参考に，乱れの速度𝑉×乱れの空間スケール𝐿でモデル化する．
(28)

𝜈𝑡 = 𝐶𝑉 ∙ 𝐿

ここで，𝐶はモデル定数である．また熱の渦拡散係数𝛼𝑡 についても，その次元が同じのため同様のモデル化を行うが，温度場の場合には乱れの空間スケールの概念が定義しづらいため，乱れの寿命τ (= 𝐿⁄𝑉 )を用い，以下のようにモデル化する．
(29)

𝛼𝑡 = 𝐶𝑉 ∙ 𝑉τ

勾配拡散型モデルの分類は，この乱れの速度𝑉と乱れの空間スケール𝐿（もしくは乱れの寿命τ）を解くための方程式の数で分類され，速度場の運動量の渦拡散係数に対する代表的なモデルとして以下のモデルがある．
(1) 0 方程式モデル：プラントルの混合距離理論（Prandtl’s Mixing Length Theory (Prandtl, 1925)）
乱れの空間スケール𝐿を，実験的な知見から得られる混合距離ℓで与え，乱れの速度𝑉は，̅ ⁄𝜕𝑦)で与える（𝑈̅は主流方向平均速度，𝑦は壁面垂直方向
平均速度勾配と混合距離の積ℓ(𝜕𝑈
座標）．そのため，渦拡散係数を求めるために新たな方程式を解く必要がなく，0 方程式モデルと呼ばれる．
̅ ⁄𝜕𝑦) ∙ ℓ = 𝐶ℓ2 (𝜕𝑈̅ ⁄𝜕𝑦)
𝜈𝑡 = 𝐶𝑉 ∙ 𝐿 = 𝐶ℓ(𝜕𝑈(2) 1 方程式モデル：乱流エネルギー輸送方程式モデル，渦拡散係数輸送方程式モデル(Boldwin
and Barth, 1990)
乱れの速度𝑉を，𝑘 = ̅̅̅̅̅
𝑢𝑖 𝑢𝑖 ⁄2で定義される乱流エネルギー（Turbulent Kinematic Energy）の- 20 -

JSME TED Newsletter, No.86, 2018
輸送方程式で解き（速度スケールは𝑉 = √𝑘），乱れの長さスケールは混合距離ℓで与える．
また，渦拡散係数そのものの輸送方程式を解くモデルも 1 方程式モデルに分類される．
𝜈𝑡 = 𝐶𝑉 ∙ 𝐿 = 𝐶√𝑘 ∙ ℓ(3) 2 方程式モデル：𝑘 − 𝜀モデル(Launder and Spalding, 1974)，𝑘 − 𝜔モデル(Wilcox, 2008)𝑘 − 𝜀モデルは，乱れの速度𝑉に乱流エネルギーで定められる√𝑘を，乱れの空間スケール𝐿に
乱流エネルギーの散逸率𝜀（Dissipation rate of Turbulence Energy）を用い，𝐿 = 𝑘 3⁄2⁄𝜀で与える．乱流エネルギーとその散逸率は，それぞれの輸送方程式から算出される．
𝜈𝑡 = 𝐶𝑉 ∙ 𝐿 = 𝐶√𝑘 ∙ (𝑘 3⁄2⁄𝜀) = 𝐶 𝑘 2⁄𝜀𝑘 − 𝜔モデルは，乱れの速度𝑉は𝑘 − 𝜀モデルと同様であるが，乱れの空間スケール𝐿ではな
く，乱れの寿命τ (= 𝐿⁄𝑉 )を用い，乱れの寿命τを乱れの特性周波数𝜔でτ = 1⁄𝜔としたモデルである（𝑘 − 𝜀モデルとは，𝜔 = 𝜀 ⁄𝛽 ∗ 𝑘の関係がある（𝛽 ∗はモデル定数））．
𝜈𝑡 = 𝐶𝑉 ∙ 𝐿 = 𝐶𝑉 ∙ 𝑉τ = 𝐶√𝑘 ∙ (√𝑘⁄𝜔) = 𝐶 𝑘 ⁄𝜔一方，温度場の熱の拡散係数に対する代表的なモデルとして以下のモデルがある．
(1) 0 方程式モデル：プラントル数（粘性係数と熱伝導率の比，𝑃𝑟 = (𝜇⁄𝜌)⁄(𝜆⁄𝜌𝑐𝑝 ) = 𝜈⁄𝛼））の考えを，運動量の渦拡散係数と熱の渦拡散係数の比の関係に拡張した乱流プラントル数（𝑃𝑟𝑡 = 𝜈𝑡 ⁄𝛼𝑡 ）によって与えられるモデル．通常，乱流プラントル数は定数で与えられる
ため，温度場 0 方程式モデルと呼ばれる．
𝛼𝑡 = 𝜈𝑡 ⁄𝑃𝑟𝑡
(2) 1 方程式モデル：熱の渦拡散係数輸送方程式モデル(裴他, 1990)
熱の渦拡散係数そのものを輸送方程式を解くモデル．
(3) 2 方程式モデル：温度場 2 方程式モデル(Nagano and Kim, 1988)
乱れの速度𝑉に乱流エネルギーで定められる√𝑘を，乱れの寿命τに温度乱れ𝑡̅2とその散逸率𝜀𝑡 から定まる温度場特性時間スケール𝜏𝑡 (= 𝑡̅2⁄𝜀𝑡 )を導入し，速度場の特性時間スケール

𝜏𝑢 (= 𝑘⁄𝜀 )とのタイムスケール比𝑅 (= 𝜏𝑡 ⁄𝜏𝑢 )を用いてτ = 𝜏𝑢 𝑓(𝑅)として与える（Nagano andKim, 1988 の場合は，𝑓(𝑅) = √2𝑅で与えらえる）．温度乱れとその散逸率は，それぞれの輸送方程式から算出されるが，別途速度場タイムスケールを出すための方程式が必要である．
𝛼𝑡 = 𝐶𝑉 ∙ 𝑉τ = 𝐶√𝑘 ∙ √𝑘𝜏𝑢 𝑓(𝑅) = 𝐶(𝑘 2 ⁄𝜀)𝑓(𝑅)

以上により，レイノルズ応力と乱流熱流束に対する，勾配拡散型モデルを用いた基本的なモデル化の概念を説明した．次回は，これらの更なる具体的なモデリング方法について説明をする．
文献
Boussinesq, J., Essai sur la théorie des eaux courantes, Mémoires présentés par divers savants à l'Académiedes Sciences, Vol. 23, No. 1 (1877), pp.1-680.
Baldwin, B.S. and Barth, T.J., A one-equation turbulence transport model for high Reynolds number
wall-bounded flows, NASA TM 102847 (1990).
Durst, F., Melling, A. and Whitelaw, J. H., The application of optical anemometry to measurement in
combustion systems, Combustion and Flame, Vol.18 (1972), pp. 197-201.
Hinze, J. O., Turbulence (2nd Edition), McGraw-Hill (1975).
Kim, J., Moin, P., and Moser, R., Turbulence statistics in fully developed channel flow at low Reynolds
number, Journal of Fluid Mechanics, Vol. 177 (1987), pp. 133-166.
Morinishi, Y., Lund, T.S., Vasilyev, O.V. and Moin, P., Fully conservative higher order finite difference- 21 -

JSME TED Newsletter, No.86, 2018
schemes for incompressible flow, Journal of Computational Physics, Vol. 143, No. 1 (1998), pp. 90-124.
Launder, B. E. and Spalding, D. B., The numerical computation of turbulent flows, Computer Methods in
Applied Mechanics and Engineering, Vol. 3, No. 2 (1974), pp. 269-289.
Nagano, Y. and Hattori, H., Direct numerical simulation and modelling of spanwise rotating channel flow
with heat transfer, Journal of Turbulence, Institute of Physics, Vol. 4-010 (15pages) (2003).
Nagano, Y. and Kim, C., A two-equation model for heat transport in wall turbulent shear flows, ASME. J.
Heat Transfer, Vol. 110, No. 3 (1988), pp. 583-589.
Nakajima, T., Utsunomiya, M., Ikeda, Y. and Matsumoto, R., Simultaneous measurement of velocity andtemperature of water using LDV and fluorescence technique, Applications of Laser Techniques to Fluid
Mechanics, Springer-Verlag, (1990) pp. 34-53.
Prandtl, L., Z. angew, Math. Mech. Vo. 5, No. 1 (1925), pp. 136-139.
Reynolds, O., An experimental investigation of the circumstances which determine whether the motion ofwater shall be direct or sinuous, and of the law of resistance in parallel channels, PhilosophicalTransactions of the Royal Society of London, Vol. 174 (1883), pp. 935-982.
Reynolds, O., On the dynamical theory of incompressible viscous fluids and the determination of thecriterion, Philosophical Transactions of the Royal Society of London, A, Vol. 186 (1895), pp. 123-164.
Wilcox, D. C., Formulation of the 𝑘 − 𝜔 Turbulence Model Revisited, AIAA Journal, Vol. 46, No. 11(2008), pp. 2823–2838.
Willert, C.E. and Gharib, M., Digital particle image velocimetry applied to an evolving vortex ring, Proc. of
the Fifth Intern. Symp. on Application of Laser Techniques to Fluid Mechanics and Workshop on the Useof Computer in Flow Measurements, 18.1 (1990).
佐々木達治郎, 玉野光男，空気力学実験法，河出書房(1941).
田川正人，下地敏央，太田安彦，二線式熱電対による変動温度測定法の精度評価と改良，日本機械学会論文集 B 編，Vol. 64, No. 625 (1998), pp. 3077-3083.
田川正人，永谷真一，下地敏央，太田安彦，高温乱流場における速度と温度の高精度同時測定，
日本機械学会論文集 B 編，Vol. 65, No. 634 (1999), pp.2157-2164.
田川正人，エネルギー式を巡る，日本伝熱学会誌「伝熱」，Vol. 43, No. 178 (2004), pp. 26-31.
菱田幹雄，長野靖尚，田代真一郎，速度変動と温度変動の同時測定，日本機械学会論文集，Vol.43,
No.365 (1977), pp.225-232.
裴昌謙，服部博文，長野靖尚,新しい低レイノルズ数型１方程式モデルによる熱･流体解析（第２報：温度場１方程式モデルの構築），Thermal Science & Engineering，Vol. 7, No. 5 (1999)，pp. 31-39.

- 22 -