２つの乱流スキーム（LESとMYNN）の基本設計	
中西 幹郎	
防衛大学校地球海洋学科	
JAMSTEC RIGC

2012年11月6日	

今日の内容	n 平均流に対する方程式	
► Averaged Navier-Stokes方程式	

n 乱流を計算するスキーム	
► 乱流スキームがモデル化する領域	
► 乱流スキームの計算上の分類と名称	
► ２つの乱流スキーム（LESとMYNN）の決定的な違い	

n LESの詳細と計算例	
n MYNNの詳細と計算例 （一部、新しい結果）	

LES： Large Eddy Simulation

MYNN： Mellor-Yamada-Nakanishi-Niinoスキーム	

2012年11月6日	

平均流に対する方程式	
Navier-Stokes方程式 から	
Averaged Navier-Stokes方程式	

2012年11月6日	

平均流と乱流の分離	
n Navier-Stokes方程式	
∂ui ∂ui u j
1 ∂p g∂ ⎛⎜ ∂ui ∂u j ⎞⎟+=−+ (θ − θ 0 )δ i 3 +ν+
⎜∂t
∂x j
ρ0 ∂xi θ 0∂x j ⎝ ∂x j ∂xi ⎟⎠	
ここで、浮力項（右辺第２項）にはBoussinesq近似を使っている	n 流れを平均流と乱流に分離	
ui = ui + ui'n 移流項（非線形項、左辺第２項）は	
uiu j = ( ui + ui' )( u j + u 'j ) = ui u j + ui u 'j + u j ui' + ui'u 'j

2012年11月6日	

Averaged Navier-Stokes方程式	
n ２つの平均操作	
= ui ,

ui' = 0,ui u j= ui u j ,ui u 'j = 0	：	 ui ≠ ui ,ui' ≠ 0,ui u j≠ ui u j ,ui u 'j ≠ 0► アンサンブル平均	：	 ui► 空間平均	

気象分野では、空間平均の操作もアンサンブル平均と同様とすることが多い	n 方程式に平均操作を加える。平均操作した移流項は	ui u j = ui u j + ui u 'j + u j ui' + ui'u 'j = ui u j + ui'u 'jn 平均流に対するNavier-Stokes方程式	
∂ ui
∂t+

∂ ui u j∂x j' '
1 ∂ pg∂ ⎛⎜ ∂ ui ∂ u j ⎞⎟ ∂ ui u j=−+ ( θ − θ 0 )δ i 3 +ν+−⎜⎟
ρ 0 ∂xi θ 0∂x j ∂x j
∂xi
∂x j
⎝⎠

未知の項 ← ‘渦’粘性項	
2012年11月6日	

乱流を計算するスキーム	
乱流スキームの主目的は渦粘性項を求めること	
・乱流スキームがモデル化する領域	
・乱流スキームの計算上の分類と名称	・LESとMYNNの決定的な違い	

2012年11月6日	

乱流スキームがモデル化する領域	
風速変動のスペクトル	

慣性小領域	

σ 2 = ∫ E dk = ∫ kE

・等方性乱流	・E =αk−5/3dk= kE d (ln k )
k ∫

散逸領域	
MYNN

エネルギー保有領域	
LES

密度
×
周波
数	
4日	

100秒	

1日	
周波数 ⇔ 波数	2012年11月6日	

物理空間で見ると	
大気境界層上端	

エネルギー保有領域の渦	アンサンブル平均≒水平平均：MYNN

空間平均：LES

Δ以下の渦	
2012年11月6日	

Δ
慣性小領域に入るサイズの格子間隔	

乱流スキームの計算上の分類と名称	
n １次乱流クロージャースキーム	
► ２次乱流量を１次の量（平均量）でパラメタライズ	⎛ ∂ u ∂ w ⎞⎟u ' w' = − K M ⎜⎜+∂x ⎟⎠⎝ ∂z

