第 10 章 乱流と乱流解析
この章では乱流の基本的性質と乱流のシミュレーションについて述べる．乱流の計算は，最近 コンピュー
タの性能が飛躍的に向上したので，直接数値シミュレーションによってもある範囲のものを扱うことができる．しかし通常は計算量を減らすために乱流モデルを用いて行われる．乱流モデルは乱流を記述する数式のことである．各種のモデルが提案されているが，現状はまだ完成の域に達していない．そのために乱流を計算する場合には，適切なモデルを選択することとその限界を知ることが必要である．ここでは乱流の基本的性質について簡単に触れた後に，乱流のモデルと計算法について概観する．

10.1

乱流の性質

乱流は，大小さまざ まの渦塊 (eddy，単に渦ともいう) を含み非常に複雑な流れである．それでも乱流は，基本的に Navier-Stokes 方程式によって支配される．非圧縮性流れの Navier-Stokes 方程式は次のように書くことができる．

duu @uu +u ru = ; 1 rp +  r2u(10.1)dt @t

ただし u は流速，p は静圧， は動粘性係数 (kinematic viscosity) である．この方程式の意味を復習すれば，左辺は対流項 (convection term)，右辺第 1 項は圧力の項，第 2 項は分子粘性による拡散項 (viscous di usionterm) である．また d=dt @=@t+ur は実質微分 (substantial derivative) の演算子で，時空間 x t における流跡線に沿う内微分 (interior derivative) である．式 (10.1) を１つの流跡線に沿って積分すれば，Bernoulliの定理に相当のものが得られるが，流跡線に沿って時間間隔 t の間 積分すれば次式が得られる．

u(x + xx t + t) = u(x  t) +x=
ただし xZ t+t;t; 1 rp +  r2u t(10.2)

R t+t ut，またこれらの式の積分路は流跡線上に取られる．式 (10.2) をことばで言えば，流tx だけ移動するが，その流速はこの間に圧体素片 ( uid particle) は時間 t の間に１つの流跡線に沿って x

R t+tf;(1=)rp +  r2ugt だけ変化するということになる．なお拡散

力勾配と粘性拡散の作用をうけ， t
項の

r2u は u の平均曲率で，その１つの成分 r2 u が正ということはその点の u の値がそのまわりの u の平均値よりも小さいということである．このとき式 (10.2) により u の値は増加することになる． は増加の程度を示す拡散係数で正数である．拡散項はこのように流速の遅いところを加速し，速いところを減速し均一化する働きを持つが，漠然とした言い方になるが，全体の流速の平均値を変えるものではない1 ．流れは一般に高レ イノルズ数 (Reynolds number) で乱流になる．レ イノルズ数は Re = UL= のように定義さ1 拡散の意味は塩水の濃度拡散の場合には分かり易い．拡散によって濃度は次第に一様化するが，濃度の平均値は溶解している塩の質量が不変であるから一定に保たれる．

1

2

乱流と乱流解析れ，式 (10.1) の 対流項]/ 粘性項] の比を表している2 ．U は代表速度で例えば一様流中に置かれた円柱周りの流れでは一様流の流速，L は代表長さでこの流れでは円柱の直径である．一般に遅い流れ，微小な流れ，
粘性の小さくない高粘度の流れは乱流にはならない．また極超音速 (hypersonic) 流れも別の理由で乱流にはならない．

Navier-Stokes 方程式の左から rを演算し連続の条件 ru = 0 を用いれば，次の渦度輸送方程式 (vorticity
transport equation) が導かれる．

ddt@ +(u r) = ( r)u +  r2
@t(10.3)

ここではこの方程式の導き方は省略し，解釈についてだけ述べる．この式の左辺は渦度の対流項，右辺第 1
項は渦度の発生項 (generation term，生成項)，第 2 項は渦度の拡散項である．第 2 項が発生項であること

は，式 (10.3) で拡散項を無視した場合に，この式の解釈上そうならざ るを得ないということはあるが，渦度の発生がなぜ (r)u のように表せるのかについては多少説明を要するところである．

図 10.1: 渦度の輸送

x 空間内のある領域に定義される連続関数を (x ) とする．点 x 0 における関数値は (x 0 )，この点から空間ベクトル a だけ離れた点 x 0 + a における関数値は Taylor 展開により (x 0 + a ) = (x 0 )+(a r) (x 0 )+(1=2!)(a r)2 (x 0 )+    となる．式をたてる際には通常 1 次の微少量まで考えれば十分である．図 10.1 に示すように点 A の流速を u ，渦度を  とする．点 A から  だけ離れた点 B の流速は u +( r)u ，また点 Aから流跡線に沿って u だけ離れた点 C の渦度は  + @ =@t +( r) =  + d =dt となる．今式 (10.3) で拡散
項を無視すれば ，この式は d =dt = ( r)u となるから，図に見るようにベクトル u +( r)u とベクトル
 +d =dt の先端は一致し１つの点 D 上にくる．これより拡散のない流れでは渦度は流体に凍結され (frozenin)，流体と共に移流することが分かる．また初めに AB にあった流体が単位時間後に CD まで移動するとき，点 A と点 B の流速が ( r)u だけ異なれば ，AB にあった流体は移動に伴って伸縮し 同時に渦度ベク
トルも伸縮しその強さが増減することになる．渦度の強さは流れの中に分布する渦糸 (vortex laments) の密度で表すことができる3 ．初め AB にあった流体が移動に伴って伸長し断面積が縮小すれば，渦糸は流体に凍結しているので，その断面の単位面積当たりの渦糸の本数，すなわち渦度は増加することになる，というように説明することもできる．２次元の岐点流れでは，上流側で等方乱れ (isotropic turburence) であっても岐点に近付くにつれて異方性 (anisotropy) が出てくる．それは流線に平行な渦度ベクトルは岐点に向かうときに縮小し，垂直なベクトルは伸長し，紙面に垂直なベクトルは不変であるからである．乱れのどの成分が増幅しどの成分が不変なのかは，渦度ベクトルの伸縮を考えれば容易に分かることである．

2 次元解析すれば， 対流項]
U 2=L， 粘性項] U=L2 である．これらの比を取ればレ イノルズ数になる．また微分方程式の各項は同じ次元を持つべきであるから，レ イノルズ数は無次元量で，動粘性係数  は LU の次元を持つ．
3 この渦糸は非粘性流れではその循環が一定で (constancy of circulation)，また閉じているかその端点は境界上にある．

3

乱流と乱流解析
完全な２次元平面流れでは渦の発生項 (r)u は常にゼロになる．それは，この流れでは  r は紙面に垂直方向の微分で，流速 u はこの方向に不変であるからである．２次元平面流れでは，渦度方程式に発生項がないので，物体表面で生じた渦度は，表面を離れて強くなることはなく，拡散によって次第に弱まるだけである．このような状況下では流れは乱流にならないはずである．しかしながら実在の平板や翼面上の２次元
平面的な流れは乱流に遷移する．この乱流遷移は次のようなプロセスを経て実現することが知られている．
境界層流れ (boundary layerow) は始め単純な剪断流であるが，まず横渦の TS 波 (Tollmien-Schlichtingwave) が発生する．この段階までの流れは２次元的である．これとほぼ同時に固体表面に遅い流れの縦縞，
ストリーク構造が現れる．横渦は，縦縞の遅い部分にのし上げると，その部分が外側の速い流れよって先へ引きずられ，ヘアピン渦 (hairpin vortex) に変形する．ヘアピン渦の渦度は渦の引き伸ばしのため強く，その先端 (head) は境界層の外縁に達し，またその基部 (leg) には小規模ではあるが強力なバースト渦 (bursting

vortex) が現れる．バーストは壁面近くの遅い流体の吹き上げ (ejection) とそれに続く速い流体の壁への吹き付け (sweep) からなる現象で，乱流の乱れの大部分はこのようにして作られる．ストリーク構造，ヘアピ

