第 4 章 数値解法 82
4.4 乱流モデル
4.4.1 渦粘性モデル
(a) レイノルズ応力
単位体積あたりのx方向の運動量がρ(u±u0) = ρu(y ±l)である流体塊が，y 軸に直角な単位 面積あたりに体積(∓v0)だけy軸の負の方向(あるいは正の方向)に移動する．その結果，平均的に E[ρ(u±u0)·(∓v0)] =E[−ρuv0−ρu0v0]の運動量がy方向に輸送される．ここに，v0は流体塊のy方 向への輸送速度であり，Eは統計的平均を意味する．したがって，運動量の法則によりy軸に直角な 単位面積には，輸送された運動量に等しいτyxが運動量成分の方向(x方向)に作用するものとしてみ なしうる．すなわち，E[v0] = 0を考慮すれば，
τyx=−E[ρ(u+u0)v0] =−E[ρu0v0]
=−ρu0v0
このような流れの乱流運動に伴う応力はレイノルズにより導かれた．この応力は，レイノルズ応力 (Reynolds stress)Ri,j と呼ばれる(76)．このレイノルズ応力は，乱流現象の理論の基礎をなすもので ある．上の説明ではτyx成分のみ導いたが，他の成分についても同様にレイノルズ応力が存在する．
すでに，前節において，NS方程式(前節の全成分の運動量の式)からレイノルズ応力を導き出して おり，平均流運動方程式と全成分の運動量方程式を区別する唯一の量としてあらわれている．この平 均流運動方程式は，レイノルズ平均NS方程式と呼ばれ，レイノルズによって導き出された．このレ イノルズ平均NS方程式のことをRANS(Reynoldes averaged Navier-Storks equations)やレイノル ズ方程式と呼ぶ場合もある．
これら，上記の関係は非圧縮性乱流に対してレイノルズが導いた．圧縮性乱流に対しても適当な仮 定のもと前節で示したように方程式が求めることができる．その方程式を遷音速・超音速流れの計算 流体力学研究者は，時間平均NS方程式(Time averaged Navie-Stokes equations，TANS)と呼ぶこ とが多い(73)．また，Morkovinの仮説より，圧縮性乱流において，−ρ¯ug00iu00j =−ρu0iu0jと扱うことが 出来る．したがって，−ρ¯ug00iu00j を圧縮性流体での乱流運動による応力(レイノルズ応力)Rijとして取 り扱う．
本研究で乱流モデルとして採用したRANS(Reynoldes averaged Navier-Storks equations)では，レ イノルズ応力Ri,jの取り扱いに関して2つの方法がある．第1の方法では，Ri,jを渦粘性の概念(渦 粘性モデル)を用いて近似する方法．第2の方法では，Rij の輸送方程式(レイノルズ応力方程式)を 直接取り扱う方法．本研究では，計算が容易な第１の渦粘性モデルによりレイノルズ応力を導出する. レイノルズ応力を具体的に表すには関係式が必要となる．最も単純な考え方はレイノルズ応力を平 均速度分布と関連付ける方法である．渦粘性の概念の下に平均速度および乱流を特徴付けるスカラー 量を用いてレイノルズ応力を近似する方法を，渦粘性(型)モデルと総称する(68)．
ここでは，渦粘性モデリングとして最も簡単なブシネスク(Boussinesq)近似と，その基となったプ ラントル(Prandtl)の混合長理論(混合距離理論)について説明する．後述の乱流モデルの
Baldwin-Lomaxモデルでは混合長理論を適応し，付録の中で述べる修正k−Rモデルにブシネスク近似を適
応している．
(b) ブシネスク近似
まず，考えの基本となった非圧縮性流体に対しての理論を述べる．その後，圧縮性流体の対応を述 べる．
・非圧縮性乱流
乱流渦によって発生する応力を分子粘性応力との類推から与える考え方がある．ここで分子粘性応 力は粘性係数と速度勾配の積で表現される事が実験的に確認されており，次式のように表される．
τij =ρν(∂uj
∂xi + ∂ui
∂xj)
分子速度の平均をv¯，平均自由行程をlF としたとき，分子動粘性係数νは次式と表される．
ν = 1 3vl¯F
乱流応力も同様に，渦による拡散を表す係数νtと平均速度勾配の積であると考えれば，
−ρu0iu0j =ρνt(∂uj
∂xi + ∂ui
∂xj)−2
3ρkδij (4.87)
これをレイノルズ応力のブシネスクの渦粘性近似と呼ぶ．νtは分子動粘性係数νに対応する渦運動効 果を表すもので，渦動粘性係数(kinetimatic eddy viscosity)または乱流動粘性係数とよばれる．νと は異なりνtは流体の性質では定まらず，流体運動の場所の関数となる．分子動粘性係数と同じく[長 さ]×[速さ]の次元をもち，次のように表すことができる．
νt=utl
ここでutは乱流渦の特徴的な速度であり，lはその特徴的な長さスケールである．ただし，分子動粘 性係数にみられる係数1/3に対応する定数はl の定義に含まれるものとして省略されている．
また，式(4.87)の右辺第二項2/3·ρkδijは両辺の縮約をとったときの恒等関係から必要となる項で あり，レイノルズ応力の等方成分を表している．
式(4.87)が高い精度で成立するためには平均場と変動場におけるスケールの分離と変動場の高い等
方性が必要となり，工学分野で実際に取り扱われる複雑な乱流場では十分に満足されない事が多い．
・圧縮性乱流
上記の関係を圧縮性流体に適応すると，せん断応力に体積膨張の効果が加わり分子粘性応力は次式 で示される．
τij =µl(∂ui
∂xj +∂uj
∂xi) +λδij∂uk
∂xk ここで第二粘性係数λを−2/3·µlとする．
τij =µl³∂ui
∂xj + ∂uj
∂xi −2 3δij∂uk
∂xk
´
=µlSij 乱流応力(レイノルズ応力)では，
Rij =−ρ¯ug00iu00j =−ρu¯ 0iu0j =µt³∂uei
∂xj
+∂fuj
∂xi − 2 3δij∂ufk
∂xk
´− 2
3ρkδ¯ ij =µtSij −2 3ρkδ¯ ij
このモデルを導入することで，前節のRANS方程式の粘性項では，次式のようになる．ここで分子 粘性応力の速度は，前章で導出したようにファーブル平均である．
τij +Rij = (µl+µt)Sij− 2 3ρkδ¯ ij
ここで Sij = ∂uei
∂xj +∂ufj
∂xi −2 3δij∂ufk
∂xk と簡単な形となり，計算が容易である．
次に問題となるのは，レイノルズ応力内の渦粘性係数µtをどう評価するかである．この評価方法 に3つの方法がある．µtを求める際に必要な乱れの特性を的確に反映する物理量の数つまり，追加さ れる方程式の数で呼び方が決められる．
• 0方程式モデル …… 新たな方程式を解かなくて良い．
• 1方程式モデル …… 1つの方程式を解いて求める．
• 2方程式モデル …… 2つの方程式を解いて求める．
(c) 混合長理論
ここでは，レイノルズ応力の近似法だけでなく，渦粘性係数の導出方法まで述べる．ただし，流れ 場は，考えが容易な非圧縮乱流で，かつ単純せん断流(クエット流)として述べる．
ある点に生じる流速uの変動分u0 は流体粒子の移動距離l1とその平均時間勾配にd¯u/dyに比例す る．v0も同じように平均速度勾配に比例すると考えられる．
u0 =l1d¯u
dy, v0 =l2d¯u dy
流速の変動分を時間(もしくは空間か集合)の二乗平均を行い，レイノルズ応力を求める．ここでα はu0v0 の相関係数で，さらにl2 =αu0v0 とする．
−ρu0v0≈ρu0rmsvrms0 ≈ρα(l1d¯u dy)(l2d¯u
dy) =ραl1l2(d¯u dy)(d¯u
dy) =ρl2|d¯u dy|d¯u
dy
上式の絶対値記号は，レイノルズ応力−ρu0v0と∂u/∂y¯ が同符号となるように考慮したためである．
ここで，前節で行ったように，乱流渦によって発生する応力(レイノルズ応力)を分子粘性応力と の類推から与える．ここで分子粘性応力は，次式(ニュートンの粘性抵抗則)でありこの式と対応さ せる．
τxy =ρνdu dy
したがって，レイノルズ応力Rij と乱流動粘度νtは次式となる．
Rij =−ρu0v0 =ρνtd¯u dy νt=l2|d¯u
dy|
この理論は，運動量が混合距離の間で保存されて運搬されると考えるので運動量輸送理論，またはプ ラントルの混合距離理論(混合長理論)(mixing length theory)と呼ばれ，lは混合距離または混合長 (mixing length)と呼ばれている．
通常はνtよりも，µt=ρνtが用いられ，µtを乱流動粘度(eddy kinematic viscosity)，もしくは乱 流拡散係数(eddy diffusion coefficient)と呼ぶ．
Rij =−ρu0v0 =µtd¯u dy
この混合長理論を，本研究で対象とする一般的な圧縮性乱流に対応させる．分子粘性応力とレイノル ズ応力は，次式で示される．
τij =µl³∂ui
∂xj
+ ∂uj
∂xi −2 3δij∂uk
∂xk
´
=µlSij
Rij =−ρ¯ug00v00 =−ρu¯ 0v0 =µt³∂uei
∂xj +∂fuj
∂xi −2 3δij∂ufk
∂xk
´
=µtSij
このモデルを導入することで，前節のRANS方程式の粘性項では，次式のようになる．ここで分子 粘性応力の速度は，前章で導出したようにファーブル平均である．
τij +Rij = (µl+µt)Sij
ここで，
Sij = ∂uei
∂xj +∂fuj
∂xi −2 3δij∂ufk
∂xk
このように，簡単な形となり，前節のブシネスク近似よりも計算が容易である．
また，渦粘性係数は次式で示される．速度勾配成分にひずみ速度テンソルを用いた場合，
µt=ρ²=ρl2|∂uei
∂xj
+∂fuj
∂xi|=ρl2
s³∂uei
∂xj
+∂fuj
∂xi
´2
速度勾配成分に渦度を用いた場合．本研究では2次元より，
µt=ρ²=ρl2|∂u˜
∂y − ∂˜v
∂x|=ρl2 s³∂u˜
∂y −∂x˜v
∂x
´2
ここで混合長lは，壁面からの距離ynに比例すると考えられ，
l=κyn または l=κynf
である．ここでKarman定数κ(= 0.41)，f はVan-Driestの減衰関数で，
f = 1−exp−y+ A
である．ただし，定数A(= 26)である．減衰関数f は境界層の粘性底層内で乱れが減衰していく状 態をしめす関数で壁面上で0，外層で１となる関数である．y+は，後述する壁法則で壁面近くの流速 分布決定する無次元量である．
以上のように，この理論を使うことによって，新たな方程式を解かず流れ場のみの情報で渦粘性係 数νtを求め，レイノルズ応力Rij を決める事ができる
本章では，0方程式モデルの代表といえる混合長理論を応用したBaldwin-Lomaxモデルについて 説明する．また，ブシネスク近似を使った2方程式モデルの修正k−Rモデルについて付録で説明 する．