\underline{2 重周期領域順圧モデル  2 支配方程式系  1}

2 次元 2 重周期境界領域での順圧流体モデルの定式化

竹広真一, SPMODEL 開発グループ

平成 20 年 10 月 10 日

1 はじめに

この文章では, 2 次元 2 重周期境界領域でのβ面順圧流体モデルの定式化とテスト実験の問題設定について述べる.

2 支配方程式系

2.1 支配方程式系

x, y 方向に各々 L_{x}, L_{y} の拡がりを持つ矩形領域での 2 次元非圧縮流体を考える. 支配方程式は, 2 次元β面上での運動方程式と連続の式である:

\frac{∂u}{∂t} + u\frac{∂u}{∂x} + v \frac{∂u}{∂y} − f(y)v = −\frac{1}{ρ} \frac{∂p}{∂x} + F_{x}^{visc},  (1)
\frac{∂v}{∂t} + u\frac{∂v}{∂x} + v \frac{∂v}{∂y} + f(y)u = −\frac{1}{ρ} \frac{∂p}{∂y} + F_{y}^{visc},  (2)
\frac{∂u}{∂x} + \frac{∂y}{∂y} = 0.  (3)

ここで u, v は速度の x, y 成分, p は圧力, ρ は密度であり定数である. f(y) がコリオリパラメター, F_{x}^{visc}, F_{y}^{visc} は粘性散逸項であり, それぞれ具体的な形はここでは定めず一般的に定式化する.

境界条件は x, y 両方向に周期的であるとする.

u(0, y) = u(L_{x}, y), v(0, y) = v(L_{x}, y), u(x, 0) = u(x, L_{y}), v(x, 0) = v(x, L_{y}). (4)

2.2 渦度方程式 (ポテンシャル渦度保存則)

(1) に y 微分, (2) に x 微分を施して差引くことにより, 圧力を消去できて渦度の時間変化の式が得られる.

−\frac{∂}{∂y} [\frac{∂u}{∂t} + u\frac{∂u}{∂x} + v \frac{∂u}{∂y} − f(y)v] + \frac{∂}{∂x} [\frac{∂v}{∂t} + u\frac{∂v}{∂x} + v \frac{∂v}{∂y} + f(y)u]

= −\frac{∂F_{x}^{visc}}{∂y} + \frac{∂F_{y}^{visc}}{∂x} ,

\frac{∂}{∂t} (\frac{∂v}{∂x} − \frac{∂u}{∂y}) + u\frac{∂}{∂x} (\frac{∂v}{∂x} − \frac{∂u}{∂y}) + v \frac{∂}{∂y} (\frac{∂v}{∂x} − \frac{∂u}{∂y})

−\frac{∂u}{∂y} \frac{∂u}{∂x} − \frac{∂v}{∂y} \frac{∂u}{∂y} + \frac{∂u}{∂x} \frac{∂v}{∂x} + \frac{∂v}{∂x} \frac{∂v}{∂y} − f(y) (\frac{∂u}{∂x} + \frac{∂v}{∂y} ) + \frac{df(y)}{dy}v

= −\frac{∂F_{x}^{visc}}{∂y} + \frac{∂F_{y}^{visc}}{∂x} ,

(\frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} ) (\frac{∂v}{∂x} − \frac{∂u}{∂y}) + (\frac{∂v}{∂x} − \frac{∂u}{∂y}) (\frac{∂u}{∂x} + \frac{∂v}{∂y} )

= −\frac{∂F_{x}^{visc}}{∂y} + \frac{∂F_{y}^{visc}}{∂x} ,

(\frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} ) (\frac{∂v}{∂x} − \frac{∂u}{∂y}) + \frac{df(y)}{dy}v = −\frac{∂F_{x}^{visc}}{∂y} + \frac{∂F_{y}^{visc}}{∂x} ,

途中, 連続の式 (3) を用いた. したがって, 渦度 ζ(x, y, t) = \frac{∂v}{∂x} − \frac{∂u}{∂y} を定義すると,

(\frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} ) ζ(x, y, t) + \frac{df(y)}{dy}v = F_{ζ}^{visc},  (5)