KMは渦粘性係数という	► ほとんどのLESで、このスキームが用いられている	

n ２次乱流クロージャースキーム	
► ２次乱流量の予報式（乱流ui'の予報式にuj' を掛けて平均）で計算	

d ui'u 'jdt= − ui'uk'

∂ uj
∂xk

− u 'j uk'

∂ ui
∂xk

+  −

∂ ui'u 'j uk'∂xk⎛ ∂ u ' v' ∂ u ' w' ∂ v' w' ⎞	ただし、 u ' v' w' = − K q ⎜⎟ のようにパラメタライズ	++⎜ ∂z∂y∂x ⎟⎠⎝► MYNNで用いている	
2012年11月6日	

2次乱流クロージャーの分類と名称	
n 乱流の等方性の程度に応じて分類 （Mellor & Yamada 1974）	
► レベル

１、 ２、 2.5、 ３、 ４	

• レベル４を除き、境界層近似（水平微分は無視）の併用により、 	２次乱流量が１次乱流クロージャーの形に帰着できるのが特徴	
• 本来	d ui'u 'j
dt= − ui'uk'

∂ uj
∂xk

− u 'j uk'∂ ui∂xk+  −	
を解くべきところ	⎛ ∂ u ∂ w ⎞∂u⎜⎟u ' w' = − K M ⎜+⎟ ≈ − K M ∂z
∂z
∂
x⎝
⎠
	のように求めることができる	
► 見かけ上のLESとMYNNの違いはない？	
2012年11月6日	

∂ ui'u 'j uk'∂xkLESとMYNNの決定的な違い	
n 乱流の長さスケールにある	
► LES• 格子間隔Δの中の渦だけをパラメタライズする（空間平均）	
• 乱流の長さスケールは、格子間隔で代表できる	
► MYNN• エネルギー保有渦をもパラメタライズする（アンサンブル平均）	
• 乱流の長さスケールは、格子間隔とは無関係	
– 地上からの距離に依存する長さ 	
Ls– 大気境界層の厚さに依存する長さ 	
Lt– 安定度に依存する長さ		Lb
1/L = 1/Ls + 1/Lt + 1/Lbn 渦粘性係数の‘係数’にある（定数または関数）	
2012年11月6日	

LES（Large Eddy Simulation）	
・支配方程式	
・渦粘性係数	

2012年11月6日	

支配方程式	n 連続の式	

∂u i=0∂xi⎛ ∂u ∂υ ∂w⎞⎜⎜++= 0 ⎟⎟⎝ ∂x ∂y ∂z⎠

オーバーバーが空間平均を表す	

n 運動方程式 （右辺第１項はコリオリ項、分子粘性項は省略）	

∂ui'u 'j∂ui ∂ui u j
1 ∂p g+= −2ε ijk Ω j uk −+ (θ − θ 0 )δ i 3 −∂t∂x jρ0 ∂xi θ 0∂x jn 熱力学方程式（温位保存式）	

∂θ 'u 'j∂θ ∂θ u j+=−
∂t
∂x j
∂x j
n

ui'u 'j と θ 'u 'j が平均操作で現れた２次乱流量	

2012年11月6日	

渦粘性係数	n ほとんどのLESで、１次乱流クロージャースキームを使用		

ui'u 'j = −2K M Sij（厳密には ui'u 'j = −2 K M Sij + 2 eδ ij ）	
321 ⎛⎜ ∂ui ∂u j ⎞⎟S ij =+⎜
2 ⎝ ∂x j ∂xi ⎟⎠θ 'u 'j = −⎛ K ⎞e = ⎜⎜ M ⎟⎟⎝ Ck l ⎠ l は長さスケール	K M ∂θPr ∂x jn 乱流の局所平衡（生成＝散逸）を仮定すると	
1/ 21 g ∂θ ⎞2 ⎛⎟⎟
K M = (CS l ) ⎜⎜ 2Sij Sij −Pr θ 0 ∂x3 ⎠⎝
2012年11月6日	

