地衡流乱流のエネルギー散逸とスヶ
–
リング則
九大・応力研
伊賀啓大
(Keita Iga)
Research
Institute
for Applied
Mechanics,
Kyushu Univ.
1
はじめに大気や海洋の大規模な運動は準地衡流渦位方程式にょってよく表されることが知られて
いるが,その方程式に支配される乱流運動は
(
準)
地衡流乱流と呼ばれる. その基本的な 性質は2 次元ナビエーストークス方程式に従うぃゎゆる
2
次元乱流と似ており, エネルギーの逆カスケードなどの現象が見られる
.
しがし, ロスビーの変形半径という特徴的な スケールの存在のために,逆カスヶードの速さなど定量的性質は異なってぃる
.
この準地衡流乱流の性質を詳しく調べたのが
Watanabe
et
al.
(1997, 1998)
である.(
彼らの研究はもともと磁場に対して垂直な面内の準
2
次元静電ポテンシャルを調べるこ とを動機としている.しかし彼らが用いた基礎方程式はチャーニー
-長谷川
- 三間方程式 と呼ばれ,準地衡流渦位方程式と全く同等なものである
.
)
特にWatanabe
et
al. (1998)
は,準地衡流渦位方程式で表される減衰性乱流を調べたものであるが
,
スヶーリング則による議論によってエネルギースペクトルが従う幕則を理論的に導出し
,
数値実験にょって この理論の妥当性を確かめている. しかし,彼らは数値実験の結果を説明する理論を構築する際に全エネルギーの減少のし
かたを表すパラメータを用いてぃるが,
その値の決定には実際に数値計算をして得られ
た結果を使っている. また,理論的にはエネルギースペクトルが,
ロスビー変形半径の逆 数(以下 $\lambda$ と記す) より十分小さな波数領域では波数の-5
乗に, $\lambda$ より十分大きな波数 領域では波数の-3
乗に比例することが予想されてぃるのに,
数値計算の結果では $\lambda$ との 大小に関わらす波数の-6
乗に比例するエネルギースペクトルになってしまってぃる
.
本稿では,準地衡流乱流のスケーリング則につぃて
,
特に, エネルギー散逸に注目して 調べた結果を述べる.2
準地衡流乱流のスケーリング則
大気や海洋の大規模な流れを記述する準地衡流渦位方程式
$\frac{\partial}{\partial t}(\nabla^{2}\psi-\lambda^{2}\psi)+J(\psi, \nabla^{2}\psi)=D$(2.1)
数理解析研究所講究録 1285 巻 2002 年 170-177
を考える. ここで$\psi$ は流線関数, $\lambda$ はロスビー変形半径の逆数, $D$ は非常に高波数の領 域で働く無限小の係数を持った散逸項を表す. まず, この方程式で記述される準地衡流乱 流のスケーリング則を簡単にまとめる. 詳しい導出は
Watanabe
et
al.
$(1997, 1998)$ を参 照のこと. エネルギーカスケード領域におけるエネルギースペクトルのスケーリング則は $E(k)\propto\{$ $\lambda^{2}k^{-\frac{11}{3}}\epsilon^{\frac{2}{3}}$,for
$k\ll\lambda$, $k^{-\frac{5}{3}}\epsilon^{\frac{2}{3}}$,
for
$k\gg\lambda$(2.2)
(ただし, $\epsilon$ はエネルギー輸送率) となり, エンストロフィーカスケード領域におけるエネ ルギースペクトルのスケーリング則は $E(k)\propto\{$$\lambda^{2}k^{-5}\eta^{\frac{2}{3}}$
,
for
$k\ll\lambda$,$k^{-3}\eta^{\frac{2}{3}}$,
for
$k\gg\lambda$(2.3)
(ただし, $\eta$ はエンストロフィー輸送率) である. 減衰性乱流の場合, 渦の特徴的な波数
km。とそれに対するエネルギースペクトルのピーク値 $E_{\max}$ は $k\mathrm{m}\text{。}\ll\lambda$ において,
$k\mathrm{m}\text{。}\propto E^{-\frac{1}{8}}\lambda^{\frac{3}{4}}t^{-\frac{1}{4}}$, $E\mathrm{m}\text{。}\propto E^{\frac{9}{8}}\lambda^{-\frac{3}{4}}t^{\frac{1}{4}}$
(2.4)
のようにそれぞれ$t^{-\frac{1}{4}}$ および$t^{\frac{1}{4}}$ に比例して時間変化する. 但し, $E$ は系の全エネルギー である. この理論は彼らの数値計算の結果をかなりよく説明したが, それでもまだ若干の ずれも見られた. 理想的には全エネルギー $E$ は保存するはずであるが, 数値計算では安 定に時間積分を行なうために, 小さくても有限の大きさを持つ散逸項を入れるので, $E$ はゆっくりと減衰していく. 全エネルギー $E$ が $E\propto t^{-\theta}$
(2.5)
の形で減衰するとすると, km。と Em。の時間変化のしかたは少し修正されて $k\mathrm{m}\text{。}\propto\lambda^{\frac{3}{4}}t^{-\frac{2-\theta}{8}}$ , $E\mathrm{m}\text{。}\propto\lambda^{-\frac{3}{4}}t^{\frac{2-9\theta}{8}}$(2.6)
となる.
Watanabe
et
al.
(1998)
ではこの $\theta$の値を定めるのに彼らの数値計算の結果を用 い全エネルギーの時間変化のグラフから読みとることによって $\theta=0\cdot 05$ としたが, こ の値を用いて表した時間変化
(2.6)
は彼らの数値計算の結果をよく表していたのである.3
エネル$+^{\backslash \backslash }$ー散逸量の見積り
$3\cdot 1$ 高波数域へのエネルギー輸送Watanabe et al. (1998)
で数値計算の結果を用いて経験的に定めたパラメータ $\theta$を理論 的に見積もるための手がかりとして, この系においてエネルギーがどのように散逸するか を考察する.
2
次元乱流や準地衡流乱流においてはエネルギーは全て高波数側から低波 数側に逆カスケードすることを前提にしていた. しかし,Fjortoft (1953)
で議論されて いるように, エネルギーとエンストロフィーの両方が保存するという制約からは, 「たと171
え大部分のエネルギーが逆カスケードしても, 零でないエネルギーが高波数側への輸送さ れなければならない. 」 ということが導かれる. 高波数側に輸送されたエネルギーは散逸 項の影響を大きく受けて散逸されるであろうから, この高波数側に輸送されたエネルギー が全エネルギーの散逸になるとする仮定が (図
1),
これからの議論の出発点である. 図 1: 本理論で仮定するエネルギー輸送・散逸の様子. エネルギーの大部分は元の波数より少し低波数 側にずれるが, わずかな部分は高波数側に輸送されてエネルギーとエンストロフィーの両方が保存さ れる. そして, 高波数側に輸送されたエネルギーが散逸される. エネルギー輸送によってエネルギーもエンストロフィーも変化しない条件より $E$ が満 たす微分方程式は $\frac{dE}{dk_{\max}}\cong\frac{2k_{\max}E}{k_{\eta}^{2}-k_{\max}^{2}}$(3.7)
となるが, この微分方程式は簡単に解けて, 初期におけるエネルギー, エネルギースペク トJレのピーク波数をそれぞれ$E0,$ $kf$ とすると $E \cong E_{0}\frac{k_{\eta}^{2}-k_{f}^{2}}{k_{\eta}^{2}-k_{\max}^{2}}$(3.8)
となる.(3.8)
式は $E$ を km。の関数として表しているが, さらに(2.4)
を用いることに よって, $E$ を時間の関数として表すこともできる.k\eta \gg km。の場合には (3.7)
の分母 の km。を無視して, 全エネルギーの散逸率dE/dt\uparrow ま $\dot{E}_{a}\equiv\frac{dE}{dt}=\frac{dE}{dk_{\max}}\frac{dk_{\max}}{dt}\cong 2k_{\max}k_{\eta}^{-2}E\frac{dk_{\max}}{dt}$(3.9)
となるが,(
$3.9\rangle$ に(2.4)
を代入して $\dot{E}_{a}\cong-\frac{1}{2}C_{ks}^{2}E^{\frac{3}{4}}\lambda^{\frac{3}{2}}k_{\eta}^{-2}t^{-\frac{3}{2}}$(3.10)
となる. ただし, $C_{ks}$ は(2.4)
のkm
。の式の明示されていない比例定数である.
さて, $\theta$ の値という元の問題に戻ろう. このようにして求めた $\dot{E}$ の形から, もはや$E$ は時間の幕乗の形では表されるものではないことがわかる. そのため, $E$ の時間変化を 単純に一つのパラメータ $\theta$ を用いて, $E\propto t^{-\theta}$ と表現するのは不適当である. しかし, 各時刻における $\theta$ に対応する量を一$d\log$$E/d\log t$ で評価することはできる. この値は(3.10)
を用いて, $\theta=-\frac{d1\mathrm{o}\mathrm{g}E}{d1\mathrm{o}\mathrm{g}t}=-\frac{t}{E}\dot{E}_{a}\cong\frac{1}{2}C_{ks}^{2}E^{-\frac{1}{4}}\lambda^{\frac{3}{2}}k_{\eta}^{-2}t^{-\frac{1}{2}}$(3.11)
172
と見積もられる.
この結果には変数$E$ を含んでいる. 散逸の小さい理想的な状態では$E$ はほとんど保存
量であるとして, 初期のエネルギーの値$E_{0}$ を用いればよいであろうが, 実際の数値計算
では, $E$ が元の値$E_{0}$
とは異なってきてしまうことも多
$\text{く}$ ,(3.11)
の正確な見積もりのためには $E$ の値も見積もっておかなければならない. 十分に時間が経つて $k\mathrm{m}\text{。}\ll k_{\eta}$ に
なったとすると, $E$ の値は
(3.8)
から$E \cong E_{0}\frac{k_{\eta}^{2}-k_{f}^{2}}{k_{\eta}^{2}}$
(fo.r
$k_{\max}\ll k_{\eta}$)
(3.12)
と見積もることができる.
3.2
スケーリ $\grave{}$グ則の係数 理論的に見積もった(3.10)
式や(3.11)
式には, 問題設定で与えるパラメータからだけ では決められない量がまだ二つ含まれている. 一つはエネルギー散逸が起こる波数 $k_{\eta}$ で あるが, これは第33
節で具体的な散逸項の形が与えられてから見積もりを行なう.
もう 一つは,Watanabe
etal. (1998)
のスケーリング理論では決められていない比例係数で ある. この節で簡単な仮定に基づいてその値を見積もることを試みる.
$k\mathrm{m}\text{。}\ll\lambda$ の場合の(2.4)
の係数だけでなく, $k\mathrm{m}\text{。}\gg\lambda$ の時の係数についても同時に 考察する.(2.4)
と同様にしてkm
。が大きな極限での $k_{\max}$, Em
。の時間依存性の形も求 めることができ, 両極限での km。’ Em。は km。$\underline{\sim}-\{$ $.C_{ks}E^{-\frac{1}{8}}\lambda^{\frac{3}{4}}t^{-\frac{1}{4}}$, $k\mathrm{m}\text{。}\ll\lambda$, $C_{kl}E^{-\frac{1}{2}}t^{-1}$, $k\mathrm{m}\text{。}\gg\lambda$,
(3.13)
$E\mathrm{m}\text{。}\cong\{$ $C_{Es}E^{\frac{9}{8}}\lambda^{-\frac{3}{4}}t^{\frac{1}{4}}$, $k\mathrm{m}\text{。}\ll\lambda$, $C_{El}E^{\frac{3}{2}}t$, $k\mathrm{m}\text{。}\gg\lambda$(3.14)
と書ける. 但し $C_{ks},$ $C_{kl},$$C_{Es},$$C_{El}$ は比例定数で, この節ではこの係数の見積りを行なう ことになる. まず, km。’Em
。の漸近形は本来$k\mathrm{m}\text{。}\gg\lambda$ (または$k\mathrm{m}\text{。}\ll\lambda$)
の時に成り立つものであるが, これらが $k\mathrm{m}\text{。}>\lambda$ (または $k\mathrm{m}\text{。}<\lambda$
)
でありさえすれば, $\lambda$に近い領域でも 成り立つと仮定する. km。が $\lambda$ になる時刻を$t_{\lambda}$ とする (km。$(t\lambda)=\lambda$
).
$k\mathrm{m}\text{。}<\lambda$ では$k\mathrm{m}\text{。}\cong C_{ks}E^{-\frac{1}{8}}\lambda^{\frac{3}{4}}t^{-\frac{1}{4}}$ が成り立つから, 時刻 $t=t_{\lambda}$ における関係式から $t_{\lambda}\cong C_{ks}^{4}E^{-\frac{1}{2}}\lambda^{-1}$
となるが, 一方, $k\mathrm{m}\text{。}>\lambda$ では$k\mathrm{m}\text{。}\cong \mathrm{Q}_{l}E^{-\frac{1}{2}}t^{-\mathrm{I}}$ だから時刻 $t=t\lambda$ における関係式か
ら $t_{\lambda}\cong C_{kl}E^{-\frac{1}{2}}\lambda^{-1}$ となり, 両者が一致する条件より
$C_{kl}\cong C_{ks}^{4}$
(3.15)
の関係式を得る.
同様にして
(3.14)
式を用いると, 時刻 $t=t_{\lambda}(\cong C_{ks}^{4}E^{-\frac{1}{2}}\lambda^{-1})$ における E。へに関して, $k\mathrm{m}\text{へ}<\lambda$ の関係式から導かれる Em。$(t\lambda)\cong C_{Es}C_{ks}E\lambda^{-1}$ と $k\mathrm{m}\text{。}>\lambda$ の関係式か
ら導かれる E。へ$(t_{\lambda})\cong C_{El}C_{ks}^{4}E\lambda^{-1}$ が一致する条件から次の関係式を得る.
$C_{El}\cong C_{Es}C_{ks}^{-3}$
.
(3.16)
次にエネルギースペクトルと全エネルギーとを関係づける拘束条件につぃて考えよう
.
$k\mathrm{m}\text{。}\gg\lambda$ の場合, $k\mathrm{m}\text{。}<k<\lambda$ {こおけるエネノレギースペクトノレは(2.3)
で示したよう {こ $k^{-5}$ に比例する. 一方,k<km
。におけるエネルギースペクトルの形は指数の正確な値
についてはまだ議論の余地があるが,Watanabe
et al.
(1998)
の結果にょると $k^{4}$ に比例 するようである. k=km。でのピーク値が Em。であるので, このエネノレギースペクト ルの形を積分することにより全エネルギー$E$ は,$E= \int_{0}^{\infty}dkE(k)\cong \mathit{4}^{k_{\mathrm{m}\mathrm{m}}}dkk^{\mathit{4}}k_{\max}^{-4}E_{\max}+\int_{k_{\mathrm{m}\mathrm{r}}}^{\infty}dkk^{-5}k_{\mathrm{m}}^{5}$
axEm
。
$= \frac{9}{20}$km
。
Emax\cong --290
$C_{ks}C_{Es}E$(3.17)
と計算され, $C_{ks}C_{Es} \cong\frac{20}{9}$(3.18)
という関係を得る. 係数を決定するにはもうーっの条件がいる.
ここでは初期における乱流場形成過程に注 目し,「コヒーレントな初期の擾乱のエネルギーが再配分されるのにがかる時間は
,
乱流 の平均速度の流体粒子が波の1
波長を進むのにかがる時間である」
と仮定する. この仮 定により, 流体粒子が1
波長$l_{f}$進むのにかかる時間 $t_{0}$ は$t_{0} \cong\frac{l_{f}}{v_{0}}=\frac{C_{0}}{k_{f}v_{0}}\cong C_{0}k_{f}^{-1}(2E)^{-\frac{1}{2}}$
,
但し $C_{0}=2\pi$(3.19)
と見積もられる. ここで, 初期に与える擾乱の波数$k_{f}$ ま $\lambda$
より十分に大きく, 全エネル
ギー $E$ の大部分は運動エネルギー$v_{0}^{2}/2$ が占めるもの仮定とした.
(3.19)
を$k\mathrm{m}\text{。}\gg\lambda$ で成り立つ $k\mathrm{m}\text{。}\cong C_{kl}E^{-\frac{1}{2}}t^{-1}$
に代入すると $kf\cong C_{kl}E^{-\frac{1}{2}}\cdot C_{0}^{-1}kf2^{\frac{1}{2}}E^{\frac{1}{2}}$ となり,
もうーっ の条件を得る. $C_{kl}\cong 2^{-\frac{1}{2}}\cdot C_{0}$
.
(3.20)
以上で4
つの未知定数に対して4
っの関係式が揃った.(3.15), (3.16), (3.18), (3.20)
から4
つの係数$C_{ks},$ $C_{kl},$$C_{Es},$$C_{El}$ は次のように決定される.$C_{ks}=2^{-\frac{1}{8}}\cdot C_{0}^{\frac{1}{4}}$
,
$C_{kl}=2^{-\frac{1}{2}}\cdot C_{0}$,
$C_{Es}= \frac{20}{9}2^{\frac{1}{8}}\cdot C_{0}^{-\frac{1}{4}}$,
$C_{El}= \frac{20}{9}2^{\frac{1}{2}}\cdot C_{0}^{-1}$.
(3.21)
3.3
超粘性率を用いたエネルギー散逸量の見積り
この系での全エネルギーの散逸量を,
高波数側に輸送されるエネルギーとして見積もっ
てきたが,
実際に数値実験を行なう時にはこの高波数側に輸送されるエネルギーが切断波
数付近に溜ってしまうので,
(2.1)
式の$D$ として, 小さくても有限の大きさの散逸項を 考えることになる.Watanabe
et
al.
(1998)
でも $D$ として小さな超粘性項を用いた$\frac{\partial}{\partial t}(\nabla^{2}\psi-\lambda^{2}\psi)+J(\psi, \nabla^{2}\psi)=\nu(-1)^{p+1}\nabla^{2(p+1)}\psi$
(3.22)
を時間積分している. この基礎方程式
(3.22)
には散逸項が明示されてぃるので, もしエネルギースペクトルの形がわかれば直接エネルギーの散逸量を計算することもできる.
まず, $k_{\eta}$ は全エネルギー散逸と全エンストロフィー散逸の比から $k_{\eta}^{2}\equiv$ $\dot{Q}/\dot{E}$ と計算 されるが, 散逸は主に切断波数 k。に近い高波数側で起き, $k\gg\lambda$ を満たすこの範囲では エネノレギースペクトノレは $E(k)\propto k^{-3}$ になることから, $k_{\eta}^{2} \cong\frac{p-1}{p}k_{c}^{2}$
(3.23)
となる. さて, これまでの結果を用いて(3.11)
で与えられた全エネルギー減衰率の計算を行なうことができる.
(3.11)
の中で用いられる $E$ [ま(3.12)
により $E\cong E_{0}(k_{\eta}^{2}-k_{f}^{2})/k_{\eta}^{2}$ で,$k_{\eta}$ は
(3.23)
G こより $k_{\eta}\cong[(p-\mathfrak{h}/p]^{\frac{1}{2}}$k。で見積もり,Watanabe et al. (1998)
で用$\mathrm{A}\mathrm{a}$ た $E_{0}=0.5,$ $\lambda=40,$ $k_{f}=50,$ $k_{c}=85,$ $p=2$ を代入し, 彼らが評価に用いたと思われる時 亥IJとして $t=15$ をとること{こより, $\theta\cong 0.03$
(3.24)
となる. このようにして理論的に導いた結果はWatanabe
et
al.
(1998)
が数値計算から 見積もった $\theta=0.05$ にかなり近い値を与えた. ところで, エネルギースペクトルの具体的な形が理論的に求まっているので, エネル ギー散逸は超粘性率の値から直接見積もることができる. $k>\lambda$ の範囲では$E(k)$ は$k^{-3}$ に比例して $E(k)=k^{-3}E(\lambda)\lambda^{3}\cong C_{ks}^{5}C_{Es}E^{\frac{1}{2}}\lambda t^{-1}k^{-3}$(3.25)
となるので, エネルギー散逸は主に切断波数 k。に近い高波数側で起こることを考え, そ の領域での具体的な形(3.25)
を用いて,$\ovalbox{\tt\small REJECT}\equiv\frac{dE}{dt}$ $\cong$ $-2 \nu\int_{0}^{k_{\mathrm{c}}}k^{2p}\cdot C_{ks}^{5}C_{Es}E^{\frac{1}{2}}\lambda t^{-1}k^{-3}dk$
$- \frac{1}{p-1}\nu C_{ks}^{5}C_{Es}k_{c}^{2(\mathrm{p}-1)}E^{\frac{1}{2}}\lambda t^{-1}$
(3.26)
を得るのである.
この $\dot{E}_{b}$
と
(3.10)
の E。とは本来同じものになるべきのものであるが, 実際にWatan-abe et al. (1998)
で数値実験の数値を代入すると, $E=0.5\cross[(85^{2}/2)-50^{2}]/(85^{2}/2)\cong$$0.15,$ $\lambda=40,$ $k_{\eta}=85/2^{\frac{1}{2}}\cong 60,$$t=15$ として $\dot{E}_{a}\cong-3\cross 10^{-4}$
なの(こ対して, $E\cong 0.15,$ $\lambda=40,$ $k_{c}=85,$ $p=2,$ $\nu=3.0\cross 10_{\ovalbox{\tt\small REJECT}}^{-8}t=15$ として
$\dot{E}_{b}\cong-2\cross 10^{-3}$ となり, 後者の方が絶対値が
1
桁大きく, 与えた超粘性率が大き過ぎることを示してい る. 逆に,(3.10), (3.26)
の両者が等しくなるという条件から, ちょうどよい超粘性率の 値を定めると, $\nu\cong\frac{p-1}{2}C_{ks}^{-3}C_{Es}^{-1}E^{\frac{1}{4}}\lambda^{\frac{1}{2}}k_{\eta}^{-2}k_{c}^{-2(p-1)}\cdot t^{-\frac{1}{2}}$(3.27)
と計算される.175
4
数値計算
この節では, 数値計算を行なった結果を示す
.
数値積分する式は(3.22)
である. 切断波数$k_{c}=85$, 変換格子
256
$\cross 256$ の擬スペクトル法を用いて数値計算を行なった.
システムの領域は周期境界を持つ $2\pi \mathrm{x}2\pi$ で, ロスビー変形半径の逆数は $\lambda=40$
とする.
散逸は $p=2$ の超粘性を用いる. この時, 超粘性係数の適当な値を
(3.27)
から見積もって $\nu=2\cross 10^{-8}\cdot t^{-\frac{1}{2}}$
とした. 与える初期擾乱はエネルギースペクトルが
$E(k) \propto\frac{k^{30}}{(k+k_{f})^{60}}$
,
$E= \int E(k)dk=0.5$.
(4.28)
となるようにしてあるが, このエネノレギースペクトノレは$k=k_{f}$ [こ鋭1 ピークを持っ. こ こでは $k_{f}=50$ という値を用$\mathrm{V}^{\mathrm{a}}$, 各波数成分には乱数で位相を与えた
.
以上の条件は本 質的に超粘性率が異なるだけで,Watanabe et al.
(1998)
とほぼ同じ条件である. 図2
に流線関数場の時間発展の様子を示す.
パターンが時間とともにゅっくりと大き くなっていき,エネルギーの逆カスケードが起こってぃることが明らがにゎかる
.
この時 のエネルギースペクトルを図3
に示す.Watanabe et al.
(1998)
ではエネルギースペク トルが $k^{-6}$に比例してしまっていたのに対して, $\nu=2\cross 10^{-8}\cdot t^{-\frac{1}{2}}$
の超粘性率を用いた 計算は理論的に予想されるエネルギースペクトル
(2.3)
をよく実現してぃる. さらに, こ の時間発展に伴う諸量の振舞いを見てみよう.
図2: 流線関数の時間変化の様子. 時間とともに 図3: $p=2,$$\nu=2\mathrm{x}10^{-8}\cdot t^{-:}$ のを用いた数値 パターンのスケールがゆっくりと大きくなって 計算で得られるエネルギースペクトル. $k<\lambda$ レ$\mathrm{a}$る. では$k^{-5}$ [ニ, $k>\lambda$ では$k^{-3}$ [こ比例してぃる. 図4
は km。の時間変化を示してある. 理論で予想されるように $t^{-\frac{1}{4}}$ に比例して時間変 化しているが, それだけでなく, 実線で示されている(3.13)
の理論結果は絶対値として もかなりいい見積もりとなっていることがわかる. また, 図5
は全エネルギーの時間変 化を示す. 全エネルギーはゆっくりと減衰しているが, 時刻 $t=10\sim 20$ 付近では $t^{-0.05}$176
に比例しており, この指数の値は