ン渦，バーストは３次元的現象である．なおレ イノルズ数が低い場合にはヘアピン状ではなく馬蹄形の渦(horseshoe vorex) になる．このように２次元境界層内の流れは線形波の成長，非線形構造の出現，乱雑化の一連のプロセスを経て乱流に遷移するのである．
乱流ははじめ統計的に乱れた流れであると考えられ，気体論 (kinetic theory) で顕著な成果を出した統計的手法を基に理論的研究が進められた．分子には乱子 (eddy)，平均自由行程には混合長 (mixing length)，   が対応付けられ理論が展開されたが，次第に高度化するにつれて理論のための理論の観を呈するようになった．それは乱流が単なる統計的な乱れでなく，組織構造を持つことを見過したために，実在の多くの乱流を説明できなかったからである．乱流の組織構造 (turbulence structure) を研究したものには，古くはRosenhead(1931) による混合層の離散渦法による安定解析，Townsend(1949) による円柱後流の二重構造の可視化などがある．しかしながらその多くは最近の流れ可視化技術の進歩4 ，プローブで得られたデータの処理技術の開発5 ，更に決定的なのは直接数値シミュレーション (DNS, direct numerical simulation) によるものである．乱流の DNS によって得られた結果は乱流データベース (turbulent

ow database) といわれ，そこからは非定常３次元の流速分布はもとより，実験では得難い豊富な情報，圧力，渦度，ヘリシティ(helicity) 6の分布，更にはレ イノルズ応力成分，乱れのエネルギーとその散逸など 各種の統計量も引き出すことができる．乱流データベースの最初で有名なものは Stanford 大学のグループが Illiac IV コンピュータを使い，一様剪断乱流を 64  64  64 格子点を用い Navier-Stokes 方程式をスペクトル法で解いて作ったものである．このデータベースを基に 1982 年頃に各種の乱流モデルの検証が行われた．また NASA Ames 研
究所の Moin-Kim によって壁面剪断流のヘアピン渦のシミュレーションも行われた．その直後に Illiac IV
コンピュータは運用を停止し，当時の文献によればこのような研究は当分行えないだろうとのことであったが，今では多くのデータベースが作られ利用できるようになっている．
工学で重要な物体表面の境界層やその後流などの，いわゆる薄剪断層乱流 (turbulent shearow) は，乱流構造に等方的乱れが重なり合ったものとみることができる．また乱れエネルギーの大部分は小スケール乱interferogram(密度), Schlieren method(密度勾配), shadowgraph(密度の 2 階微分) がある．
5 乱流境界層の中ではほぼ同じ乱れ現象が繰り返される．1 組の熱線プローブを基準位置に固定し ，他の 1 組の熱線プローブを移動する．移動プローブのデータを，基準プローブから得られる情報のある特定の瞬間に合わせて測定すれば，ある瞬間の流速分布す
なわち乱流構造を求めることができる．
6 無次元ヘリシティ H = u  = u  ，２次元流れでは H = 0，また H = 0 は３次元的縦渦の発生またはその兆候を意味する．
4 煙などで可視化するトレーサ法，水素気泡法，また高速流ではj  j j jj j

6

4

乱流と乱流解析