-5/3乗則と室内実験の結果を利用して求めることができる。
ほぼ定数 ≒ 0.18

LESの計算例	
・Wangara実験（昼間） 	
： 不安定成層	
・Kelvin-Helmholtz波 	： 安定成層	・台風境界層

2012年11月6日			
： 中立成層	Wangara実験の計算条件		
： 250×250×50n 格子間隔 	
： 40m×40m×40mn 時間間隔 	
： 2.0 s （9時～16時）	
n 格子数n 側面境界 	：

周期境界条件	n 下面境界 	： sinカーブで時間変化する顕熱・潜熱輸送を与える		 		： 摩擦速度は粗度長1cmとして計算	n 上面境界 	
： フリースリップ、Rayleigh摩擦層無し	
n 初期値n 地衡風		

2012年11月6日		
： Wangara実験33日目の0900LSTの観測値		
： = 水平方向の平均圧力勾配		
： 鉛直分布を持つが、時間変化はしない	

鉛直流のパターンと高比湿の分布	

2012年11月6日	

大気境界層内の対流	

2012年11月6日	

Kelvin-Helmholtz波の計算条件		
： 270×135×200n 格子間隔	
： 0.1×0.1×0.1 （無次元）	
n 時間間隔	
： 0.01 （無次元時間0～150）	
n 格子数n 側面境界	：

周期境界条件	n 下面境界	： フリースリップ、Rayleigh摩擦層無し	
n 上面境界	： フリースリップ、Rayleigh摩擦層無し	
n 初期値

2012年11月6日		
： U(z) = tanh(z)： Θ(z) = tanh(z) + Γz + Θ(0)： シア層中心（z=0）で、Ri = 0.062012年11月6日	

台風境界層構造の観測例	
水平の平均風向	

風速差	30m/s

Wurman and Winslow (1998)より．	
強風域の平均の間隔は約600m．	
彼らは、鉛直シア流中の熱対流、そうでな	ければ、ヘアピン渦（ストリーク）と考えた	
2012年11月6日	

台風境界層の計算領域	

2012年11月6日	

http://agora.ex.nii.ac.jp/digital-typhoon/より	

台風境界層の計算条件	
n 方程式n ケース	
： 円筒座標に変換し、遠心力の項は近似		
： 目からの距離 100km および 40km	
： 500×500×100n 格子間隔 	
： 40m×40m×40mn 時間間隔 	
： 0.1 s （定常計算、4時間30分積分）	
n 格子数n 側面境界 	：

周期境界条件	n 下面境界 	： 地面からの顕熱輸送は無し（熱対流は考えない）		 		： 摩擦速度は粗度長10cmとして計算	
n 上面境界 	： フリースリップ、Rayleigh摩擦25層	n 初期値2012年11月6日		
： 全層に傾度風35m/s： 地上温位300Kで始まり、温位勾配0.004K/m

方位角方向の変動風v'の水平分布	
R = 40km

R = 100km
方位角方
向	動径方向	

最大風速差はともに約26m/s ← Wurman and Winslow (1998)に近い値	
2012年11月6日	

TKE収支	∂ u '2 / 2∂t+= − u ' w'

∂t

∂z

−

∂ w' u '2 / 2
∂z

− u'

∂ (τ 1 j − τ 1 j )∂x jp' ∂u '+ f u 'υ ' +u 'υ '
ρ 0 ∂x
R= − υ ' w'

∂υ
∂z

2−

∂ w'υ ' / 2∂z− υ'∂ (τ 2 j − τ 2 j )∂x jυ
p' ∂υ '+− f u 'υ ' −u 'υ 'ρ 0 ∂yR∂ w'2 / 2
∂t=g

θ0

< θ ' w' > −p' ∂w' ∂ w' p'+−ρ 0 ∂z∂z
2012年11月6日	