ここで F_{ζ}^{visc} = −\frac{∂F_{x}^{visc}}{∂y} + \frac{∂F_{y}^{visc}}{∂x} は渦度の散逸項である. あるいは左辺をまとめて,

(\frac{∂}{∂t} + u\frac{∂}{∂x} + v \frac{∂}{∂y} ) [ζ(x, y, t) + f(y)] = F_{ζ}^{visc}.  (6)

散逸がない場合に ζ + f(y) は保存量となる. q(x, y, t) ≡ ζ(x, y, t) + f(y) は系の絶対渦度あるいはポテンシャル渦度である. これに対して ζ(x, y, t) は相対渦度, f(y) は座標系のもつ渦度^{1}と呼ばれる.

1気象学では惑星渦度 plbaro˙ teishiki.tex  平成 20 年 10 月 10 日 (竹広真一)

2.3 流線関数を用いた表現

連続の式 (3) から速度場を流線関数 ψ(x, y, t) を用いて表すことができる. すなわち,

u = −\frac{∂ψ}{∂y} , v = \frac{∂ψ}{∂x}  (7)

とすると, 相対渦度は

ζ = \frac{∂v}{∂x} − \frac{∂u}{∂y} = \frac{∂^{2}ψ}{∂x^{2}} + \frac{∂^{2}ψ}{∂y^{2}} = ∇^{2}ψ.  (8)

ここで ∇^{2} = \frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}} は 2 次元のラプラシアン演算子である. 流線関数を用いて渦度方程式 (5) を表すと

\frac{∂ζ}{∂t} − \frac{∂ψ}{∂y} \frac{∂ζ}{∂x} + \frac{∂ψ}{∂x} \frac{∂ζ}{∂y} + \frac{df(y)}{dy} \frac{∂ψ}{∂x} = F_{ζ}^{visc}.  (9)

すなわち,

\frac{∂ζ}{∂t} + J(ψ, ζ) + β(y)\frac{∂ψ}{∂x} = F_{ζ}^{visc}, ζ(x, y, t) = ∇^{2}ψ(x, y, t).

ここで J(f, g) ≡ \frac{∂f}{∂x} \frac{∂g}{∂y} − \frac{∂f}{∂y} \frac{∂g}{∂x} はヤコビアン演算子, β(y) ≡ \frac{df(y)}{dy} はコリオリ

パラメターの y 微分である. ポテンシャル渦度保存則 (6) は

\frac{∂q}{∂t} + J(ψ, q) = F_{ζ}^{visc},  (10)
q(x, y, t) = ζ(x, y, t) + f(y) = ∇^{2}ψ(x, y, t) + f(y).  (11)

境界条件は x, y 両方向に周期的であるので,

ζ(0, y) = ζ(L_{x}, y), ζ(x, 0) = ζ(x, L_{y}),  (12)
ψ(0, y) = ψ(L_{x}, y), ψ(x, 0) = ψ(x, L_{y}).  (13)

あたえるコリオリパラメターも y 方向に周期的でなければならない. しかしながら (2.3) には f(y) でなくその微分 β(y) しか方程式にあらわれないので

β(0) = β(L_{y}),  (14)

であれば解くことができる^{2}.

2f(y) が周期的でない場合には圧力を求めようとすると破綻する. この場合ポテンシャル渦度保存則を用いて計算を行なえない.

plbaro˙ teishiki.tex  平成 20 年 10 月 10 日 (竹広真一)

\underline{2 重周期領域順圧モデル 3 数値計算 : 2 重フーリエ級数展開と時間積分 4}

2.4 粘性散逸の表現

粘性散逸項として, 超粘性型のものを考えることにする.

F_{x}^{visc} = (−1)^{p+1}ν_{2p}∇^{2p}u, F_{x}^{visc} = (−1)^{p+1}ν_{2p}∇^{2p}v,  (15)

p は超粘性の次数を表す整数パラメターである. このとき渦度に対する散逸は,

F_{ζ}^{visc} = (−1)^{p+1}ν_{2p}∇^{2p}ζ.  (16)

2.5 支配方程式 (まとめ)

以上, 支配方程式をまとめると,