れによるもので，大スケールの乱流構造は輸送に大きい役割を果たす．例えば 部分加熱面上を過ぎ る流れを考えれば ，層流境界層では，熱は分子拡散によるので境界層の外縁近くに達するのは 容易でないが，乱流境界層で境界層厚さ規模の渦があれば ，熱は境界層外縁近くまで簡単に輸送されることになる．乱流構
造すなわち組織的 (coherent) に運動する渦塊は，熱などのスカラー量だけでなく小規模渦も輸送し ，更にその周囲の流れにも大きい影響を及ぼす．乱流構造は乱流の中で重要な役割を果たすものであるが，その運動は遅く激しく変動する小スケール渦に隠され多くの流れでは発見が遅れた．
次に DNS によって解明された組織構造のいくつかを紹介する．滑らかな壁面上の乱流遷移についてはすでに述べたが，折れ曲がりのある壁面ではどのようになるのであろうか．例えば流れに平行に置かれた長方形柱周りの２次元的流れでは，ある実験研究によれば流れは角から剥離しその下流側で再付着しここに剥離泡が形成される．しかし宇宙研の桑原ら (1984{ ) による同様の流れの DNS の結果は，角から横渦が次々に発生し ，始めは渋滞気味であるが，少し 下流にいくと車間距離が開いて速く走れるようになるというものである．同じ乱流境界層内の横渦の転がりでも，渦間距離の大小によって時間平均流が剥離したりしな
かったりするのである．なお この流れは２次元計算で得られたもので乱流の DNS ではない．次に一様流中に置かれた円柱周りの流れについて述べる．Re = 104

 105 の流れは，初期の計算では抗力係数が実験とかなり良く一致していたが，コンピュータが高速化し 計算点が増えスキームの精度が高くなるにつれて抗力係数が過大に計算されるようになった．これは縦渦の発生を無視したためで，一方 初期の計算では，格子が粗く 1 次上流差分または大きい人工粘性を用いたため，その数値粘性が縦渦による拡散の代役を務め，後流の広がりを抑え，怪我の功名でたまたま良い結果が得られたと言うことである．３次元計算で縦渦まで詳細にシミュレーションすればその抗力係数は実験と良く一致する．更に Re = 4

 105 の臨界レ イノル

ズ数近辺の微妙な抗力の変化も実験と同じように出すことができる．凹面上の境界層内には G
ortler 渦の形
成されることが知られている．円柱の下流側に凹面はないが，薄剪断層の流れが凹面上と同じように湾曲して流れるので縦渦構造が現れるのであろう．
次に自由剪断乱流 (free turbulent shear layer) について述べる．流速の異なる２つの並進流れの境界には混合層 (turbulent mixing layer) が形成される．この境界面の不安定性，渦塊 (roller) の形成，渦塊のペアリングによる合体成長は離散渦法で２次元計算されている．差分法では更に２次元的横渦に絡まる縦渦(rib) や横渦の変形が３次元計算されている．圧縮性流れの混合層で両側の流速差が大きくなり，そのマッ
ハ数が 2 に近付くと混合層は急に発達しなくなる．これは，渦塊から上下に衝撃波 (shock waves) が伸び，流れの中の情報は超音速域を越えて遡ることができないので，渦の相互干渉が止むからである．円形ジェット流の表面は混合層でこの流れでは横渦は環状になる．この混合層は平面混合層のように安定ではなく，環状渦はすぐ に非対称になり，崩壊してしまう．鍋の中の水や日の出後の地表の大気では，Benard セルや上昇気流という乱流構造によって熱は平均温度の低い方から高い方へ運ばれ，熱の逆勾配乱流輸送が起きる．
このような浮力乱流に対しては特別のモデル (sophisticated model) が必要である．
次節以降においては乱流の計算とその際必要になる乱流モデルについて述べるが，乱流を計算するには
乱流そのものをある程度知ることが必要である．
問 環状流路の流れについて，断面積一定，外径が流れ方向に増加する場合に，この流れの乱れはどのように変化するか考えよ．

5

乱流と乱流解析

10.2

乱流の計算

以上見てきたように乱流は，その乱れエネルギーの大部分を占め乱雑に運動している小規模渦と，輸送に重要な役割を果たす乱流構造からなる．統計理論は乱流構造を記述することが困難であり，一方 直接数値
シミュレーション (DNS) は膨大な計算を必要とし しかも扱える流れは限られる．現時点で乱流を計算する

一般的方法は，古くから研究さられてきた統計的手法によるものである．統計的手法では，NS 方程式の代

わりに，流れの変数をその時間平均値と変動成分に分け，これを NS 方程式に代入することによって得られる時間平均 NS 方程式が解かれる．この式には新たな変数としてレ イノルズ応力 (Reynolds stress) が現れるので，これを何らかの方法で決定しなければならない．簡単に代数式で与えるものから，レ イノルズ応力の輸送方程式を解くものまで各種の戦略がある．簡単なものは適用限界が厳し く，一方普遍性に富むものは複雑で計算量も多くなる．レ イノルズ応力など の輸送方程式は NS 方程式から導かれる．しかしながら これらの式にはまた新たな変数が現れるので，未知変数の数を式の数に合わせ方程式系を完結させることが必要である7 ．物理数学的考察を基に未知変数の数を減らし ，またその際に導入される係数を決定することが必要である．このようにして作られた乱流の計算に用いられる数式のことを乱流モデル (turbulencemodel) といい，また乱流モデルを作ることを乱流のモデリング (modeling) という．
非圧縮性流れの NS 方程式は次のように書くことができる．
@ui + @ u u = ; @p + F + @ n u + u ; 2 u oij ji 3 ij kk@t @xj j i@xi i @xj

(10.4)

ただし uij = @ui =@xj である．ここでは Boussinesq 近似を置くことにする．すなわち密度変化は浮力 Fiとして考慮するが，流れは非圧縮性とする．変数を時間平均値と変動成分の和すなわち

ui ui + ui

(10.5)

0のように表し ，式 (10.4) に代入し 時間平均を取れば ，次のレ イノルズ方程式 (Reynolds equation) が導か
れる．



@
@
@p
@  ;
@t ui + @xj uj ui = ; @xi + F i + @xj uij + uji ; uiuj
0

ただし

0

(10.6);ui uj は応力の次元をもつ量でレ イノルズ応力 (Reynolds stress) と呼ばれる．この式の誘導に際00

しては，変動成分の平均値がゼロすなわち u0i = uj u0i = 0，と連続の条件 ukk = 0 が用いられた．また非定常流れの場合には時間平均を取ることができないので，アンサンブル平均 (ensemble average) 8が取られる．一般に変数を式 (10.5) のように分解することをレ イノルズ分解 (Reynolds decomposition)， の付いている量をレ イノルズ平均 (Reynolds average) という．レ イノルズ方程式 (10.6) は NS 方程式と類似のもので，左辺は対流項で平均流速の平均流による輸送を表し ，右辺は圧力項，浮力など の力の項，拡散項を示している．拡散は分子粘性とレ イノルズ応力によるが，乱流では後者が支配的である．このレ イノルズ
応力の拡散項は，NS 方程式の対流項に由来するものであるが，物理的には乱れによる拡散を表すものである．レ イノルズ応力を何らかの方法で決定できれば ，NS 方程式の代わりにレ イノルズ方程式を使って平均

close

second-moment closure

7 完結するはで，例えばのように使われる．
8 例えば バルブを開けて管路に水を流す場合に，バルブを開けてからある時間後のある点における流速のアンサンブル平均値は，バルブを開け管路に水を流し流速の瞬時値を測定し ，これを繰り返し平均値を取ることによって求めることができる．

6

乱流と乱流解析

流を求めることができる．レ イノルズ応力を決定する方法は渦粘性型モデルとレ イノルズ応力輸送方程式
モデル (普通 略して応力方程式モデルといわれる) に大別できる．

まず渦粘性モデル (eddy viscosity model) について説明する．歪み速度テンソル (rate of strain tensor)

S 9 は膨張圧縮に関係する等方 (isotropic) 成分と粘性拡散に関係する非等方 (anisotropic) 成分に分けることができる．

2
3
2u11 u21 + u12 u31 + u13
6
7
2S = 64u12 + u21 2u22 u32 + u2375

u13 + u31 u23 + u32 2u33
3
2
3 2
1 0 0
2u11 ; (2=3)ukkuu21 + u12
31 + u13
7
6
7 6
= 32 ukk64 0 1 0 75 + 64 u12 + u21
2u22 ; (2=3)ukk
u32 + u23 75
2u33 ; (2=3)ukk
0 0 1
u13 + u31
u23 + u32で表せば ，式 (10.4) の粘性項は

この式の非等方成分をはが定数で連続の条件r

のようにベクトル表示できる．この粘性項r u = 0 が満足されれば r2u になる．レ イノルズ応力テンソル R もその等方成分と非等方成分に分けることができる．

2 2
3
;
u
;
u
u
;
u
u
1
1
2
1
3
7
6
R = 64;u2 u1 ;u22 ;u2 u375
0

0

0

0

0

0

0

0

0

0

;u3 u1 ;u3 u2 ;u32
3
2
3 2 2
1 0 0
;
u
+(2
=
3)
k
;
u
u
;
u
u
1 2
1 3
7
6
7 6 1
= ; 32 k64 0 1 0 75 + 64 ;u2u1;u22 +(2=3)k;u2 u3 75;u3u2;u32 +(2=3)k
0 0 1;u3u1
0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

ただし k = (1=2)u0k u0k は乱れの運動エネルギー (turbulent kinetic energy) である．Boussinesq の渦粘性近似は，レ イノルズ応力の非等方テンソルが歪み速度の非等方テンソルに比例するとするものである．

;uiuj = t (uji + uij ) ; 23 ij k(10.7)

t は渦粘性係数 (eddy viscosity) と呼ばれる比例定数である．なおこの式では非圧縮性 ukk = 0 とし= ( + t) ;ている．いま レ イノルズ方程式 (10.6) の拡散項を r のように表せば ，応力テンソルは

ただし(2=3)kII となる．ただし I は identity である．またij = ( + t) uji + uij ; 32 ij kの成分は(i j = 1 2 3)(10.8)となる．これは層流の場合とあまり変わらないコンパクトな表現といえる．乱流境界層は層流境界層に比べ発達 (境界層厚さの増加) が著しく，後流やジェットは通常乱流でその発達も著しい．このことは層流の拡散係数に比べ著しく大きいt を用いれば表現でき，その意味で渦粘性近似は一義的に正しい．しかし

9 変形テンソル U は対称の歪み速度テンソル S と反対称の回転テンソル  に分けることができる．

2
u1 1
4u2 1
u3 1

u1 2
u2 2
u3 2

3

2

u1 3
1 2u1 1
u2 3 5 = 4u2 1 + u1 2
2 u3 1 + u1 3
u3 3

u1 2 + u2 1
2u2 2
u3 2 + u2 332

u1 3 + u3 101
u2 3 + u3 2 5 + 4u2 1 ; u1 2
2 u3 1 ; u1 3
2u3 3

u1 2 ; u2 10

u3 2 ; u2 33

u1 3 ; u3 1
u2 3 ; u3 2 507

乱流と乱流解析

例えば Franke-Rodi-Schonung(1989) また最近の数値的研究によれば，実験データや乱流データベースから逆にt を求めるとその値は成分ごとに異なり，かなり広い範囲で負になることが指摘されている．円柱の後流で，平均流として時間平均流に周期変動 (渦列の影響) を加えた位相平均流に対しても同じことが言える．このように渦粘性係数近似にはかなりの難点もあるのである．
圧縮性流れでは，渦粘性型モデルの応力テンソル
置かれる．



とこのモデルに相当の熱流束 q の成分は次のように

ij = ( + t) uji + uij ; 32 ij ukk ; 32 ij k(i j = 1 2 3)(10.9)
 @c2(i = 1 2 3)(10.10)qi = ;  ;1 1 Pr + Prt @xti
ただし  は比熱比，c は音速で，空気に対しては Prandtl 数 Pr = 0:72，乱流 Prandtl 数 Prt = 0:9 である．
渦動粘性係数 (eddy kinematic viscosity) t (= =) は 速度] 長さ] の次元を持つ．Prandtl の混合長理論 (Prandtl's mixing length theory) では，この係数は速度勾配に比例するものとして (勾配拡散モデル )，次のように置かれた．

 
t = l2  @u@y (10.11)

ただし u は平均流の速度，y はそれに直交する方向の距離である．また l は混合長ないしは乱流を特徴づけ
る特性長さで，van Driest(1956) は境界層に対し次のように置いている．

n + ol = y 1 ; exp ; yA

(10.12)ただし  = 0:41 は Karman 定数，y + = yu = は壁面からの無次元距離を表す壁座標 (wall coordinate)，pu = w = は摩擦速度，w は壁面剪断応力，また通常 A = 26 である．この式は Prandtle-van Driest の

f g は，壁面近くの粘性底層 (viscous sublayer) と移行層 (緩和層 bu er layer) に対する van Driest の補正で，慣性底層 (対数領域 inertial sublayer, logarithmic region) ではこの式は l  y式といわれる．この式のとなる．乱流遷移点は別途与えられ，この乱流モデルはその下流側に適用される．
渦粘性係数 t は，乱れの特性速度 q =ppk と特性長さ l に共に比例する量で次のように置くことができる．

t = C k l
広く用いられている k(10.13); " モデルでは，乱れの特性長さ l に代わるものとして，乱れエネルギーの散逸率(dissipation rate) " = uik uik が用いられ，渦粘性係数は次のように置かれる．
00t = C k2="

(10.14)

このモデルでは k と " の値は k と " のそれぞれの輸送方程式を解いて決定される．このように t の値を２つの輸送方程式を解いて決定するものを２方程式乱流モデル (two-equation turbulence model) といい，また１つの輸送方程式を解くものを 1 方程式乱流モデル (one-equation turbulence model)，代数式から計算するものを代数モデル (algebraic model) または０方程式モデルという．kに次の壁法則 (law of the wall) が併用される (図 10.2)．

u+ = 1 ln y+ + 5:15

; " モデルでは，壁面のところ(10.15)

8

乱流と乱流解析

ただし u+ = u=u は無次元速度である．このモデルでは壁面に隣接する格子点は，この法則が成立する壁面から少し離れた慣性底層内に取られることになり，格子点数を節約できる．しかしながら流れによっては壁法則が成立しないことに注意しなければならない．壁面まで適用できる k型k; " モデルは低レ イノルズ数; " モデル (low Reynolds number k ; " model) と呼ばれる．このモデルでは壁面の隣接格子点を粘性底層内に取らなければならないが，境界層流れは層流でも剥離していてもよく，また上流側に小さい k と "の値を与え，バイパス遷移 (by-pass transition) を期待することもできる．バイパス遷移は，一般には主流の乱れによって層流境界層内の流れが不安定になり乱流に遷移することをいう．
問1

渦粘性係数が式 (10.14) のように置かれる理由を次元解析をもとに示せ．

問2

同じ Karman 定数  が式 (10.12) と式 (10.15) にこのように現れるのはなぜか考えよ．

図 10.2: 乱流境界層の平均速度プロフィル

9

乱流と乱流解析10.3

乱れ量の輸送方程式

レ イノルズ応力の輸送方程式は NS 方程式 (10.4) を基に導かれる．(10.4) の i 成分の式に乱れの成分 u0j

を掛け平均したもの u0j (10.4)i は次のようになる．

@ u + u @ u u = ;u p + u F + u @  (u + u )uj @t
iik kij @x k i
j i
j i
j @x
0

0

0

k

0

0

k

この式を u0j = 0 ukk = u0kk = 0 などの関係を用い書き換えれば@ u + u u u + u u u + u u u
uj @t
k j iki
j k ikj k ik@ pu +pu +u f + @ u u ; u u= ; @xji
j i @xk j ikik jk
i j
0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

(10.4) の j 成分の式に乱れの成分 ui を掛け平均したもの ui  (10.4)j も同様になるからこれらの式を加え合わせれば ，次のレ イノルズ応力の輸送方程式 (transport equation for Reynolds stress) が導かれる．
00@ u u +  @ u u u = ;u u u ; u u u + u f + u f + p (u + u )j k iki k jkj i i jji ij@t i j @xk k i j
対流 cij
生成 PijFij
圧力 歪み相関 ij; @x@ ui uj uk + p ( kj ui + ki uj ) ; @x@ ui uj ; 2 uik ujkkk
拡散 dij
散逸 ;"ij
0

0

0

0

0

0

0

]

0

0

0

0

0

0

0

]

