2. モデルと数値計算法


次の渦度方程式に支配される回転球面上の2次元非圧縮流れを考える.


ここで, ψ(λ, μ, t): 流れ関数, ω(λ, μ, t): 渦度の鉛直成分(≡ Δψ), λ: 経度, μ: サイン緯度, t: 時刻, Ω: 球自転角速度, ν2p: 超粘性係数, p: 超粘性の次数, である. また, Δは 水平ラプラシアンで以下のように定義される.


なお,

従属変数 ψ を球面調和関数で以下のように展開する:


ここで, 上の展開に基づいて, スペクトルガラーキン法により, 渦度方程式から 展開係数 ψnm(t) に対する連立常微分方程式系を導く. 非線型項は実空間で計算し, フーリエ-ルジャンドル変換で 波数空間に戻す変換法を用いる. エイリアシングを除 くために, N=682 で切断した場合には, 全球面領域を 2048(経度方向)×1024(緯度方向)の変換格子で, N=341 で切断した場合には, 1024×512の変換格子で覆う.

減衰性乱流の実験では, ある初期乱流場を与えて時間発展をみることになる. エネルギー密度は,


と表わされる. 初期のエネルギースペクトルが Cho and Polvani(1996)の用 いた分布関数,


になるように拘束して, 各帯状波数 m 成分の振幅と位相を乱数で与え る. ここで, 係数 A は無次元化した全運動エネルギー: Etotal ≡ Σn=2N E(n,t=0) が1となるように定める. n0 は中心波数 を与えるパラメータで, 10, 50 または 100 とする. また, γ はスペ クトル幅をきめるパラメータで, 大きい値ほど初期エネルギーを狭い波数帯に 限ることになる. 本研究では, 1000 とした.

自転角速度 Ω を実験パラメータとして, 0, 25, 50, 100, 200, 400の6つの値で与え る. Ω≠0の場合には, 惑星は単位時間あたり Ω/(2π) 回 転することになる.

ここで, 渦度方程式の非線型項が線型の自転効果項と同程度の大きさとなる波数 (Rhinesスケールに対応する波数): nβ≡√(πΩ/(4U))= √(πΩ/(4√(2×Etotal)) ) =√(πΩ/(4√2)), が重要な指標となる. 実際, 減衰性乱流の時間発展の様相は, nβと初期エネルギー ピークの波数n0との大小関係に大きく依存している. nβ > n0 ならば, 初期から線型項が非線型項より卓越するので, 必ずしも乱流理論で考 えたエネルギー逆カスケードの状況が出現するわけではない.

数値粘性は, p=8 とし, ν2p = 1×10-43 (N = 682のとき); 1×10-38 (N = 341のとき), で与える. 時間積分には4次精度のルンゲ-クッタ法を用い, 時間刻みは Δt = 5×10-4 (N = 682のとき), 1×10-3 (N = 341のとき); として, t=5 まで時間発展を計算する.

なお, 数値計算には, ISPACK(石岡, 1999) 中の効率的変換コードを用いている.


<<index <previous next>