\frac{∂ζ}{∂t} + J(ψ, ζ) + β(y)\frac{∂ψ}{∂x} = (−1)^{p+1}ν_{2p}∇^{2p}ζ,  (17)
ζ(x, y, t) = ∇^{2}ψ(x, y, t),  (18)

境界条件は

ζ(0, y) = ζ(L_{x}, y), ζ(x, 0) = ζ(x, L_{y}),  (19)
ψ(0, y) = ψ(L_{x}, y), ψ(x, 0) = ψ(x, L_{y}).  (20)

ただし, 与える物理パラメター β(y) は

β(0) = β(L_{y}),  (21)

を満たしていなければならない. β が一定の場合がいわゆる β 面モデルと呼ばれるものである.

3 数値計算 : 2 重フーリエ級数展開と時間積分

3.1 2 重フーリエ展開

スペクトル法による数値計算を行うために, 各物理量を境界条件を満たす函数系である 2 重フーリエ級数展開する. すなわち,

ζ(x, y, t) = ∑_{k=−K}^{K} ∑_{l=−L}^{L} ζ̃_{kl}(t) exp (i\frac{2πx}{L_{x}}k) exp (i\frac{2πy}{L_{y}} l) ,  (22)
ζ̃_{kl}(t) = \frac{1}{L_{x}L_{y}} ∫_{0}^{L_{x}} ∫_{0}^{L_{y}} ζ(x, y, t) exp (−i\frac{2πk}{L_{x}} x) exp (−i\frac{2πl}{L_{y}}y) dx dy,(23)
plbaro˙ teishiki.tex  平成 20 年 10 月 10 日 (竹広真一)

\underline{2 重周期領域順圧モデル 3 数値計算 : 2 重フーリエ級数展開と時間積分 5}

ζ̃_{kl}(t) は渦度のスペクトルの波数 k, l 成分, K, L はそれぞれ x, y 方向の切断波数である. 流線関数 ψ(x, y, t) のスペクトルも同様に ψ̃_{kl}(t) を定義する.

これらを (17) に代入し, exp (−i\frac{2πk′}{L_{x}} x) exp (−i\frac{2πl′}{L_{y}} y) をかけて全領域で積分す

ることにより, 基底関数の直交関係から各スペクトル成分の式が得られる.

\frac{dζ̃_{kl}}{dt} = − ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl} + (−1)^{p+1}ν_{2p} ˜ [∇^{2p}ζ]_{kl}.

ただし […]_{kl} はスペクトルの k, l 成分であることを表している.

\frac{dζ̃_{kl}}{dt} = − ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl} + (−1)^{p+1}ν_{2p} ˜ [∇^{2p}ζ]_{kl}.

散逸項はラプラシアンが 2 重指数函数系の固有演算子であることから

(−1)^{p+1}ν_{2p} ˜ [∇^{2p}ζ]_{kl} = (−1)^{p+1}ν_{2p} [− (\frac{2πk}{L_{x}} )^{2} − (\frac{2πl}{L_{y}} )^{2}]^{p} · ζ̃_{kl}
= −ν_{2p} [(\frac{2πk}{L_{x}} )^{2} + (\frac{2πl}{L_{y}} )^{2}]^{p} · ζ̃_{kl} = −ν^{∗}(k, l, p) ζ̃_{kl}

となる. ここで, 波数毎の超粘性係数 ν^{∗}(k, l, p) を

ν^{∗}(k, l, p) ≡ ν_{2p} [(\frac{2πk}{L_{x}} )^{2} + (\frac{2πl}{L_{y}} )^{2}]^{p}  (24)

と定義した (ν^{∗}(k, l, p) は 0 以上であることに注意). したがって,

\frac{dζ̃_{kl}}{dt} = − ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl} − ν^{∗}(k, l, p)ζ̃_{kl}.  (25)

3.2 時間積分

以下では時間に関して差分化し, 時間積分を行うための定式化を記す. Δt を時間格子間隔, 時刻 t, t + Δt, t − Δt における ζ̃_{kl} の値をそれぞれ ζ̃_{kl}^{τ}, ζ̃_{kl}^{τ+1}, ζ̃_{kl}^{τ−1} 等と表す.

3.2.1 euler スキーム (1 次精度)

euler スキームの場合は時間変化項をすべて現在時間で評価して,