0

0

0

0

0

-

0

0

0

0

0

]

]

0

(10.16)]式 (10.16) は，レ イノルズ応力が平均流によって対流輸送されるときに，その生成 (production, 生産)，圧
力-歪み相関 (pressure-strain correlation)，拡散 (di usion)，粘性散逸 (viscous dissipation) によって変化することを示している．次に各項の解釈とモデル化について簡単に述べる．
生成項 Pij： この項の意味を，先に示した渦度輸送方程式 (10.3) の渦度の発生項と対応付けて説明する．

流れの中に取られた仮想的面に作用するレ イノルズ応力はレ イノルズ応力テンソル R にこの面の外向き単

位法線ベクトル n を掛けたものである． レ イノルズ応力]= n R = R n ．渦度  の輸送方程式の発生項は( r)u であった．したがって これに相当するレ イノルズ応力 n R の輸送方程式の生成項は (n R r)u となる．n は任意べクトルであるから，レ イノルズ応力成分 u0i u0j の輸送方程式の生成項は u0j u0k uik となる．
またレ イノルズ応力の対称性により発生項には u0i u0k ujk も含まれるべきである．渦度とレ イノルズ応力は同じ現象を別の角度から見ていると考えれば ，レ イノルズ応力ベクトルも流体に凍結され流体とともに移流し，伸縮によってその強さが増減することは当然のことと言える．しかしながらレ イノルズ応力は磁場や渦度とは多少異なる性質の方程式で支配されるので全く同一に論ずることはできない面もある．この項のすべてがレ イノルズ応力の生成である必然性もないともいう人もいる．生成項はモデル化する必要がない
のであまり論じられていない．
圧力-歪み相関項

ij ： 始めにこの項の圧力変動 p を圧力の Poisson 方程式 r2 p = ;r (u ru ) を用
0

い，流速の式に置き換える．この方程式は Green の定理を適用し積分すれば次のようになる．

Z ;
p = 4 @x@ u^mu^lm dVrl10

乱流と乱流解析

ただし r は p の点から体積要素 dV までの距離，\ ^ " の付いている量は dV における値である．また積分は乱れのマクロなスケールにわたって取れば 十分である．広い範囲が関係するという反駁もあるがここの議論には関係ない．これより次式が得られる．

 Z @ 2u^l u^m (u + u ) dV +  Z 2u^ u^ (u + u ) dV
=
p
(
u
+
u
)
=
ij
lm ml ji ij r
ji ij4 @xl @xm ji ij r 4ij 1 + ij 2
0

0

0

0

0

0

0

0

0

0

(10.17)

圧力-歪み相関項はモデル化しなければならない．一般に数学モデル (mathematical model) は，物理現象を反映するものでなければならない．乱れは，剪断流れの横渦が３次元的渦構造に発達する過程で発生し ，平均流によって対流輸送される間に，渦の干渉によって生成消滅し，等方化し，拡散する．渦は互に干渉し細分化されるが，大きいスケールの渦から小さいスケールの渦へ乱れエネルギーが伝達されることをエネ
ルギー・カスケード (energy cascade) と言う．渦はまた散逸しついには熱エネルギーになる．Rotta(1951)
は乱流のこのような物理のうち乱れの等方化過程をモデルを提案している．

"2

ij 1 = ;c1  k ui uj ; 3 ij k00

ij 1 に割り当て，次の線形等方化 (return-to-isotropy)(10.18)ただし c1 は正の係数，"=k は次元を合わせるために導入されたもの，u0i u0j; (2=3) ij k はレ イノルズ応力の非等方成分である．このモデルは，レ イノルズ応力の法線成分 (ui ui  i = 1 2 3) をその平均値 (2=3)k00= uk uk =3 に近付け，また剪断成分 (ui uj  i 6= j ) を 0 に近づけ，レ イノルズ応力を等方化するものである．
一般にモデルはまた，次元，ランク，対称性，縮約 (contraction) などの条件を再現するものでなければならない．式 (10.18) は，他の項と同じ U 3L の次元を持ち，またモデル化前の式と同様に乱れ速度の３次式になっている． ij 2 のモデルの最も簡単なものは次の Naot ら (1970) の生成の等方化 (isotropizaton of
production, IP) モデルである．
 2 (10.19)ij 2 = ;c2 Pij ; 3 ij P
ただし P = Pkk =2 である．この式は式 (10.18) と類似の性質のものである． ij 2 も他の項と同じ 次元を持ち，モデル化前と同様に Pij を介して平均速度勾配の 1 次式，乱れ速度の２次式になっている．このモデル
0

0

0

0

は簡単であるが勝れたモデルである．
圧力-歪み相関項

ij をモデル化したものは再分配項 (redistribution term) といわれる．Launder は，再分配項は課税分配システムのようなもので， ij 1 は富裕税， ij 2 は所得税に譬えられるとしている．一方が高ければ他方は安くなり，多くの研究者の結果は 0:23c1 + c2 = 1 の近辺にある． ij 1 は return 項といわれる． ij 2 はかつて rapid distortion theory を用いてモデル化されたため rapid term といわれ，これに応じて ij 1 も slow trem といわれることがある．最近，DNS や LES により多くの乱流データベースが作られ，それを基にモデル化前と後の各項の値が比較され，モデルの妥当性が論じられている．上記のモデルの問
題点もいろいろと指摘されている．再分配項については壁面近傍の効果を含めるものもあり，多くのモデルが提案されている．なお 実験研究により，乱流に外部から音場をかけると 乱れの等方化が促進されることは以前より知られている．

r a のように書かれる項を持つ輸送方

拡散項 dij ： ここに集められた項はすべて発散形になっている．

11

乱流と乱流解析
程式 du=dt =dZr a はある流体のシステム V にわたって積分すれば ，

Z du

Z

Z @aa

dt V u dV = V dt dV = V r a dV = S @n dSとなる．ただし S は V の表面で，n は面積要素 dS の外向き単位法線ベクトルである．この式は，流体のシ
ステム V において，a がその内部で生成消滅することなく境界からの流入流出によってのみ変化する保存量で，したがってこの輸送方程式では u も保存されることを示している．このことは r a が拡散項であるための必要条件を満たすことになるが，その十分条件は a = ru  > 0 のように書けることである．dij の第 3 項は紛れもなく分子粘性による拡散項である．この項は乱流の中では小さく無視される．また第 2 項の圧力拡散項も小さいとされる．乱流内における乱れによる拡散過程は残る３重速度相関に割り当てられる．
この項は通常 Daly-Harlow(1970) の一般化された勾配拡散仮説によりモデル化され次のように置かれる．


