直線直角座標系における大気局地数値モデルの開発
余偉明 沢田雅洋 岩崎俊樹
東北大学大学院理学研究科地球物理学専攻
概要
東北大学地球物理学専攻流体地球物理学講座では斬新な大気局地数値モデルを開発している。それは数値計算手法の発達と強力なコンピューターの出現により、近い将来、O（100ｍ）の水平方向分解能で大気局地数値モデルの実行が現実なことになると考えられるからである。我々は、急峻な地形および複雑な物体をより適切に扱うため、適応な数値的手法を提案し、安定、効率的および正確な力学フレームの開発を目的として行った。これは、地球の表面にある急峻な地形および複雑な物体を、より高い分解能で適切に扱うことが期待されるものである。本研究では、非定常三次元可圧縮ナヴィエ･ストークス方程式 （ Navier-Stokes equations ） を 計 算 す る た め の SIMPLER （ Semi-Implicit Method forPressure-Linked Equation Revised ）アルゴリズムとともに、有限体積法（finite volume method）を用いた。我々は、慣例となっている地形追随規格化（terrain-following normalization）は行なわず、高さを垂直座標として用いるデカルト座標を採用した。地球の平均海面高度より上にある急峻な地形および全ての複雑な物体を扱うために、ブロッキング･オフ法（blocking-off method）を導入した。空間および時間に関する離散化については、高次風上対流法（higher-order upwind convection scheme）を採用し、完全時間陰解法（fullytime-implicit scheme）を利用している。いろいろなテスト計算に対して本モデルを実行した。シミュレーションの結果は、今後の大気局地数値モデル開発において本研究によって提案された数値手法の適用可能性を示すものであった。
1．はじめに大気の流れのモデルに山岳を取り入れるために、地形追随鉛直座標（シグマ）系（Phillips、1957；Gal-Chen および Somerville、1975）が広範に用いられてきたため、現在のほとんどの地域中規模気象数値モデルは、垂直座標として地形追随座標を使用している。しかし、シグマ系の原始方程式モデルの開発において、むしろ初期にも注目を集めていた問題は、運動量方程式における勾配力の 2 つの項にある取り消し不能のエラーである（Smagoringsky 他、1967）。圧力勾配力を表すためのこの 2 つの項は、急峻な地形に対して同等の大きさと逆の符号を有しており、そのため、計算中にそれらの和は重大なエラーを招く可能性がある。Mesinger と Janjic（1985）、中でも Sundqvist（1976）は、シグマ座標を用いた各モデルにおける水平圧力勾配力のコンピューター計算のエラーが、急峻な地形の近傍においてかなり大きくなる可能性のあることを見出した。このエラーを最小化するために、国家環境予測センター（NCEP; National Centers for Environmental Prediction）
中規模イータ･モデル（Mesinger 他、1988）では、イータ座標と呼ばれる段階的山型垂直座標（step-mountain vertical coordinate）を利用している。このモデルでは、地形が、不連続の段階（段階的山型）として表されている。しかし、この段階的山型による表現は、各段階の角（かど）に擬摂動をもたらすことがあり、この表現の正確さは、地形の水平方向の規模、および、垂直グリッドによる実際の地形の分解能に強く依存すると思われる（Gallus と Klemp、2000）。
近年、地形の表現法、すなわち、薄切りセル法（shaved cell approach）、および、これに関連した、高さを垂直座標として用いた海洋および大気のモデルにおける地球物理学の流体方程式に対する数値的手法が、有限体積法に基づき提案されている（Adcroft 他、1997；
Marshall 他、1997；Bonaventura、2000）。コンピューターの数値計算手法の急速な発達と、強力なコンピューターがアプリケーションに使用できることは、分解能が高くなったことにより地形がさらに正確に表現されると思われる近い将来、O（100）ｍの水平方向分解能で地域数値モデルの実行が達成可能になるかも知れないという当然の期待を持たせるものである。このような状況は、現在の高分解能モデル、並びに、さらに高い分解能で実行されることが期待されている次世代モデルでは、地球の表面上の急峻な地形および複雑な物体を扱うためのさらに適切な方法の追求を促進するものである。
本研究では、有限体積離散化に基づくいくつかの先端的な数値的手法（すなわち、完全時間陰解法と高次風上対流法、SIMPLER アルゴリズム、複雑な形状を扱うためのブロッキング･オフ法）を提示し、これらの手法を、次世代大気局地数値モデルのための安定かつ効率的な、および、正確な力学フレームのコーデイングに組み入れる。開発された力学フレームは、より高い分解能において局地規模流れのシミュレーションのために、急峻な地形および複雑な物体をより適切に扱うことが期待されるものである。
2．数値モデルの説明
2.1 支配する方程式
本研究では、粘性可圧縮ニュートン流体に対する非定常三次元ナヴィエ･ストークス方程式を用いた。
運動量方程式：∂ρu ∂ρuu ∂ρuv ∂ρuw
∂p
∂ 2u
∂ 2u
∂ 2u+++= − + ρfv + µ 2 + µ 2 + µ 2
∂t
∂x
∂y∂z
∂x
∂x
∂y∂z
∂ρv ∂ρvu ∂ρvv ∂ρvw
∂p
∂ 2v
∂ 2v
∂ 2v+++= − − ρfu + µ 2 + µ 2 + µ 2
∂t∂x
∂y∂z
∂y∂x
∂y∂z
∂ρw ∂ρwu ∂ρwv ∂ρww∂p∂2w∂2w∂2w+++=− −g+µ 2 +µ 2 +µ 2∂t∂x∂y∂z∂z∂x∂y∂z
エネルギー方程式：∂ρT ∂ρT ∂ρT ∂ρT
∂ 2T
∂ 2T
∂ 2T+++= κ 2 +κ 2 +κ 2
∂t∂x
∂y∂z
∂x
∂y∂z
連続方程式：∂ρ ∂ρu ∂ρv ∂ρw+++=0∂t∂x∂y∂z
理想気体の状態方程式：
P= RT
2.2 時間および空間に関する離散化
従属変数が全て φ によって表されるとすれば、上記の微分方程式は、一般的な形で以下のように書ける。
⎞⎛∂(ρφ ) + ∂ (ρuφ ) + ∂ (ρvφ ) + ∂ (ρwφ ) = ∂ ⎛⎜ K ∂φ ⎞⎟ + ∂ ⎜⎜ K ∂φ ⎟⎟ + ∂ ⎛⎜ K ∂φ ⎞⎟ + S∂t∂x∂y∂z∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠ここで、 φ =(u, v, w, T)