\frac{ζ̃_{kl}^{τ+1} − ζ̃_{kl}^{τ}}{Δt} = − ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ} − ν^{∗}(k, l, p)ζ̃_{kl}^{τ}.

すなわち,

ζ̃_{kl}^{τ+1} = ζ̃_{kl}^{τ} + (Δt) × \{− ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ} − ν^{∗}(k, l, p)ζ̃_{kl}^{τ}\} .  (26)

3.2.2 Leap frog and Implicit スキーム

散逸項を陰解法 (implicit スキーム), それ以外の項を Leap frog スキームで評価する. まず第一段階として散逸項以外の項に関して積分する.

\frac{ζ̃_{kl}^{τ+1} − ζ̃_{kl}^{τ−1}}{2(Δt)} = − ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ} − ν^{∗}(k, l, p)ζ̃_{kl}^{τ+1}

変形すると

[1 + 2ν^{∗}(k, l, p)(Δt)]ζ̃_{kl}^{τ+1} = ζ̃_{kl}^{τ−1} + 2(Δt) × \{− ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ}\} .

したがって

ζ̃_{kl}^{τ+1} = \frac{ζ̃_{kl}^{τ−1} + 2(Δt) × \{− ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ}\}}{1 + 2ν^{∗}(k, l, p)(Δt)}  (27)

Leap frog スキームを用いているので, 計算モードの増幅を抑制するために Asselin (1972) の時間フィルターをとあるタイムステップ間隔で適用する. フィルターの定義式は

A^{τ∗} = A^{τ} + γ(A^{τ+1} − 2A^{τ} + A^{τ−1})  (28)

ここで A^{τ∗} が修正された τ での A の値, γ はフィルターの係数である.

3.2.3 Adams-Bashforth and Crank-Nicolson スキーム (2 次精度)

散逸項に Clank-Nicolson スキーム, それ以外の項に 2 次の Adams-Bashforth スキームを適用する.

\frac{ζ̃_{kl}^{τ+1} − ζ̃_{kl}^{τ}}{Δt} = \frac{3}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ} − \frac{1}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ−1} − ν^{∗}(k, l, p)\frac{ζ̃_{kl}^{τ+1} + ζ̃_{kl}^{∗τ+1}}{2} ,

plbaro˙ teishiki.tex  平成 20 年 10 月 10 日 (竹広真一)

ただし (\frac{dζ}{dt} )_{nd} は時間変化の非散逸項を表している:

(\frac{dζ}{dt} )_{nd} ≡ −J(ψ, ζ) − β(y)\frac{∂ψ}{∂x} ,  (29)
[ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ} = − ˜ [J(ψ, ζ)]_{kl}^{τ} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}^{τ}  (30)

これを解くと,

[1 + ν^{∗}(k, l, p)(Δt)/2]ζ̃_{kl}^{τ+1} = [1 − ν^{∗}(k, l, p)(Δt)/2]ζ̃_{kl}^{τ}

+Δt × \{\frac{3}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ} − \frac{1}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ−1}\} ,

すなわち,

ζ̃_{kl}^{τ+1} = \frac{[1 − ν^{∗}(k, l, p)(Δt)/2]ζ̃_{kl}^{τ} + Δt × \{\frac{3}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ} − \frac{1}{2} [ ˜ (\frac{dζ}{dt} )_{nd}]_{kl}^{τ−1}\}}{1 + ν^{∗}(k, l, p)(Δt)/2}.

(31)

3.2.4 演算子分割処理法 (粘性項) と 4 次精度 Runge-Kutta スキーム

変数変換することによって粘性項を解析的に評価し, その他の項を高次の数値積分スキームを適用する演算子分割法をとる場合を以下に示す. 渦度のスペクトル成分について

ζ̃_{kl}(t) = ζ̂_{kl}(t)e^{−ν^{∗}(k,l,p)t}  (32)

と変数変換すれば,

\frac{dζ̃_{kl}}{dt} = \frac{dζ̂_{kl}(t)}{dt}e^{−ν^{∗}(k,l,p)t}−ν^{∗}(k, l, p)ζ̂_{kl}(t)e^{ν^{∗}(k,l,p)t} = \frac{dζ̂_{kl}(t)}{dt}e^{−ν^{∗}(k,l,p)t}−ν^{∗}(k, l, p)ζ̃_{kl}(t)