dij = cs  @x@ k" uk ul (ui uj )l
0

k

0

0

0

;



(10.20)r (k=")R rR となる．この表現は，次元を合わ

ただし cs は係数である．この式はベクトル形で表せばせるために導入された k=" を除いて考えれば ，分子粘性拡散でスカラー量  であった拡散係数が，乱れによる乱流拡散ではレ イノルズ応力テンソル R に置換えられることを示している．乱流拡散は異方性を持つが，この式はそのことをある程度説明するものである．最近，圧力拡散項は必ずしも小さくないことが分かり，この式は圧力拡散項を含めてモデル化されたものと解釈されるようになっている．拡散項に関しても多くのモデルが提案されている．
散逸項 "ij ：渦は干渉によって細分化・等方化し，散逸し遂には熱エネルギーになる．この考えでは散逸するときにレ イノルズ応力の法線成分は等しく，剪断成分は 0 になる．したがって"ij = 2 uik ujk =00

 (2=3)  u uik ik (i = j )00= 23 ij "0(i 6= j )(10.21)

ただし " =  u0ik u0ik は乱れの運動エネルギー k = ui 2 の散逸率である．この値は " の輸送方程式を解いて0決定される．
以上をまとめればレ イノルズ応力輸送方程式はモデル化の後に次のようになる．

@ u u +  @ u u u = P + F ; c  " u u ; 2 k ; c P ; 2 P 
2 ij 3 ij@t i j @xk k i j ij ij 1 k i j 3 ij
+ cs  @x@ k" uk ul (ui uj )l ; 23 ij "
0

0

0

0

0

0

k

0

0

0

0

(10.22)" の輸送方程式も NS 方程式を基に導かれる．2uil (10.4)il を作り整理すれば@@@@t " + @xj uj " = ;2 (uil ujl + uli ulj )uij ; 2 @xi pl uil +   
0

0

0

0

0

0

0

のようなかなり長い式になる．これをモデル化すれば，レ イノルズ応力の輸送方程式と共に解かれる次の "の輸送方程式が得られる．

@ " + @ u " = " (c P + c F ) + c  @  k u u "  ; c  "2
"3
" @x " i l l"2 k@t
@xi i k "1i
00

(10.23)

12

乱流と乱流解析

しかしながらこの方程式は，実際には次に示す k の輸送方程式とアナロジーを持たせて作られたものともいえる．最近の DNS の研究によれば " は剪断乱流では異方性が強く，" を決定する方法はレ イノルズ応力
輸送方程式モデルの泣き所である．

LRR(Launder, Reece and Rodi) モデル10 は典型的なレ イノルズ応力の輸送方程式モデルとして知られるもので，このモデルでは上記のレ イノルズ応力の輸送方程式 (10.22) と " の輸送方程式 (10.23) が解かれる．これらの方程式に含まれる係数 (coe!cient) は，c1  c2 がほぼ一様な剪断乱流から，c"1 と c" が壁近傍の流れから，c"2 が格子乱れの減衰と平面噴流幅の拡大率から決定され，かつまた cs ，c"1 ，c" がいくつかの基本的乱流に対して最適化される．その後の研究成果によりこれらの係数の推奨値は例えば次のようになっている．c1 = 2:0 c2 = 0:6 cs = 0:22 c"1 = 1:44 c" = 0:18 c"2 = 1:92．

k ; " モデルでは乱れエネルギー k とその散逸率 " の輸送方程式が解かれる．k の輸送方程式は上記のレイノルズ応力の輸送方程式から導くことができる．すなわち式 (10.16) の３つの法線応力の式を加えたもの

を 2 で割れば次の k の輸送方程式が得られる．

@ k + @ u k = P ; @  1  u u 2 + p u ; k  ; "
i
i
@t
@xi i
@xi 2 i j
輸送 ck
生産
拡散 dk
散逸
0

]

ただし P =

]

0

0

]

(10.24)0]; ui uj uij である．このモデルには乱れの非等方性は考慮されない．したがってこのモデルに
00

は圧力-歪み相関項がない．式 (10.24) はモデル化すれば次のようになる．

@ k + @ u k = P + @  t k  ; "(10.25)@t@xi i@xi k i
また " の輸送方程式はこの k の輸送方程式とのアナロジーから作られ，次のようになる．
@ " + @ u " = C " P + @  t "  ; C  "2
(10.26)
"1 k
"2 k
@t
@xi i
@xi " ik ; " モデルの係数は Launder-Spalding(1974) では次のように与えられる．C = 0:09 C"1 = 1:44 C"2 =1:92 k = 1:0 " = 1:3．

問

実際の計算で応力方程式モデルが敬遠されるのはその計算量の多さのためである．非定常２次元非圧縮性流れの場合に，Navier-Stokes 方程式に較べ LLR モデルと標準型 k; " モデルの計算に必要なメモリと計算量はおよそど の程度になるのか評価せよ．なお NS 方程式の式の数は運動方程式 2，圧力方程式 1 の計3 である．

Launder, B.E., Reece, G.J. and Rodi, W., Progress in the development of a Reynolds-stress turbulence closure, J. Fluid
Mdch., Vol.68(1975), 537{66.
1013

乱流と乱流解析10.4

乱流モデルとその適用限界

乱流モデルは数多く提案されており，それらは次のように分類される．

22０方程式モデル6666 渦粘性モデル―4 １方程式モデル
2２方程式モデル66 レ イノルズ平均モデル―666
666
4 代数応力モデル
66
レ イノルズ応力輸送方程式モデル
66"
66 LES ― 代数モデル
64
輸送方程式モデル
直接数値シミュレーション

この表は，概して上のものほど 計算量が少ないが，適用できる乱流は限られる．逆にラージ・エデ ィ・シ
ミュレーション (LES) や直接数値シミュレーション (DNS) は一応どのような乱流にも適用できる．しかしながら航空機の翼の流れを解析する場合に，LES や DNS は翼前縁近くの乱流遷移を局所的にシミュレー
ションすることはできても，翼まわりの流れ全体をシミュレーションすることは現実には不可能である．普通この目的には渦粘性モデルが用いられる．この節ではレ イノルズ応力など の相関量にはReynolds 平均値のを付けるが ，

は省略することにする．

渦粘性モデルは０方程式モデルから３方程式モデルまであるが，良く使われるものは０方程式モデルと２方程式モデルである．０方程式モデルは代数モデル (algebraic model) とも言われ，非圧縮性流れの境界層に対して開発された Cebeci-Smith モデル 11と圧縮性流れの Baldwin-Lomax モデル 12が良く知られている．これらのモデルは２層代数渦粘性モデルといわれるもので，境界層は内層と外層に分けられ，渦粘性係数はそれぞれ別の代数式で与えられる．

t

 ( t )in;
= min ( )  ( ) =
t in

t out

( t )out(y  ycross )(y > ycross )(10.27)

ただし y は壁面からの距離で，ycross は ( t )in = ( t )out になるところの y である．これらのモデルでは基本的に境界層の内層には Prandtl-van Driest の式 (10.11)，(10.12)，また外層には次式が用いられる．

Z 



( t )out = CCl  (ue ; u)dyFKleb(y)(10.28)0

ただし CCl = 0:0168 は Clauser 定数，ue は境界層外縁の u の値で，また u(y )0 ならば積分値は ue とKleb (y ) は Klebano の間欠係数で次式Rなる． = (1=u ) (u ; u)dy は境界層の排除厚さである．また F


e 0

e



で与えられる．

FKleb(y) = 1+5:5(y= )6] 1;

は境界層厚さで u = 0:995ue になる y である．外層では FKleb の間だけ乱流で，1(10.29);FKleb の間は外部ポテ

ンシャル流れになる．