S =(重力、コリオリの力など)
K =( µ , κ ).

(1)

(δx)w(δx)eTｔ(δz)tWΔzｅPｗ

E(δz)bｂ
BΔx

図 1. Control volume for the two-dimensional x-z situation
式（1）を図１に示された制御体積 (control volume)にわたって積分すると、次の式を得る。

ρ φ − ρ 0φ 0P P∆t()

P P ∆x∆y∆z + J − J + J − J + J − J = S + S φ ∆x∆y∆ze w n s t b
C P P(2)

ここで、制御体積面にわたって積分された総フラックスは、次のようになる。

J e = ∫ J x dydz = J x e ∆y∆z ,J w = ∫ J x dydz = J x w ∆y∆z ,J n = ∫ J y dxdz = J yn ∆x∆z ,J s = ∫ J y dxdz = J ys ∆x∆z ,J t = ∫ J z dxdy = J zt ∆x∆y,J b = ∫ J z dxdy = J zb ∆x∆yentwsb

同じ方法で、連続方程式を制御体積にわたって積分することができ、次の式を得る。

ρ P − ρ P0∆t∆x∆y∆z + Fe − Fw + Fn − Fs + Ft − Fb = 0(3)ここで、各界面にわたる質量フラックスは、次のようになる：

Fe = ( ρu )e ∆y∆z ,Fn = ( ρv )n ∆x∆z ,Fw = ( ρu )w ∆y∆z ,Fs = ( ρv )s ∆x∆z ,Ft = ( ρw)t ∆x∆y,Fb = ( ρw)b ∆x∆y
式（3）に φ P をかけて、得られた量を式（2）から引くと、次を得る。
0(φ P − φ P0 ) ρ P ∆∆x∆t y∆z + (J e − Feφ P ) − (J w − Fwφ P ) + (J n − Fnφ p ) − (J s − Fsφ p )+ (J t − Ft φ p ) − ( J b − F b φ P ) = (S C + S Pφ P )∆x∆y∆z