であるから, (36) は

\frac{dζ̂(t)_{kl}}{dt} e^{−ν^{∗}(k,l,p)t} = − ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}.  (33)

\frac{dζ̂(t)_{kl}}{dt} = e^{ν^{∗}(k,l,p)t} [− ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}] .  (34)

\underline{2 重周期領域順圧モデル 3 数値計算 : 2 重フーリエ級数展開と時間積分 8} この左辺の時間積分を Runge-Kutta スキームで計算する. すなわち,

ζ̃_{kl}^{τ+1} = ζ̂_{kl}^{τ+1}e^{−ν^{∗}(k,l,p)Δt},
ζ̂_{kl}^{τ+1} = ζ̃_{kl}^{τ} + (Δt) × (\frac{k_{1}}{6} + \frac{k_{2}}{3} + \frac{k_{3}}{3} + \frac{k_{4}}{6} ) ,
k_{1} = f[ζ^{τ}, ψ^{τ}, 0]_{kl},
k_{2} = f[ζ^{τ+1/2}, ψ^{τ+1/2}, Δt/2]_{kl}
k_{3} = f[ζ^{∗τ+1/2}, ψ^{∗τ+1/2}, Δt/2]_{kl}
k_{4} = f[ζ^{∗τ+1}, ψ^{∗τ+1}, Δt]_{kl}
f[ζ, ψ, t]_{kl} = e^{ν^{∗}(k,l,p)t} [− ˜ [J(ψ, ζ)]_{kl} − ˜ [β(y)\frac{∂ψ}{∂x} ]_{kl}] ,
ζ̃_{kl}^{τ+1/2} = e^{−ν^{∗}(k,l,p)Δt/2}[ζ̃_{kl}^{τ} + k_{1}Δt/2],
ζ̃_{kl}^{∗τ+1/2} = e^{−ν^{∗}(k,l,p)Δt/2}[ζ̃_{kl}^{τ} + k_{2}Δt/2],
ζ̃_{kl}^{∗τ+1} = e^{−ν^{∗}(k,l,p)Δt}[ζ̃_{kl}^{τ} + k_{3}Δt].

この定式化において ζ̃_{kl}^{τ} = ζ̂_{kl}^{τ} であることを用いていることに注意されたい.

3.2.5 演算子分割処理法 (線形項) と 4 次精度 Runge-Kutta スキーム β が一定の場合 (β 平面) には, 粘性項だけでなく β 項も線形となるので演算子分割処理を施すことができる. β が一定の場合には, β 項を渦度を用いて書き直すことができて,

β̃ [\frac{∂ψ}{∂x} ]_{kl} = β ˜ [\frac{∂ψ}{∂x} ]_{kl} = −iβ (\frac{2πk}{L_{x}} ) [(\frac{2πk}{L_{x}} )^{2} + (\frac{2πl}{L_{y}} )^{2}]^{−1} · ζ̃_{kl} = iω_{R}(β, k, l) · ζ̃_{kl}

ここで ω_{R}(β, k, l) はロスビー波の位相速度

ω_{R}(β, k, l) ≡ −β (\frac{2πk}{L_{x}} ) [(\frac{2πk}{L_{x}} )^{2} + (\frac{2πl}{L_{y}} )^{2}]^{−1}  (35)

である. (36) は

\frac{dζ̃_{kl}}{dt} = − ˜ [J(ψ, ζ)]_{kl} − iω_{R}(β, k, l)ζ̃_{kl} − ν^{∗}(k, l, p)ζ̃_{kl}.  (36)

線形項を消去すべく ζ̃_{kl}(t) = Z_{kl}(t) · e^{−[iω_{R}(β,n,m)+ν^{∗}(k,l,p)]t} と変換すれば,