Cebeci, T. and Smith, A.M.O., "Analysis of Turbulent Boundary Layers.\ 1974, Academic Press.
Baldwin,B.S. and Lomax, H., Thin layer approximation and algebraic model for separated trubulent ows. AIAAPaper, 78{275(1978).
11
12

14

乱流と乱流解析このモデルでは乱流から層流への遷移点を与える必要があり，CS モデルでは次の Michel の経験式が満足されるときに乱流になるものとしている．

R > 1:174 (1+22 400=Rx)Rx0:46(10.30)

ただし R = ue = ， = (1=ue2 ) 0 u(ue; u)dy は境界層の運動量厚さ，Rx = ue x= ，x は平板の前縁からRの距離である．０方程式モデルは壁面上のある位置での情報のみを含み，上流の効果が入らない．このモデルは乱れエネルギーの局所平衡 (local equilibrium) が保たれる流れ，すなわち乱れエネルギーの生成と散逸がほぼ釣合う平板境界層流れなどに良く適用できる．しかしながら CS モデルでは，van Driest の補正のA = 26 を減速域 (adverse pressure gradient) や吹出し壁面で l が大，また逆に増速域 (favourable pressure
gradient) や吸込み壁面で l が小になるように調節し ，局所平衡の制約をある程度緩和している．また層流から乱流への遷移は遷移点で急激に起きるとすることが多いが ，遷移点からある遷移長さを経て起きるよ
うに t に遷移の間欠係数 Ftr (x) を考慮することもある．また低レ イノルズ数効果を含めるために Clauser

定数 CCl = 0:0168 の修正も行われている． この部分は簡単にいえば ，CS モデルには壁面に沿う圧力勾配，
壁面の曲率，滲出し吸込み壁，遷移長さ，低レ イノルズ数効果が組み込まれている．]

Baldwin-Lomax モデルは 2 層代数渦粘性モデルで，内層には Prandtl-van Driest の式が用いられる．
( t )in = l2j j(10.31)f ;exp(;y+=A+ )g は混合長， = ru は渦度である． = 0:4 は Karman 定数，y+ は壁座

ただし l = y 1

標，A+ = 26 である．他方 外層には Clauser の式の代わりに次式が用いられる．

( t )out = 1:6CCl Fwake FKleb (y)(10.32)

ただし CCl = 0:0168 は Clauser 定数，また;Fwake = min ymax Fmax  0:25ymaxudif 2=Fmax
FKleb(y) = 1+5:5(0:3y=ymax)6 ] 1

(10.33a)(10.33b);ymax は F (y) = y j  j f1 ; exp(;y+=A+ )g が薄剪断層を横断して最大値 Fmax になるところの y である．た
だし後流では exp(;y +=A+ ) = 0 とする．また udif = ju jmax ;ju jmin は薄剪断層を横断しての流速の大きさの差で境界層では ju jmin = 0 である．FKleb は Klebano の間歇係数である．遷移については，薄剪断層を横断してt の最大値が( t )max < 14

(10.34)

1

ならば ，その横断面では薄剪断層内の流れは層流で

t = 0 とする．BL モデルは翼まわりの圧縮性流れな

jj