これは、以下のように用いるとJ e − Feφ P = a E (φ P − φ E )J n − Fnφ P = a N (φ P − φ N ) J t − Ft φ P = aT (φT − φ P )J w − Fwφ P = aW (φW − φ P ) J s − Fsφ P = a S (φ S − φ P )J b − Fbφ P = a B (φ P − φ B )三次元離散化方程式は、次のように書ける。

a Pφ P = a E φ E + aW φW + a N φ N + a S φ S + a T φ T + a Bφ B + b φ(3.A)

ここで、a P = a E + aW + a N + a S + aT + a B + a P0 − S P ∆x∆y∆zb φ = S C ∆x∆y∆z + a P0 φ P0a P0 =

ρ P0 ∆x∆y∆z∆ta E = De A( Pe ) + − Fe ,0aW = Dw A( Pw ) + Fw ,0a N = Dn A( Pn ) + − Fn ,0a S = Ds A( Ps ) + Fs ,0aT = Dt A( Pt ) + − Ft ,0a B = Db A( Pb ) + Fb ,0

導通率(conductances)は、次のように定義される。

K e ∆y∆z,(δx )eK ∆x∆z
Dn = n,(δy )nK ∆x∆y
Dt = t,(δz )t

De =K w ∆y∆z(δx )w
K ∆x∆zDs = s(δy )sK ∆x∆y
Db = b(δz )bDw =

また、ペクレ数は次のようになる。

Pe =Fe,De

Pw =FwF, Pn = n ,Dw
Dn

Ps =FsF, Pt = t ,Ds
Dt

Pb =FbDb

関数Ａ（|P|）は、希望する手法に対して表 1 から選ぶことができる（詳しくは、Patanker（1980）、Sha 他（1991）、および、Ferziger と Peric（1997）を参照）。

Scheme( )Formula for A P

Hybrid1 − 0.5 P
10,1 − 0.5 P

Power law0, (1 − 0.1 P )Exponential(exact)

P /[exp( P ) − 1]Central differenceFirst order upwind

5Higher order upwind

QUICK

Table 1 The functionＡ(|P|) for different schemes
これで、時間に関するこの式の積分には、完全時間陰解法が用いられていることが分かる。完全陰時間離散化が用いられているため、時間のステップは、物理的な基準と計算の正確さに関する配慮によってのみ決定することができる。また、空間に関する離散化は、食い違い格子(staggered grid)に対する有限体積技法によって得られ、高次風上対流法は、各制御体積壁面におけるフラックスに関連させて選ばれる。

2.3 SIMPLER の計算アルゴリズム
離散化された方程式における速度と圧力の結合系は、SIMPLER の計算アルゴリズムを用いて解かれる（Patanker、1980）。
一般の変数 φ に対する離散方程式(3.A)に類似して、速度 u、v、w に対する離散化運動量方程式は、次のようになる。

a e u e = ∑ a nb u nb + b u + ( p P − p E )Aenb

a n u n = ∑ a nb v nb + b v + ( p P − p N )An

(4)

nb

at wt = ∑ a nb wnb + b w + ( p P − pT )Atnb

また、推定された圧力場 p * に基づく不完全な速度場 u*、v*、w*に対する離散化運動量
方程式は、次のようになる。