∂ w'3 / 2∂z

シア生成項	
浮力生成項	

υ2

∂ υ' / 2∂u− w'輸送（渦粘性）項	

散逸項	
再分配項	

∂ (τ 3 j − τ 3 j )∂x j

コリオリ・遠心力	

TKE収支の鉛直分布（R=100km）	
動径方向<u'2>成分	

方位角方向<v'2>成分	

鉛直方向<w'2>成分	

S：シア生成、T：乱流輸送、D：散逸、R：再分配、C：遠心力、	
B：浮力生成、P：圧力輸送	
2012年11月6日	

水平平均の風の鉛直分布	
動径方向U成分	方位角方向V成分	

変曲点が存在	

傾度風で規格化した風速	2012年11月6日	

EOF解析による卓越モードのパターン	
重力波モード	

エントレイニングモード	

変曲点不安定モード	

カップリングモード	

2012年11月6日	

MYNN（Mellor-Yamada-Nakanishi-Niino）	
・支配方程式	
・再分配項のパラメタライズ	
・渦粘性係数	

2012年11月6日	

支配方程式	
n 水平一様の扱いから、鉛直１次元の計算で良い	n 表示の簡略化のため、<u> → U、u' → uに置き換え	n 運動方程式、温位保存式、比湿保存式	

平均圧力勾配を表す地衡風	

温度風平衡から得られる水平移流	
2012年11月6日	

2次乱流量の予報式	
n ２次乱流クロージャースキームを使用	
► 例えば、<u2>の予報式	

∂ u2 / 2∂t	= − uw∂U−∂z∂ wu 2 / 2∂z−ν∂u ∂up ∂u+∂x j ∂x j
ρ0 ∂xε/3• シア生成項（右辺第１項）は、パラメタライズ無しに求まる	• 渦粘性項（第２項）は、渦粘性係数を使ってパラメタライズする	• 散逸項（第３項）は、次元解析的に ∝ q3/L
ここで、 q2 = <u2>+<v2>+<w2> = 乱流運動エネルギーの２倍	• 再分配項（第４項）のパラメタライズに乱流モデルの差が出る	

2012年11月6日	

再分配項のパラメタライズ	p ⎛⎜ ∂ui ∂u j ⎞⎟q ⎛
1 2 ⎞+=−uu−q δ ij ⎟⎜ijρ 0 ⎜⎝ ∂x j ∂xi ⎟⎠
3 A1 L ⎝3⎠

等方化の働き	
シアの効果	

⎛ ∂U i ∂U j ⎞⎟
+ C1q ⎜+⎜ ∂x⎟∂xji⎝⎠g ⎛2⎞
− C2uθδ+uθδ−uθδ⎜ i V j 3j Vi33 Vij ⎟Θ 0 ⎝3⎠
2p ∂θ lq=−uiθ lρ 0 ∂xi
3 A2 L∂U i+ C5 ukθ l∂xk
2012年11月6日	

g− C3
θ lθV δ i 3
Θ0

浮力の効果	
追加導入：破線	

新たなクロージャー定数
C2，C3等が加わる
LESのデータベースを利用して決定した	

渦粘性係数	n 通常、レベルは2.5または３を使用	

uw = − K M∂U∂U∂U= − LqSM= − Lq ( S M 2.5 + S M' )∂z∂z∂zwθl = − K H∂Θl∂Θl∂Θ= − LqS H= − Lq ( S H 2.5 + S H' ) l∂z∂z∂zn SM と SH が、LESの定数 CS に対応する関数で、	

Φ2 や D2.5 などは、かなり複雑な関数 （Nakanishi & Niino 2009）	

2012年11月6日		

具体的にレベルとは	n 乱流の等方性の程度に応じて、予報式を診断式に簡略化	

レベル４	