ど の計算に広く用いられている．このモデルは，内層の式で平均流の @u=@y の代わりに  を用い，また外層の式で Clauser の式の ueの代わりに ymax Fmax ，後流の式のudif の代わりに 0:25ymaxudif 2=Fmax を用いている．このモデルは境界層外縁を見つける必要がなく，積分も不要で，外層の式は剥離境界層と後流にも適用できる．
上流の効果は輸送方程式を解くことによってはじめて本格的に考慮することができる．輸送方程式を解
く最も簡単なモデルは Johnson-King(1985) モデル 13 である．このモデルは他の輸送方程式モデルとは異なJohnson, D.A. and King, L.S., A mathematically simple turbulence closure model for attached and separated turbulent
boundary layers. AIAA J., Vol.23(1985), 1684{92.
13

15

乱流と乱流解析り境界層を横切っての最大レ イノルズ応力に関係する１つの量の輸送方程式のみが解かれる．このモデル
は０方程式モデルと１方程式モデルの中間に位置付けられるもので，翼まわりの圧縮性流れの解析に用いられている．１方程式モデルでは，例えば乱れの特性長さ l が代数式で与えられ，乱れの運動エネルギー kの輸送方程式が解かれる．また t や uv の輸送方程式が解かれることもある．Spalart-Allmaras(1992) のモ

デル 14 では渦粘性係数に相当の量の輸送方程式が解かれる．このモデルは局所的性格を持つので FEM や非構造格子法にも適している．
２方程式モデルは最も広く用いられている乱流モデルであり，中でも低レイノルズ数型 k; " モデルが多く; " モデルは乱れエネルギー k = (1=2)ui ui とその散逸率 " =  uik uik の輸送方程式を解くものである．k ; " モデルは Launder-Spalding(1974) 15 など 多くの研究者によって検証され信頼性の高用いられている．k

0

0

0

0

いモデルである．しかしながら k と " の方程式は不安定性を招き易く，ステップのある拡大流路のコーナーのところのように流れの中に乱流渦が放出されるところでは特別の注意を払わないと解が発散する．これはコーナーのところで k と " の値が極端に大きくなるためで，格子を細かくすればこれらの値は一層大きくな
り計算は更に困難になる．また一般に，計算の初期の段階で k や " が負になり解が求まらなくなる．乱流の乱れの特性を表す基本変数は乱れの特性速度 q =pk と特性長さ l であることを前に述べた．次元解析によれば，k  q 2  "  q 3 l;1 で，すなわち k は乱れ速度の 2 乗，" は更にその 3 乗に比例する量で，流れが激し

く乱れるところではこれらの量は極端に大きくなる．これが不安定性の根本原因である．安定化対策として
は，k と " の方程式の生成項の特別な取り扱い，初期段階における人工拡散の付加，収束半径の大きい反復法の使用などがある．２方程式モデルの２つの変数は，原理的には，q m ln で表される２つの独立な量であれば何でも使えることになる．乱れの特性時間は T l=q で，その逆数 !  q=l は比散逸率 (specic dissipationrate) といわれる．Wilcox(1988) は k ; ! モデル16を提案している．また Hutton-Smith-Hickmott(1987) は
q ; f モデル，Mohammadi(1990) は  ; モデルを提案している．f  q=l は乱れの振動数で ! と同じ次元を持つ．また   l=q  1=l2 である．後の２つのモデルはまだ普及していないが，用いられている変数はいずれも q の低次のもので，k; " モデルに見られるような不安定性は全く生じないか生じにくいということである．

k と " の輸送方程式は前節の終わりに示した式 (10.24) と (10.26) である．このモデルには浮力を考慮したものや自由剪断層に適用できるものもある．また壁法則を併用せずに，壁面近傍の低レ イノルズ数流れ域にも適用できるように拡張したものは低レ イノルズ数型 k; " モデルと呼ばれる．低レ イノルズ数型 k ; "モデルは一般に次のように表される．

t = C f k 2 ="

@ n

t o

@ k + @ u k = P +@t@xi i@xi + k ki ; " ; D
@ " + @ u " = C f " P + @ n + t " o ; C f  "2 + E
"1 1 k"2 2 k
@t
@x i
@x
 i
i

p

" = " ; 2 (@ k=@y)2

i"(10.35)(10.36)(10.37)(10.38)t = 0 になる減衰関数 f が導入される．また " の式には生成項
Spalart, P.R. and Allmaras, S.R., A one-equation turbulence model for aerodynamic ows. AIAA Paper, 92-0439(1992).
Launder, B.E. and Spalding, D.B., The numerical computation of turbulent ow. Comp. Meth. Appl. Mech. Engng.,
Vol.3(1974), 269{89.
16 Wilcox, D.C., Multiscale model for turbulent ows. AIAA J., Vol.26(1988), 1311{20.

このモデルでは， t の式に壁のところで
14
15

16

乱流と乱流解析と散逸項にモデル関数 f1  f2 が導入される．C  C"1  C"2  k  " はモデル定数である．D E は等方散逸 " を用いる場合の付加項である．計算格子は壁面近傍で十分細かく分割し ，壁面境界条件はすべりなし(non-slip) とする．上式は f = f1 = f2 = 1，D = E = 0，" = " ならば通常の k ; " モデルの式になる．
低レ イノルズ数型 k ; " モデルは数多く提案されている．Jones-Launder(1972) モデル 17 では，f = exp ; 1+2R:5=50  f1 = 1:0  f2 = ;0:3 exp(;Rt2 )tpC = 0:09 C"1 = 1:45 C"2 = 2:0 k = 1:0 " = 1:3 D = 2 (@ k=@y)2 E = 2t (@ 2 u=@y2)2 ，ただし Rt (= ql= ) = k 2 = " は乱流レ イノルズ数 (turbulent Reynolds number) である．壁面上では k = " = 0である．また Myong-Kasagi(1990) モデル 18では " の方程式が解かれ，D = E = 0，; 3:45  f = 1:0 f = 1 ; exp(;y+=70) 1+ p
1
Rt
 2o
;
n
f2 = 1 ; exp(;y+=5) 2 1 ; 29 exp ; R6t
C = 0:09 C"1 = 1:4 C"2 = 1:8 k = 1:4 " = 1:3，壁面上で " = @ 2 k=@y2 である．
Chien モデルは，一般曲線座標系の圧縮性流れに対しては次のように書き直される．
@ q^t + @ F^ti + D^ +^g = 0
t t
@t @
i

ただし

!

0

!

(10.39)1( + t )@k=@jkkAq^t = J
 F^ti = Jui  D^ t = ; @@ Jgij @i""+ t @ "=@j0g^t = ;J @ "1

P ; " ; D
A
(
C
P
;
C
f
"
;
f
D
)
1
2
3
k

(10.40)i は曲線座標，J は x 座標から  座標への変換のヤコビアン，gij = ik jk は測度テンソルの成分，P =;ui uj @ui =@xj は乱れの生成項，D = 2k=d2 は壁面におけるエネルギー散逸率，d は壁面までの最短距00離である．またモデル定数とモデル関数は次のように与えられる．

C = 0:09  C1 = 1:35  C2 = 1:8   = 1:3
f = 1 ; e C3y+  f = 1 ; 10::84 e (Re =6)2  f3 = e C4 y+
C3 = 0:0115，C4 = 0:5 である．このモデルを用いるときには，壁面の隣接格子点を粘性底層内の y+ = 1  4に取り，また上流側に k と " の小さい値を与えなければならない (何も与えないと全領域 k = " = 0 になる)．上流境界上で主流に k と " の非常に小さい値を与えれば ，これが境界層に外乱として作用しその流れが遷移する．壁面上の境界条件は k = " = 0 である．通常 安定化のため，層流として計算を始めること，kと " の方程式の生成項の大きさを散逸項と比較しながら抑えること，Navier-Stokes 方程式と k ，" の方程
;

;

t

;

式を同時に解くことが必要である．

Jones, W.P. and Launder, B.E., The prediction of laminarization with a two-equation model of turbulence. Int. J.
Heat Mass Transfer, Vol.15(1972), 301{14.
18 Myoung, H.K. and Kasagi, N., A new approach to the improvement of k ; " turbulence model for wall bounded shearow, JSME Int. J. Fluids Eng., Vol.109(1990), 156{60.
1717

乱流と乱流解析k ; " モデルなどの２方程式モデルでは渦粘性近似が置かれた．歪み速度の非等方成分がレ イノルズ応力の非等方成分に比例し ，その比例定数t が乱れの速度と組織構造の大きさによるという考えは，薄剪断層乱流の発達段階のようにレ イノルズ応力の 1 成分

;u v が支配的な流れでは妥当なものである．しかしな
00

がら 大きい縦渦のある流れや旋回流のようにレ イノルズ応力の 2 成分が大きい流れでは不適切なものになる．このような流れは，代数応力モデルまたは応力方程式モデルを解くことによってはじめて良い結果を得ることができる．応力方程式モデルは前に述べたようにレ イノルズ応力の 6 成分 u0i u0j と散逸率 " の輸送方
程式を解くものである．
ここで代数応力モデル (algebraic stress model) 19について簡単に述べる．k; " モデルは上記のように渦粘性近似のゆえに普遍性に欠け，一方 応力方程式モデルは複雑で計算量が多い．代数応力モデルは応力方; " モデルの簡便さで得ようというものである．レ イノルズ応力 ui uj と乱れエ

程式モデルに近い結果を k0

0

ネルギー k の輸送方程式は，略して次のよう書くことができる．

cij = Pij + ij + dij ; "ijck = P + dk ; "(10.41)(10.42)

これらの式はすでにモデル化されているものとする．偏微分方程式 (10.41) は対流と拡散の２つの輸送項を除けば単なる代数方程式である．今 比 u0i u0j =krij がほぼ一定になるものと仮定すればこれらの項は次のように近似できる20 ． rij は定数であるから微分演算子を自由に出入りできる]D u u = r D k = r c 
cij = Dt
ij Dt
ij k
i j




dij = cs  @x@ k" uk ul (ui uj )l = rij cs  @x@ k" uk ul kl = rij dkkk
ただし D=Dt @=@t + ui@=@xi である．これらの式から次の関係が得られる．
0

0

0

0

0

0

0

0

cij ; dij = rij (ck ; dk ) = rij (P ; ")
式 (10.41) はこの関係を用いれば代数式になり，レ イノルズ応力は結局 次式から計算できることになる．

(2=3) ij P
ui uj ; 32 ij k = (1 ; c2) k PPij;;"(1 ; c )00(10.43)1

k と " の値はそれぞれの輸送方程式を解いて求められる．ここでは簡単のため基本の k ; " モデルに対して説明したが，低レ イノルズ数効果を入れることは容易である．
圧縮性流れでは Navier-Stokes 方程式とともにエネルギー方程式，普通 岐点内部エネルギーの保存方程式が解かれる．また非圧縮性流れの熱伝達を伴う流れでは，温度の輸送方程式が解かれる．これらのエネ
ルギー方程式をアンサンブル平均 (ensemble average) した式には新たな変数として乱流熱流束 (turbulent

heat ux) の 3 成分 qti = h ui が現れる．ただし h (= cp  ) は (比) エンタルピーの変動成分である21 ．渦粘性モデルに相当の渦熱拡散型モデル (thermal eddy di usivity model) では
0

0

0

0

 ui = ;at @=@xi
00

algebraic second-moment approach
ASM
Rodi, W., A new algebraic relation for calculating the Reynolds stresses. ZAMM, Bd. (1976), T219{T221.
()
平均値すなわちファーブル平均 (Favre mean) とその変動成分の和に分けられる．詳細は省略する．

(10.44)19
モデルとも言われるが，これはの略でもある．
20
56
21 圧縮性流れのエネルギー方程式をアンサンブル平均する際には，流速と温度 内部エネルギー，エンタルピー が密度の重みつき

18

乱流と乱流解析と置かれる．ただし  は平均温度，また at は渦温度 (熱) 拡散係数 (thermal eddy di usivity) である．温度
場代数モデル (algebraic heat transfer model) では簡単に at = t =Prt と置かれる．ただし Prt は乱流プラントル数で一定である．このモデルは簡単で広く用いられているが，流れと温度場の様相が明らかに異なる場合には，Prt 一定と置くことは不適当で，輸送効果を考慮しなければならない．温度場２方程式モデル(two-equation heat transfer model) では

at = C f k (k=")n ( 2=" )m(10.45)0と置かれる．ただし次元解析から n + m = 1，また 02 は温度バリアンス，" はその散逸率で，これらの値はそれぞれの輸送方程式を解いて決定される．n = 1 m = 0 ならば上の代数モデル，Nagano-Kim(1988)のモデルでは n = m = 1=2 である．レ イノルズ応力輸送方程式モデルに相当のモデルは，乱流熱流束輸送
方程式モデル (turbulent heatux transport equation model) である．このモデルでは乱流熱流束の 3 成分

qti = h ui の輸送方程式が解かれる．この輸送方程式の散逸項 ;"i は代数式で与えられる場合とそれぞれの輸送方程式を解いて決定する場合があり，また浮力項が無視できないときには  2 と " の輸送方程式000

も解かなければならないことになる．
問1

Baldwin-Lomax モデルから渦粘性係数 ANUT(J),J=1,JE を計算する FORTRAN プログラムのサブルーチンを書け．
問2

壁面上 (y = 0) では u = @v=@y = 0 である．したがって k

げた３つの低レ イノルズ数型 k y2 ，"  y ，t  y3 である．ここに挙; " モデルがこれらの壁面漸近挙動を満足するか否かを調べよ．

19

乱流と乱流解析

10.5

LES

ここで節を改めて空間平均乱流モデルの LES(large eddy simulation) について簡単に述べる．LES は，

DNS(直接数値シミュレーション ) と同様に，適宜細かい格子を用い乱流の渦まで非定常３次元流れのコー
ド で計算するものであるが，DNS と異なる点は，計算に用いた格子では解像できない (unresolvable) 小ス
ケール渦の効果をモデル化し，これを非定常３次元計算に考慮することである．乱流渦は，その生成時には流れごとに固有の性質を持つが，相互作用によって細分化されるにつれて普遍的性質を持つようになる．そのため小スケール渦のモデル化は比較的容易で，また格子を細かくすればモデル化に伴う誤差の影響も小さくなる．ところで適宜細かい格子，小スケール渦はどのように定義されるのであろうか．乱れエネルギーをの Fourier 成分に分解し，乱れのエネルギー・スペクトル関数 (energy spectrum)E ( ) を作る．その波数

高波数成分は，乱流全体は非定常非平衡状態であっても，定常平衡状態になる．これは Kolmogorov(1941)

によって指摘されたもので，Kolmogorov の局所平衡 (Kolmogorov local equilibrium) といわれる．
ここで波数の渦に着目すれば ，そのエネルギーはより低い波数の渦から流入しより高い波数の渦に流出し，またその一部が粘性消散しバランスするということである．このエネルギー伝達を支配するパラメータは散逸率 " で，粘性消散を支配するパラメータは  である．したがって，エネルギー・スペクトルはE ( ) = "m  n F ( = d ) のように置くことができる．指数 m n は次元解析から決定され，次の Kolmogorovの第 1 仮説が導かれる．

E ( ) = "2=3  ;5=3 F ( = d)
ただし(10.46)d = ("= 3 )1=4 は Kolmogorov 波数で単に次元解析から得られたものである．またその逆数 ld = 2= d

は Kolmogorov のスケール (Kolmogorov microscale) といわれる．式 (10.46) は単なる仮説ではなくその後実験によりその正当性が証明されている．

まだ図を用意していないので，保原充・大宮司久明編，

数値流体力学，p.214, 1992, 東京大学出版会の図を参照．]

図 10.3: Kolmogorov の局所平衡

20

乱流と乱流解析

図 10.3 は式 (10.46) を円形ジェットの実験データと比較したものである．図から式 (10.46) の成り立つ波数の範囲は，高レ イノルズ数流れほど 広く低波数領域に及んでいることが分かる．この低波数領域では粘性消散  の効果は小さいと考えられるので，エネルギー成分 E を散逸率 " と波数すれば次の Kolmogorov の;5=3 則 (Kolmogorov の第 2 仮説) が得られる．

の関数と見て次元解析

E ( ) = "2=3 ;5=3

(10.47)ただし  = 1:62 は Kolmogorov 定数である． の指数がほぼ;5=3 になる渦の波数の範囲 (図のほぼ直線の部分) を Kolmogorov の慣性小領域 (inertial subrange) という．小スケール渦の上限は慣性小領域内にあるべきで，計算格子の間隔はその上限に近い波長 (Kolmogorov のスケールの 102

 103 倍) に取ることができる．なお実験で得られるのは 1 次元データで，その Fourier 成分は次のように求められる．

Z1

fi (ri ) = 12 ui (xi )ui (xi + ri ) = 22 Fi ( i ) cos( i ri )d iui
ui 0

乱流渦から小スケール渦を分けるために空間フィルタが用いられる．物理量 f (x  t) の局所平均は次のようになる．

f(x  t) =Z1
0

G(r )f (x +r  t) drr(10.48)

またこの式を Fourier 変換したものは次のようになる．

f(  t) = G( )f(  t)
ただし G(r ) と(10.49)G( ) は空間フィルタ関数 (spatial lter function) である．よく使われるフィルタ関数はガウス・フィルタ (Gaussian lter), スペクトル・カットオフ・フィルタ (spectrum cuto lter, sharp cutolter), トップハット・フィルタ (top hat lter，体積平均フィルタ) で次のように定義される．
r2i i )2 
Gaussian lter: G(ri ) = 6 1 exp ; 6ri2  G( i ) = exp ; ( 24ii(
1 (j i j  = i )
sin(
r
=
)
i
i
spectrum cuto lter: G(ri ) =

G
(
)
=
i
ri
0 (j i j > = i )(1= i (jri j  i =2)
top hat lter:
G(ri ) = G( i ) = 2 sin( i i =2)i i0(jri j > i =2)

ただし(10.50)(10.51)(10.52)i はフィルタ幅で，ガウス・フィルタでは格子間隔の 2 倍，他のフィルタでは格子間隔に選ばれる．

フィルタリングによって物理量は f (x  t) = f(x  t)+ f 0 (r  t) に分けられる．平均量 f(x  t) の状態は，これらのフィルタでは，もとの物理量 f (x  t) の凹凸をある程度ならしたものになる．スペクトル・カットオフ・

 0 = 0 @f=@x = @ f=@x
 でフィルタでは f0 = 0 f = f であるが一般にはそうはならない．また一般に fg

6ある．

Navier-Stokes 方程式 (10.1) をフィルタリングすれば次式が得られる．
@ ui = 0@xi
@ ui +u @ ui = ; 1 @ p +  @ 2ui ; @ (u u ; u u )@t j @xj@xi @xj2 @xj i j i j(10.53)(10.54)

21

乱流と乱流解析

図 10.4: ガウス・フィルタ，スペクトル・カットオフ・フィルタ，トップハット・フィルタ
この式の最後の項はフィルタリングによって現れた SGS(subgrid-scale) 項で，() の中は次のように分けられる．

ui uj ; uiuj = (ui uj ; uiuj ) + (ui u0j + u0i uj ) + (u0i u0j ) = Lij + Cij + Rij
Lij  Cij  Rij はそれぞれ Leonard 項，cross 項，Reynolds 応力項と呼ばれる．
Leonard 項： スペクトル・カットオフ・フィルタでは Lij = 0 である．この項はモデル化せずに見積もることができる．

Cross 項： この項はモデル化する必要がある．Lij +Cij の値は必ずしも小さくないがともに無視されることが多い．

Reynolds 応力項： SGS モデルは普通 Reynolds 応力項をモデル化したもののことである．この SGS のレ イノルズ応力は，SGS の乱れによる拡散効果を表すものと解釈され，渦粘性モデルの式 (10.7) にならい
次のようにモデル化される．

u0i u0j = ;t 2D ij + 13 ij u0k u0k 

2Dij = uij + uji(10.55)

ただし t は SGS 渦粘性係数である．このモデル化されたレ イノルズ応力項は，t > 0 ならば，分子粘性項と同様に平均流速を拡散する．しかし LES では，時間平均乱流モデルとは異なり，境界層や混合層の厚さを増す働きをするのは主に非定常３次元流れとして計算された渦で，この項はこの渦を壊す働きをすると考えられる．よく使われる SGS モデルは Smagorinsky モデル 22で，このモデルでは Lij + Cij は無視され，
22

Smagorinsky, J., General circulation experiments with the primeitive equations. I. The basic experiment. Mon. Weather, Vol.91(1963), 99{164.

Rev.

22

乱流と乱流解析t が次のように置かれる．

pt = (CS )2 2D

D = (D ij D ij )1=2

ただし CS は Smagorinsky 定数で，普通 0:1(10.56)

 0:2 の範囲に取られる．壁乱流では小さく，上昇気流を伴う= ( 1 2 3 )1=3 ，また壁の近くでは壁に垂直方向のフィルタ幅が取られる．壁近くではまた低レ イノルズ数効果を入れるために t に van Driest 型

大気の流れでは大きくなる．

はフィルタ幅で直方体格子ではの減衰関数が考慮されることもある．

Smagorinsky モデルは定数が確定しないことが問題である．Germano(1991) らの Dynamic SGS モデル
は，フィルタ幅の異なる２つのフィルタ関数を用い，SGS に近い GS(grid-scale) の情報から Smagorinsky
定数を決定するものである．また Bardina(1983) にはじまるスケール相似則モデルは (最大スケールの)SGS
成分は最小スケールの GS 成分と相似であるとして SGS レ イノルズ応力をモデル化するものである．最近のモデルは Kolmogorov 慣性小領域よりもかなり大きい格子に適用できるものになっている．LES は乱流のデータベースの作成に用いられ，既存の乱流モデルの検証が行われている．この場合には計算領域の形
状はほとんどが直方体で，計算にはスペクトル法という特別の方法が用いられている．しかしながら実用的問題では流れ場の形状は一般に単純でなく，曲線座標格子を用い差分法で，非構造格子を用い有限体積法
で，あるいは任意形状メッシュを用い有限要素法で計算することになる．このような場合の計算例もあるにはあるが，まだその信頼性は十分とはいえない．
問1

式 (10.46)，式 (10.47)，Kolmogorov 波数の式を次元解析によって求めよ．

問2

Smagorinsky モデルを用い SGS レ イノルズ応力 (REY(I,J),I=1,3),J=1,3 を計算する FORTRAN

プログラムのサブルーチンを書け．