a e u e * = ∑ a nb u nb * + b u + ( p P * − p E *)Aenba n v n * = ∑ a nb v nb * + b v + ( p P * − p N *)An

(5)

nb

at wt * = ∑ a nb wnb * + b w + ( p P * − pT *)Atnb

ここで、正しい圧力と速度が以下より得られることを提案する。

p = p * + p' ,u = u * +u ' ,v = v * +v' ,式（4）から式（5）を引くと、次のようになる。

a e u e ' = ∑ a nb u nb ' + ( p P '− p E ')Aenb

a n v n ' = ∑ a nb v nb ' + ( p P '− p N ')Annbw = w * + w'

at wt ' = ∑ a nb wnb ' + ( p P '− pT ')Atnb

上記方程式の最初の各項を落とすと、次のようになる。

⎛A ⎞⎜⎜ d e = e ⎟⎟ ;ae ⎠⎝⎛A ⎞⎜⎜ d t = t ⎟⎟at ⎠⎝u e ' = d e ( p P '− p E '),wt ' = d t ( p P '− pT '),v n ' = d n ( p P '− p N '),⎛A ⎞⎜⎜ d n = n ⎟⎟an ⎠⎝従って、速度補正式は次のように書ける。

u e = u e * + d e ( p P '− p E ')v n = v e * + d n ( p P '− p N ')(6)

wt = wt * + d t ( p P '− pT ')連続方程式を積分した形は方程式（3）の場合と同じになるが、再び示す。

(ρ − ρ )∆x∆y∆z + {(ρu ) − (ρu ) }∆y∆z + {(ρv) − (ρv) }∆x∆z + {(ρw) − (ρw) }∆x∆y = 0P0Pe∆twnstb

速度補正式で得られた速度の構成要素を、上記離散化連続方程式に代入し、整理すると、p'（圧力補正式）に関する以下の離散化方程式を得る。
a P p P ' = a E p E '+ aW pW '+ a N p N '+ a S p S '+ aT pT '+ a B p B '+b(7)ここで、
a E = ρ e d e ∆y∆z ,aW = ρ w d w ∆y∆z

a N = ρ n d n ∆x∆z ,a S = ρ s d s ∆x∆z

aT = ρ t d t ∆x∆y,a B = ρ B d B ∆x∆ya P = a E + aW + a N + a S + aT + a B
b=(ρ − ρ )∆x∆y∆z + {(ρu *) − (ρu *) }∆y∆z + +{(ρv *) − (ρv *) }∆x∆z + {(ρw *) − (ρw *) }∆x∆y0PPw∆tensbt

さらに、離散化運動量方程式（4）の u 構成要素は、次のように書ける。

∑ a u + b + d (p − p )u =unbenbnbeaePE(8)ここで、擬速度を次のように定義する。

∑ a u +buˆ =nbenbunbae( 9)式（8）は次のようになる。

u e = uˆ e + d e ( p P − p E )同様に、ほかの速度成分も次のように書ける。

v n = vˆn + d n ( p P − p N )wt = wˆ t + d t ( p P − pT )式（7）のように、圧力 p についての離散化方程式が得られ、次のように書ける。

a P p P = a E p E + aW pW + a N p N + a S p S + aT pT + a B p B + b

(10)

ここで、
a E = ρ e d e ∆y∆z ,aW = ρ w d w ∆y∆z

a N = ρ n d n ∆x∆z ,a S = ρ s d s ∆x∆z

aT = ρ t d t ∆x∆y,a B = ρ B d B ∆x∆ya P = a E + aW + a N + a S + aT + a Bb=(ρ P0 − ρ P )∆x∆y∆z + {(ρuˆ ) − (ρuˆ ) }∆y∆z + +{(ρvˆ) − (ρvˆ) }∆x∆z + {(ρwˆ ) − (ρwˆ ) }∆x∆y∆twensbt

流れの場を計算するために、収束する解が得られるまで、式（5）、（6）、（7）、（9）、および（10）を順次に解かなくてはならない。