∂ ui u j∂t− Duu = − ui uk∂( ui u j − δ ij q 2 / 3)

0

∂t

∂U j
∂xk− u j uk− Duq = − ui w

∂U ig+uiθV δ j 3 + u jθV δ i 3 + ∂xk Θ0(∂U j∂z− u jw

∂U i 2∂U l+ δ ij ul w+
∂z 3∂zレベル３	

∂q 2∂U∂Vg−0Dqq = −2 uw− 2 vw+2wθV − 2ε∂t∂z∂z
Θ0

レベル2.5∂ θ l2∂t

0 − Dθl = −2 wθl ∂Θl − 2ε θl∂z

∂ qw2

レベル２	∂t0 − Dqw = −2 wqw∂ θ l qw

0 − Dθq = − wθl∂t
2012年11月6日	)∂Qw− 2ε qw∂z∂Qw∂Θl− wqw− 2ε θq∂z∂z

MYNNの計算例	
・Wangara実験（昼間） － LESと同様の条件	
・新潟・福島豪雨	
・Wangara実験（昼夜）	

2012年11月6日	

Wangara実験（温位の鉛直分布）	
MYNN Level 3

赤線はLES2012年11月6日	

MY Level 3

MYNN Level 2.5

Wangara実験（TKE収支）	
MYNN Level 3

赤線はLES2012年11月6日	

MY Level 3

新潟・福島豪雨（気象庁、原さん）	
従来のモデル	

2012年11月6日	

MYNN Level 3の導入後	

Radar-AMeDAS

Hara (25 June 2007; J3.5 AMS会議)より	

Wangara実験（昼夜）の計算条件	
n 物理過程 	
：放射過程を考慮（Katayamaスキーム）		
： 100 （鉛直方向のみ）	
n 格子間隔 	： 20m
n 時間間隔 	
： 10.0 s （9時～翌朝6時）	
n 格子数n 下面境界 	：

地上1.2mで測定された気温を与えて、				： 顕熱・潜熱輸送を計算する		 		： 摩擦速度は粗度長1cmとして計算	
n 上面境界 	： フリースリップ	n 初期値n 地衡風
2012年11月6日		
： Wangara実験33日目の0900LSTの観測値		
： 鉛直分布を持ち、時間変化もする	温位の鉛直分布（昼間）	
Wangara （観測）	

MYNN Level 2.5

下降流の影響？	
2012年11月6日	

比湿の鉛直分布（昼間）	
Wangara （観測）	

2012年11月6日	

MYNN Level 2.5

風速の鉛直分布（昼間）	
Wangara （観測）	

2012年11月6日	

MYNN Level 2.5

温位の鉛直分布（夜間）	
Wangara （観測）	

MYNN Level 2.5

乱流の散逸が大きすぎる？	
MYNNの課題	
水平移流の影響？	

2012年11月6日	

比湿の鉛直分布（夜間）	
Wangara （観測）	

2012年11月6日	

MYNN Level 2.5

風速の鉛直分布（夜間）	
Wangara （観測）	

2012年11月6日	

MYNN Level 2.5

まとめ	n 乱流スキーム	
► n次乱流クロージャースキーム	

• n＋1次の乱流量をパラメタライズするスキーム	
► レベルn• ２次乱流クロージャーに対し、乱流の等方性の程度に応じた簡略化	

n LES► -5/3乗則を利用して慣性小領域の乱流をパラメタライズ	

• 通常、１次乱流クロージャースキーム	

n MYNN► すべての乱流をパラメタライズ	

• ２次乱流クロージャースキームで、通常、レベルは2.5または３	n Terra incognitaの問題	
2012年11月6日	

Terra incognitaの問題	
格子間隔1km程度	

風速変動のスペクトル	

MYNN

LES

密度×周
波
数	
4日	

100秒	

1日	
周波数 ⇔ 波数	2012年11月6日	

ご清聴ありがとうございました	

2012年11月6日