\frac{dζ̃_{kl}}{dt} = \frac{dZ_{kl}}{dt} e^{−[iω_{R}(n,m)+ν^{∗}(n,p)]t} − [iω_{R}(β, k, l) + ν^{∗}(k, l, p)]Z_{kl}(t)e^{−[iω_{R}(β,k,l)+ν^{∗}(k,l,p)]t}
= \frac{dZ_{kl}}{dt} e^{−[iω_{R}(n,m)+ν^{∗}(n,p)]t} − [iω_{R}(β, k, l) + ν^{∗}(k, l, p)]ζ̃_{kl},
plbaro˙ teishiki.tex  平成 20 年 10 月 10 日 (竹広真一)

\underline{2 重周期領域順圧モデル  4 テスト実験  9}

したがって

\frac{dZ_{kl}}{dt} = e^{[iω_{R}(β,k,l)+ν^{∗}(k,l,p)]t}G[ψ, ζ]_{kl}

ここで G[ψ, ζ]_{kl} は非線形項

G[ψ, ζ]_{kl} = − ˜ [J(ψ, ζ)]_{kl}

を表す. ISPACK のスペクトル成分の格納のし方にあわせて配列を実部と虚部にわけて表すと, k > 0 に対して Re[Z_{kl}] = ξ_{kl}ならば Im[Z_{kl}] = ξ_{−k−l} であるから

\frac{dξ_{kl}}{dt} = e^{ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t]Re[G[ψ, ζ]_{kl}]
−e^{ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t]Im[G[ψ, ζ]_{kl}]
= e^{ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t]G[ψ, ζ]_{kl} − e^{ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t]G[ψ, ζ]_{−k−l}, (37)
\frac{dξ_{−k−l}}{dt} = e^{ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t]Im[G[ψ, ζ]_{kl}]
+e^{ν^{∗}(k,l,p)]t} sin[ω_{R}(β, k, l)t]Re[G[ψ, ζ]_{kl}]
= e^{ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t]G[ψ, ζ]_{−k−l} + e^{ν^{∗}(k,l,p)]t} sin[ω_{R}(β, k, l)t]G[ψ, ζ]_{kl} (38)

時間積分のための ξ_{kl}, ξ_{n}^{−m} の各積分ステップでの初期値は ξ_{kl} と Z_{kl} の関係において t = 0 を適用すればよく, ξ_{kl} = ζ̃_{kl}, ξ_{−k−l} = ζ̃_{−k−l} となる. 上の式を用いて数値積分し, あたらしい時間の ξ_{kl}, ξ_{−k−l} が求まれば, 対応して ζ̃_{kl} がつぎのように求められる.

ζ̃_{kl} = Re[Z_{kl}e^{−[iω_{R}(β,k,l)+ν^{∗}(k,l,p)]t}]
= Re[Z_{kl}]e^{−ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t] + Im[Z_{kl}]e^{−ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t]
= ξ_{kl}e^{−ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t] + ξ_{−k−l}e^{−ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t], (39)
ζ̃_{−k−l} = Im[Z_{kl}e^{−[iω_{R}(β,k,l)+ν^{∗}(k,l,p)]t}]
= Im[Z_{kl}]e^{−ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t] − Re[Z_{kl}]e^{−ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t]
= ξ_{−k−l}e^{−ν^{∗}(k,l,p)t} cos[ω_{R}(β, k, l)t] − ξ_{kl}e^{−ν^{∗}(k,l,p)t} sin[ω_{R}(β, k, l)t]. (40)

4 テスト実験

4.1 テスト 1 : 線形項のテスト

単一モード波は非線形項が 0 となるので, 線形方程式の分散関係にしたがう解析解となる (ロスビー波解).

ψ = A_{0} exp(ik̃x + il̃y − ωt),  (41)

を考える. ここで

k̃ = \frac{2πk}{L_{x}} , l̃= \frac{2πl}{L_{y}}

は水平波数であり, k, l は整数である.

ζ = ∇^{2}ψ = −(k̃^{2} + l̃^{2})A_{0} exp(ik̃x + il̃y − ωt) = −(k̃^{2} + l̃^{2})ψ,

であるから, 非線形項が

J(ψ, ζ) = J(ψ, Aψ) = 0,

となる. よって, 残りの線形項に解を代入して分散関係が得られ,

ω(β, k, l) = −\frac{βk̃}{k̃^{2} + l̃^{2}} − iν_{2p}(k̃^{2} + l̃^{2})^{p}

となる. 即ち解は