2.4 計算領域にある不規則な形状の物体の扱い
ここでは、ブロッキング･オフ法によって、任意の形状を扱う方法を述べる（Patanker、1980）。これは、残りの不活性制御体積が、希望する不規則領域を形成するように、規則的なグリッドによって特定の制御体積をブロック･オフすることによって行なわれる。図2 にこの例を示す。同図で、影を付けた部分は、不活性制御体積を示す。任意の形状は、一連の四角形グリッドによって近似されている。

図 2. Blocked-off regions in regular gridブロック･オフ操作の概念は、不活性制御体積内の適切な φ の既知の値を確立することから構成されている。不活性制御体積における希望の値を得るには離散化方程式に大きなソース（source）の項を設定することを含む簡単な方法である。例えば、内部グリッドの点について（すなわち、固い内部において）、方程式 3A の S C と S P を以下のように設定する。

S C = 10 30 φ p,desiredS P = − 10 3030

ここで、 10 は、離散化方程式の他の項が無視できるほど大きい数を示す。その結果、次式を得る。
S C + S Pφ p ≈ 0 ,φP = −SC= φ P ,desiredSP

この手続きは、このような内部境界条件を挿入することによって、計算領域にある不規則な形状の物体の表現に容易に使用できることに注意されたい。

3．結果
まずはテスト計算として、本モデルを平面に設置した立方体を覆う流れの直接数値シミュレーション（Direct Numerical Simulation）に用いた。それは作られたモデルスキームにはバグがないかをチェックするためであった。図 3 は、均一上流速度および立方体の高さに基づいた、レイノルズ数が 800 の中立流れの結果である。図には、障害物の上流から発して、到来する流れによって側面に沿って下流に伸びる馬蹄形の渦を特徴とする流れのバターンを示している。立方体の周辺に数値不安定による流れの発生は見られず、シミュレーションの結果は予想した通りであったので、この力学フレームに関してはバグがないと考えられる。
次に、モデルの性能やスキームの保存性をチェックするため、非線形・非定常性の強いCold Bubble ベンチマークの数値実験を行った。図 4 には Straka et al. (1993) との比較を示している。両者は定量的に異なる部分が見られるが、それぞれの物理量の対応は力学的な矛盾はないことが示された。これらの結果から本モデルにおいて 2 次元重力流の力学特性を精度よく再現できたことがわかった。数値モデルにおいて基本物理量を正しく再現するほかに、計算領域の中で質量やエネルギーが保存することも重要である。図 5 に初期のエネルギーからの変動の時系列を示す。全エネルギーの時間変化を見ると、積分時間 900s で初期の全エネルギーから、およそ 0.01%の差の精度でエネルギーを保存していることが確認された。本モデルのスキームはエネルギー保存性にも優れていることが言える。
図 6 は３次元重力流の先端部付近で起こるロープ・クレフトの渦構造を示している。室
内実験と現地観測で捉えるのに非常に困難とされる重力流の３次元微細構造や力学を明らかにするため、本非静力高精度モデルを用いることによって始めて重力流の３次元微細構造（ローブ・クレフト）を再現することができた。その渦構造はフロント面に沿って周期的に曲がりくねってオメガ型（ω）の渦形状を成していることがわかった。また、その結果について色々な側面から定量的な検証を行った結果、その成長過程と発達メカニズムについても解明された。
複雑な物体を扱え例として、図 7 には都市建物群の形状及び計算格子を示している。こ
こでブロッキング･オフ法（blocking-off method）によって全ての物体形状が簡単に表現することができる。図 8 の計算結果に建物群内の地表近くの水平速度が示されている。図には、物体の角周辺に数値不安定による流れの発生は見られず、シミュレーションの結果は予想した良いものであったので、この力学フレームにおいてそのブロッキング･オフ法に関して有望なものであると言える。現在、該当モデルによる実際の都市建物内の流れに関する計算を進めている。