ψ = A_{0} exp{i[k̃x + l̃y − ω(β, k, l)]t} exp[−ν_{2p}(k̃^{2} + l̃^{2})^{p}t]
= ψ(t = 0) exp[−iω(β, k, l)t] exp[−ν_{2p}(k̃^{2} + l̃^{2})^{p}t],

である. ただし ψ(t = 0) = A_{0} exp[i(k̃x + l̃y)] は初期値流線関数である. 実数関数で書くならば,

ψ = exp[−ν_{2p}(k̃^{2} + l̃^{2})^{p}t] · [Re[ψ(t = 0)] cos ω(β, k, l)t + Im[ψ(t = 0)] sin ω(β, k, l)t]
+i[Im[ψ(t = 0)] cos ω(β, k, l)t − Re[ψ(t = 0)] sin ω(β, k, l)t]}

振幅 A_{0} は任意であるが, 粘性がない場合に振幅が大きくなるとロスビー波解が不安定となることが知られている (e.g. Hoskins and Hollingworth, 1972).

4.2 テスト 2 : 移流項のテスト

まずは β = 0 の場合に, 渦が形を変えずに移流される場合の解を正しく計算できるかチェックする.

そのための解析解として, 「定常な渦を平行流あるいは剛体回転流に重ねた場」を考える. このような場では定常な渦がその形を変えず平行流あるいは剛体回転流によって流されて行くだけである. 以下そのことを示しておく.

今, 定常な渦の流線関数を ψ′ とすると, 定常解であることから

J(ψ′, ∇^{2}ψ′) = 0,

である. 一方, 基本流の流線関数を ψ_{0} とすると, 平行流あるいは剛体回転流なので ∇^{2}ψ_{0} =const. である. したがって, これらを重ねた場の時間発展は

\frac{∂}{∂t} (ψ_{0} + ψ′) = J[ψ_{0} + ψ′, ∇^{2}(ψ_{0} + ψ′)]
= J(ψ_{0}, ∇^{2}ψ_{0}) + J(ψ_{0}, ∇^{2}ψ′) + J(ψ′, ∇^{2}ψ_{0}) + J(psi′, ∇^{2}ψ′),
\frac{∂}{∂t} (ψ_{0} + ψ′) = J(ψ_{0}, ∇^{2}ψ′)

左辺は ψ′ の渦の流れ ψ_{0} による移流を表している. 基本流が平行流あるいは剛体回転ならば ψ_{0} の流れによって移流されても渦の形が変わらないので, 上式がつねに成立ち, 渦が形を変えずに流されていくだけとなる.

局所的な渦で定常なものは軸対称な流線および渦度分布を持つものである. そのような渦の流線と渦度は渦の中心からの距離 r だけで表すことができて

ψ′ = f(r), ∇^{2}ψ′ = g(r).

これに対する移流項のヤコビアンは

J(ψ′, ∇^{2}ψ′) = \frac{∂ψ′}{∂x} \frac{∂(∇^{2}ψ′)}{∂y} − \frac{∂ψ′}{∂x} \frac{∂(∇^{2}ψ′)}{∂y}
= (\frac{∂ r}{∂x} )_{y} \frac{df}{dr} (\frac{∂ r}{∂y} )_{x} \frac{dg}{dr} − (\frac{∂ r}{∂y} )_{y} \frac{df}{dr} (\frac{∂ r}{∂x} )_{y} \frac{dg}{dr} = 0.

したがって定常である.

4.2.1 平行流とガウス型の渦

平行流の基本流に局所的な軸対称のガウス分布の渦を初期値として与える. しかしながら, 今適用しているスペクトル法では, 渦度および流線関数の (k, l) = (0, 0) スペクトル成分は 0 としなければならない. そこで平行流の代わりに領域の半分づつが互いに反平行な流れとなっている状況を与えることにする.

すなわち

U_{x} = \{ \begin{matrix}U_{0} tanh (\frac{y − L_{y}/4}{σ}) & (0 ≥ y ≥ \frac{L_{y}}{2} ) , \\ −U_{0} tanh (\frac{y − 3L_{y}/4}{σ}) & (\frac{L_{y}}{2} ≥ y ≥ L_{y}) , \end{matrix} U_{y} = 0.

U_{0} が流れの振幅, σ が反平行な流れの境界層の厚さを表す. 対応する流線関数は,

ψ_{0}(x, y) = \{ \begin{matrix}−U_{0}σ [log |cosh (\frac{y − L_{y}/4}{σ})| − log |cosh (\frac{L_{y}/4}{σ})|] & (0 ≥ y ≥ \frac{L_{y}}{2} ) , \\ U_{0}σ [log |cosh (\frac{y − 3L_{y}/4}{σ})| − log |cosh (\frac{L_{y}/4}{σ})|] & (\frac{L_{y}}{2} ≥ y ≥ L_{y}) , \end{matrix}

一方, 重ねる局所渦は

ψ′(x, y, 0) = A_{0}[1 − (x − x_{0})^{2} + (y − y_{0})^{2})] exp (−\frac{(x − x_{0})^{2} + (y − y_{0})^{2}}{a^{2}} ) , (42)

A_{0} が渦の振幅, a が渦の特徴的な大きさを表すパラメターである. 第 2 項目は一様渦度成分を 0 とするための定数である. y_{0} を中心付近にとり, 渦の大きさを十分に小さくすれば, 渦は一様流 U_{0} によって流される. このときの解は

ψ′(x, y, 0) = A_{0}[1−(x−x_{0}−U_{0}t)^{2}+(y−y_{0})^{2})] exp (−\frac{(x − x_{0} − U_{0}t)^{2} + (y − y_{0})^{2}}{a^{2}} ) ,

(43)

となる.

注意: 基本流がそもそも不安定 (順圧 or K-H 不安定) な状況ので, 基本場が大きくくずれるまでの間だけ渦が形を変えずに移流される様子が観察できれば OK とする.

4.3 テスト 3 : モドン解

β 一定の無限領域ベータ平面では孤立渦の解析解であるモドンが存在する. その形は渦の中心を原点とする平面極座標 (r, θ) を用いて次のように表される.

ζ(r, θ) = ∇^{2}ψ = −\frac{cq^{2}}{a} sin θ × \{ \begin{matrix}\frac{J_{1}(kr/a)}{J_{1}(k)} & (r < a) \\ \frac{K_{1}(qr/a)}{q_{1}(q)} & (r ≥ a) \end{matrix}  (44)

ここで a c はそれぞれ孤立渦の半径と伝播速度, パラメター q = \sqrt{β/c}, k は分散関係に相当する関係式

−\frac{J_{2}(k)}{kJ_{1}(k)} = \frac{K_{2}(q)}{qK_{1}(q)}  (45)

から定まるパラメターである.

q が実数であることから伝播速度が β と同符号でなければならないことがわかる. β > 0 のとき c > 0 であり, このときは常に東向きに伝播する.

パラメターのとりうる値の一つとして, a = c = β = 1 のとき q = 1, k = 3.9226 の場合がある (McWilliams et al., 1981, 辻村 1993).

\underline{2 重周期領域順圧モデル文献 13}

文献

Hoskins, B. J., Hollingsworth, A., 1972: On the simplest example of the barotropic instability of Rossby wave motion. J. Atmos. Sci., 30, 150-153.

McWilliams, J. C., Flierl, G. R., Larichev, V.D., Reznik, G. M., 1981: Numerical studies of barotropic Modons. Dyn. Atmos. Oceans, 5, 219–238.

辻村豊, 1993: モドンの数値実験, 気象研究ノート, 179, 2.5 章, 215–242.

\underline{2 重周期領域順圧モデル謝辞 14}

謝辞

本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用の実験の一環として

http://www.gfd-dennou.org/arch/spmodel/

において公開されているものである (\begin{matrix}c \\ ◯ \end{matrix}地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2005). 本資源は, 著作者の諸権利に抵触しない (迷惑をかけない) 限りにおいて自由に利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめることをお願いする (無保証無責任原則).

本資源に含まれる元資源提供者 (図等の版元等を含む) からは, 直接的な形での WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるものと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいている. 本資源の利用者には, この点を理解の上, 注意して扱っていただけるようお願いする. 万一, 不都合のある場合には

spmodel@gfd-dennou.org

まで連絡していただければ幸いである.