修士学位論文

高精度有限体積法による非静力学惑星大気循環モデルの開発

専攻名地球惑星科学専攻
学籍番号 114s407s
氏名河合佑太

神戸大学大学院理学研究科博士課程前期課程

平成 24 年 3 月 21 日

要旨

本論文では, 空間高精度で計算効率が良く, 保存変数の保存性を満たす, 高精度有限体積法による全球完全圧縮非静力学モデルを提案する. まず, 全球格子点モデルの困難の一つである極の特異性は, 方程式系を立方球面座標系で記述することによって回避している. また, 長時間積分を行う全球大気モデルにおいては系の全エネルギーや質量の保存が重要であるため, 支配方程式系として保存形式の密度, 運動量, 全エネルギーの時間発展式で構成される完全圧縮非静力学方程式系を採用し, 空間離散化には有限体積法を適用する. この支配方程式系の数値解法には, 数値流体分野における低マッハ数圧縮性流れに対する比較的新しい数値解法を応用することができる. 本研究ではそういった手法を積極的に取り入れることにする. 具体的には, 数値流束関数として, 低マッハ数圧縮流れに対しても精度良く安定である AUSM^{+}-up スキームを用いる. 再構築には CENO(Central Essentially Non-Oscillatory) 法を用いることで, 数値精度と計算効率のバランスが取れた任意精度の空間離散化を実現する. また時間積分に関しては, 時空間精度を維持しつつ, 狭い鉛直格子幅による時間刻み幅の厳しい制約を取り除くために, 演算子分割 Runge-Kutta-Rosenbrock スキームを用いて, 水平方向の運動を陽的, 鉛直方向を運動を陰的に取り扱う. 次に, 開発した数値モデルの性能を検証するために, 非静力学コアに対する標準テストケースを三種類行った. テスト計算の結果によって, 本モデルに適用された数値手法は系の質量と全エネルギーを丸め誤差の範囲で保存することが示された. また, 高精度な再構築と低マッハ数流れのために改良された数値流束関数を用いることにより, 流れの場の中に存在する急勾配を維持できることが確認された. 最後に, 次世代非静力学コアに対する数値解法として数値流体分野の手法の適用性や, 本モデルによる数値シミュレーションの可能性について述べる.

\underline{1}

目次

第 1 章序論 3

第 2 章立方球面格子 7

第 3 章非静力学コアの定式化 11 3.1 基礎方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 立方球面座標系の導入 . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 浅い大気近似を施した場合の定式化 . . . . . . . . . . . . . . . . . . 16

第 4 章数値解法 17 4.1 有限体積法による空間離散化 . . . . . . . . . . . . . . . . . . . . . 18 4.2 再構築 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4.2.1 Central ENO 法による水平方向の再構築 . . . . . . . . . . . 21 4.2.2 鉛直方向の再構築 . . . . . . . . . . . . . . . . . . . . . . . . 24

4.3 検査セル境界における数値流束の評価 . . . . . . . . . . . . . . . . 25 4.3.1 直交化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3.2 AUSM^{+}-up スキーム . . . . . . . . . . . . . . . . . . . . . . 27 4.3.3 鉛直方向の数値流束のための AUSM^{+}-up スキームの凖線形化 29

4.4 時間積分法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.5 パッチ境界の取り扱い . . . . . . . . . . . . . . . . . . . . . . . . . 31

第 5 章数値実験 34 5.1 テストケース 1: 音波 . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.2 テストケース 2: 内部重力波 . . . . . . . . . . . . . . . . . . . . . . 39 5.3 テストケース 3: 傾圧不安定 . . . . . . . . . . . . . . . . . . . . . . 42

第 6 章結論 48

謝辞 50

付録 A 球面立方座標系の幾何的な性質 51 A.1 計量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 A.2 球面立方座標系における微分演算子 . . . . . . . . . . . . . . . . . . 52 A.3 地理座標系と球面立方座標系間の座標変換 . . . . . . . . . . . . . . 53

\overline{main.tex 2013/04/02}

\underline{2}

A.3.1 赤道パッチにおける座標変換 . . . . . . . . . . . . . . . . . 54 A.3.2 極パッチにおける座標変換 . . . . . . . . . . . . . . . . . . . 54

付録 B 幾何モーメントの計算法について 55

\overline{main.tex 2013/04/02(河合佑太)}

第1章序論

全球的な数値天気予報や気候予測において, 大気大循環モデル (Atmospheric General Circulation Model: AGCM) は重要な道具である. 継続的な計算機技術の発展は, 全球規模の気象シミュレーションにおいても水平解像度が数 km 以下の高解像度計算を可能にしつつある. 伝統的な AGCM は静水圧平衡近似を用いるプリミティブ方程式を使用してきたが, 数 km 以下の水平スケールの運動に対しては静水圧平衡の仮定が成り立たないため, 近年開発される AGCM において, 静水圧近似を仮定しない非静力学方程式が使われる傾向にある. 研究領域における非静力学モデルの歴史は 30 年ほどあるが, 現業の数値予報分野で非静力学モデルが使われるようになったのは 10 年前ほどからであり, さらに領域モデルに限られていた. このような全球非静力学モデルとして, 日本では, 東京大学気候システム研究センター (CCSR) と地球環境フロンティア研究センター (FRCGC) で共同開発された, 非静力学正二十面格子大気モデル NICAM(Tomita and Satoh, 2004) がある. 他の正二十面格子を使った全球非静力学モデルには, アメリカ海洋大気庁 (NOAA) による, 非静力学正二十面格子モデル NIM や海洋 - 陸 - 大気モデル OLAM などが挙げられる. また, 立方球面格子を用いた全球モデルでは, 地球流体力学研究所 (GFDL) とアメリカ航空宇宙局 (NASA) のゴダード宇宙飛行センター (GSFC) による, GEOS FV CUBED(Putman and Lin, 2008) がある. 一方, ヨーロッパ中期予報センター (ECMWF) はセミラグランジュ・スペクトル変換法による全球非静力学モデル IFS(Wedi et al., 2009) を開発している.

伝統的にほとんどの AGCM の数値解法には, 格子点法に比べ高い数値精度を持つ, 球面調和関数展開に基づくスペクトル変換法が採用されてきた. しかし, 高性能計算 (High Performance Computing: HPC) で実現されるような高解像度計算では, ルジャンドル変換の O(n^{3}) (n は緯度経度格子における一次元的な格子数) の計算量や大規模並列計算機におけるノード間データ転送がネックになることが指摘されている. 一方, アメリカ国立大気研究センター (NCAR) の WARF(Skamarock and Klemp, 2008), イギリス気象庁 (UK Met Office) の UM(Davies et al., 2005) など現業の多くの領域非静力学モデルは, 空間離散化に格子点法の一種である有限差分法を用いる. 有限差分法は効率的であるが, 位相誤差, 領域形状に対する柔軟性, 並列計算におけるスケーラビリティの問題が存在する.

これらの欠点を克服するために, 近年開発されている領域および全球非静力学モデルでは, 有限体積法 (Finite Volume Method: FVM), 有限要素法 (Finite Element Method: FEM), スペクトル要素法 (Spectral Element Method: SEM), 不連続 Galerkin 法 (Discontinuous Galerkin Method: DGM) などがしばしば用いられる. これらはもともと数値流体力学 (Computational Fluid Dynamics: CFD) 分野で活発に研究されてきたものであり, 気象あるいは地球流体力学分野での応用は 10 年前ほどからである. 気象分野での FVM の応用例は, 正二十面格子全球非静力学モデル NICAM などが挙げられる. 全球モデルの難しさの一つは極の特異性であり, 緯度経度格子の極における格子線の収束は時間積分の時間スッテプの大きさに厳しい制約を与える. NICAM では六角型の非構造格子を水平格子に採用し, 有限体積法を用いて空間離散化することで, 極の特異性を回避している.

次世代の非静力学モデル開発における数値解法と関連した問題は,

1. 空間高精度化,

2. 計算領域形状に対する柔軟性,

3. 優れた分散特性,

4. 並列計算機におけるスケーラビリティ

などである. 1 の空間高精度化は, 計算コストを抑えながら, 大気の流れの詳細な構造を解像するためには不可欠である. また, 2 は急峻な地形を始めとする複雑な境界の取り扱いが必要な場合に重要である. 1,2 と関連することとしては, 静的あるいは動的な格子の局所的な細分化が挙げられる. 3 の分散特性については, 低解像度計算においても大気波動の分散関係を正しく再現できる性質が望まれる. 4 と関係することとしては, (LINPACK ベンチマークによる) スーパーコンピュータ Top500 が反映しているように, 1990 年代後半からの HPC の傾向として, ベクトル型よりもスカラ型のスーパーコンピュータ (例えばクラスタ) の開発に重点が置かれている事実である. ベクトル型と比較したスカラ型のスーパーコンピュータの特徴は, 安価な汎用マルチコアプロセッサを大量に用いて並列計算を行うところである. 例えば, 2012 年に完成した理化学研究所の京 (けい) の総コア数は 50 万個以上に及ぶ. このような, 膨大な数のプロセッサコアを持つ計算機の利点を生かすには, 一プロセッサに対する演算量が多く, ノード間でのデータ通信の少ない, より局所的な数値解法が望ましいであろう. したがって, このような HPC の動向も考慮すると, 新しく作られる全球非静力学モデルに採用される数値解法には, (メモリ空間上で) 局所的なデータに対して多くの演算を行うことにより高精度化する, コンパクトステンシルを用いる高精度 FVM や FEM, SEM, DGM といった手法がやはり適していると思われる. また, これらの手法は 1, 2, 3 の問題も得意とするところである.

前述したように, 地球流体力学的問題に対して FVM, FEM, SEM, DGM を適用した研究は, 2000 年代に入ってから見られるが, 当初は全球浅水モデルであったり, 全球非静力学モデルであっても古典的な FVM (空間 2 次精度) に限られていた. 特に, 領域または全球非静力学モデルに対する適用は数年前からであり, まだ広く研究されていない. これらの研究例には, 高精度 FVM による全球非静力学モデル Mcore(Ullrich and Jablonowski, 2012b), SEM または DGM による領域非静力学モデル (Giraldo and Restelli, 2008) などが挙げられる. このため, 本研究では CENO 法を用いた高精度 FVM による全球非静力学モデル DECAMO(DEnnou Cubed sphere grid point Atmospheric MOdel) を提案する.

ここからは, 本モデルの特徴を簡単に述べることにする. まず, 全球モデルに伴う極問題を解決するために, 格子系には立方球面格子を用いる. なお, Taylor et al. (2008) や Putman and Lin (2007) によって, 立方球面格子を用いた数値モデルは, 超並列計算機においてスケーラビリティが保たれることが示されている. また, 本モデルの支配方程式は, 保存形式の密度, 運動量, 全エネルギーの時間発展式で構成される完全圧縮方程式系であり, これらの空間離散化には有限体積法を用いる. これにより, 本質的には系の全質量, 全運動量, 全エネルギーは保存される. ただし, 非線形方程式の反復計算の残差や浮動小数演算の丸め誤差のために, 非常に小さな変動は存在する.

有限体積法の高精度化において, セル界面での場の値を正確に与えるための再構築 (reconstruction) の方法と, セル界面における数値流束の評価方法の選択が重要となる. 再構築において, 検査体積内のサブグリッドスケールの分布は, 近傍にある複数の検査体積 (ステンシル) の体積積分平均値の情報を用いて生成された, 補間多項式で近似される. 並列計算の効率性の観点から, 再構築においてできるだけコンパクトなステンシルが設定されることが望ましい. CFD 分野においてこの再構築の技術は古くからよく研究されており, 近年では WENO 法 (Weighted Essentially Non-Oscillatory Method)(Liu et al., 1994) などが有名である. これらの手法には, 単調性を保証する工夫と一緒になっている場合が多い. 本研究では, WENO 法よりもより効率的で安定な Central ENO 法 (Ivan, 2010; Ivan et al., 2012) を用いる. これらの手法は, 原理的には, より高次の補間多項式を生成することによって任意空間精度の有限体積法を実現することができる.

次に, 検査体積の界面における数値流束の評価方法の選択について述べる. 数値流束の評価する方法としては, 近似リーマン解法 (例えば, Roe スキーム (Roe, 1981)) や比較的新しいものには風上分離法がある. 伝統的にはこれらの方法は, 大気のような低マッハ数の流れよりもむしろ, 遷音速や超音速のレジームにある流れを模倣するために設計されている. そのため, 低マッハ数の流れに対して近似リーマン解法を用いた場合に方程式系が stiff になるため, 陰解法における反復計算が収束しなかったり, また陽解法でも数値精度が悪いなどの問題がある. しかし, 1990 年代後半から, 近似リーマン解法や風上分離法は低マッハ数流れに対しても適用できるように拡張されており, 前処理を用いた近似リーマン解法 (Weiss and Smith., 1995) や低マッハ数流れ用の風上分離法 (Liou, 2006) などが挙げられる. 本研究では, Ullrich and Jablonowski (2012b) と同様に, 計算効率と数値精度の観点から, Liou (2006) で提案された AUSM^{+}-up スキームを数値流束の評価に用いる.

全球モデルでは, 依然として鉛直格子幅は水平格子幅に比べ, 10^{−2} ∼ 10^{−3} 倍ほど小さい. 特に, 速い位相速度で鉛直方向に伝播する音波や重力波を含む非静力学モデルでは, この相対的に小さな鉛直方向の格子幅により CFL 条件 (Courant-Friedrichs-Lewy Condition: CFL condtion) を満たす時間ステップ幅は, 総観規模の運動がもつ時間スケールに比べて非常に小さくなる. したがって, 多くの全球モデルでなされるように, 本研究でも鉛直方向の運動に対する時間離散化は陰的に行う. 時空間精度を保持しながら, 水平方向の運動に対する陽的な積分と鉛直方向の運動に対する陰的な積分を組み合わすために, Ullrich and Jablonowski (2012b) で採用された演算子分割 Runge-Kutta-Rosenbrock スキームを用いる. 今回示す数値解法では, 水平方向に関しては空間精度は 3 次精度, 鉛直方向に関しては空間 2 次精度の再構築を行う. 時間精度に関しては, 水平方向の陽的な時間積分には, 4 次の Runge-Kutta 法を用いる. 演算子分割 Runge-Kutta-Rosenbrock スキームは, 時間積分一ステップあたりに陽解法と陰解法を複数回組み合わせることで, 次元分割していても最高 3 次精度の時間積分を可能にする. 今回は Ullrich and Jablonowski (2012b) に習って, 全体としての空間・時間精度が 2 次精度となるような方法を適用した. そのため, 結果的に全体としては水平空間精度の 3 次収束は保証されないが, Ullrich and Jablonowski (2012b) が示すように, 惑星波から中速度の重力波の振動数を持つ水平伝播する波の分散関係の表現が改善される. なお, 鉛直方向に対してもより高次の再構築を用い, かつ全体で 2 次精度以上を保つような演算子分割 Runge-Kutta-Rosenbrock スキームを適用すれば, 時空間 2 次精度以上を実現できる.

開発したモデルの性能を検証するために, 非静力学モデルの力学コアに対する標準テストケースである, 音波や内部重力波のテストケース, 傾圧不安定のテストケースを行う. これらの数値計算の結果は, 高精度 FVM と低マッハ数流れ向けの数値流束の評価方法のおかげで, 古典的な FVM による非静力学モデルに比べて場の急勾配をより捕捉できる特性を示している.

最後に, この修士論文の構成について説明する. 第 2 章は, 立方球面格子を導入する. 第 3 章では, 立方球面格子座標系において本モデルの支配方程式の定式化を行う. 第 4 章において, 支配方程式を解くために用いられる数値解法について述べる. 第 5 章では, 開発したモデルの数値精度や特性を検証するために行った, 非静力学コアのための標準テストケースの結果を示す. 第 6 章では, 結論および今後の開発方針や本モデルを用いた数値シミュレーションの可能性を議論する.

第2章立方球面格子

本章では, 本モデルにおいて格子系として採用される立方球面格子について述べる.

立方球面格子は, Sadourny (1972) や Ronchi et al. (1996) によって考案された格子である. 立方球面格子では, 立方体を内接する球面に立方体の各辺を投影することによって, 球面を 6 個の同一の矩形領域 (本論文ではこれをパッチと呼ぶ) に分割し, それぞれのパッチには同一の局所座標系を導入する. そのため, 6 個のパッチ間で支配方程式系の見かけは変わらない. 図 2.1 は, この投影によって生成された立方球面格子を表す模式図である. 立方球面格子は, 全球モデルで長年採用されてきた緯度経度格子が持ついくつかの欠点を回避できるため, 正二十面格子や陰陽格子とともに近年注目されている. 緯度経度格子の欠点は経線が両極で集中するこであり, そのためこの極における特異点に対して特別な取り扱いが必要となる. 一方, 立方球面格子では取り扱いの難しい両極の特異点は存在せず, 代わりに, それよりも幾分取り扱い易い 8 個の特異点がもとの立方体の三面の交点で現れる.

図 2.1: 左図は心射投影による立方球面格子の球面上の分布を表している. 一方, 右図は立方球面格子の展開図を示している.

立方球面格子の各パッチに導入される局所座標系には大きく分けて 2 種類が存在し, Sadourny (1972) によって考案された心射投影 (gnomonic projection) によるものと, Ronchi et al. (1996) によって考案された正角投影 (conformal projection) によるものがある. 心射投影には, 解像度が高くなるにつれて検査セルの最小最大面積比が 1/\sqrt{2} に収束する良い特性があり, 正角投影に比べて球面上でより一様な格子を生成することができる. しかし, 心射投影で導入される局所座標系はパッチ境界で不連続であるといった欠点がある. また, その局所座標系は非直交であり, そのため支配方程式は直交座標系において記述した場合に比べて複雑となる. 一方で, 正角投影による局所座標系はパッチ間で連続的に接続し, また局所的に直交するので支配方程式はより簡単に記述される. 正角投影の欠点は, 局所座標系がパッチの頂点の周辺で格子線が集中する点と, 解析的なメトリック項が存在しない点である. なお, Adcroft et al. (2004); Putman and Lin (2007) は, 正角投影による立方球面格子の非一様性を改善する方法を説明している.

本研究では, 再構築には, 近接する検査体積セルの積分平均値を用いて補間多項式を生成する方法を用いるため, 心射投影におけるパッチ間の格子線の不連続はあまり問題にはならない. 一方で, 正角投影における格子の非一様性は時間積分の時間刻み幅に制約を与えるので, 本研究では心射投影による立方球面格子を採用する.

心射投影による立方球面座標系は, 等角座標を用いる場合, 一般的に成分ベクトル (r, α, β, n_{p}) を使って与えられる. ここで, r は球の中心から放射方向の距離, α, β ∈ [−\frac{π}{4}, \frac{π}{4}] は各パッチに対する水平方向の座標であり, 心射等角投影における各方向の中心角である. また, n_{p} ∈ {1, 2, 3, 4, 5, 6} はパッチのインデックスを表している. 図 2.2 は, パッチ 2 における心射投影の関係を表している. この図において, 点 P は球面上の点であり, P′ は点 P と対応する立方体の面上の点である. さらに, 等角座標と近い関係にある心射座標を用いた場合の座標系, (r, X, Y, n_{p}) を導入する. ここで,

X = tan α Y = tan β.  (2.1)

心射座標系は, 心射投影において直線が球面上で必ず大円の弧を形成するという性質を持つ. また, 心射座標系の座標変数を用いた方が, 幾何的な項をより簡単に書ける利点がある. このため, Urllich et al(2012) と同様に, 二種類の座標系 (r, α, β, n_{p}), (r, X, Y, n_{p}) を相互に用いることにする. また, 立方球面格子と関係するメトリックの計算において, 次の定義,

δ = \sqrt{1 + X^{2} + Y^{2}}  (2.2)

をよく用いる.

等角立方球面座標において, (α, β) ∈ [−\frac{π}{4}, \frac{π}{4}] 平面を各方向に一様に N_{c} 分割して生成される格子を考える. このとき, 各パッチの球面上に存在する各方向の格子線は一定の角度 Δ の円弧として分割される. この (α, β) 平面の分割数 N_{c} と Δ の関係は,

Δ = \frac{π}{2N_{c}}  (2.3)

\underline{9}

{図 2.2: パッチ 2 に対する心射投影の関係. 図中の α, β は心射投影における中心角, λ, θ は経度と \atop 緯度である.} \overline{main.tex 2013/04/02(河合佑太)}

である. 歴史的に, 立方球面格子の解像度は c⟨N_{c}⟩ と簡潔に記述する. ここで, N_{c} は各パネルにおける各水平方向の格子セル数と一致する. したがって, ある高度 r における全格子セル数は, 6 × N_{c}^{2} である. なお, 放射方向の格子幅は水平格子に依存せず, 任意に選択することができる.

心射投影による立方球面格子の計量や座標変換などの詳細に関しては, 付録 A の「球面立方格子の幾何的な性質」において説明する.

第3章非静力学コアの定式化

本章では, 非静力学モデルの定式化を行う. 初めに, 座標系に依存しないベクトル形式で基礎方程式系を記述した後に, 立方球面座標系を導入し, 各方向成分に対する方程式に書き下す. DECAMO のオプションとして, 支配方程式系に浅い大気近似を適用した場合の数値シミュレーションも可能にするために, 浅い大気近似を適用した支配方程式系も同時に導く.

3.1 基礎方程式系

DECAMO では基礎方程式系として, 密度, 運動量, 全エネルギー密度といった保存変数を使った次の非静力学方程式系,

\frac{∂ρ}{∂t} + ∇ · (ρu) = 0,  (3.1)
\frac{∂(ρu)}{∂t} + ∇ · (ρu ⊗ u + pI) = −ρ∇Φ − 2Ω × (ρu),  (3.2)
\frac{∂(ρe_{t})}{∂t} + ∇ · (ρu (e_{t} + p)) = 0  (3.3)

を用いる. ここで, ρ は密度, u は 3 次元速度ベクトル, e_{t} は単位体積・単位質量当たりの全エネルギー, p は圧力, Φ は重力ポテンシャル, Ω は惑星自転角速度ベクトル, ⊗ はテンソル積, I はランク 3 の単位行列である. 全エネルギー密度 e_{t} は, 単位体積・単位質量当たりの内部エネルギー, 運動エネルギー, 位置エネルギーの和で表される. すなわち, 単位体積・単位質量当たりの内部エネルギーを e とすれば,

e_{t} = e + \frac{1}{2}u · u + Φ.  (3.4)

また, 運動量方程式 (3.2) に現れる圧力は, 理想気体を仮定すれば, 他の予報変数を使って次のように書ける.

p = (γ − 1) (ρe_{t} − \frac{1}{2}ρu · u − ρΦ) .  (3.5)

ただし, 理想気体の状態方程式

p = ρR_{d}T  (3.6)

により, 内部エネルギー密度 ρe が

ρe = ρC_{v}T = \frac{p}{γ − 1}  (3.7)

と書けることを用いた. ここで, T は温度, R_{d}, C_{v}, C_{p} はそれぞれ乾燥大気に対する気体定数, 定積比熱, 定圧比熱であり, γ(= C_{p}/C_{v}) は比熱比である. (3.1)-(3.3) に加えて, 診断方程式 (3.5) を用いることで, この方程式系は閉じられる.

ここで, 完全圧縮系を用いる非静力学モデルにおける予報変数の組み合わせについて述べておく. 本モデルでは, 予報変数の組み合わせとして, 密度, 運動量, 全エネルギーを用いている. この組み合わせは, CFD 分野の圧縮性流れの数値解法で用いられる支配方程式系としては自然な選択である. しかしながら, 現業の領域完全圧縮系非静力学モデルにおいて, これらの予報変数の組み合わせが使われることは少ないように思われる. 例えばアメリカ大気研究局 (NCAR) を中心にして開発された気象解析予報モデル WRF(Skamarock and Klemp, 2008) に見られるように, 多くは全エネルギーの代わりにエントロピー密度 (温位に密度を掛けたもので, 理想気体において圧力と等価) が用いられている. これは既存のサブグリッドスケールの湿潤物理過程のパラメタリゼーションが温位を利用するためであり, 全エネルギーを予報変数として用いた場合には, 温位を診断する余分なステップが必要となる. なお, 全エネルギーの代わりにエントロピー密度を予報変数とするモデルでは, 系の全エネルギーは本質的には保存しない.

近年, 全球非静力学モデルが開発され, 気候シミュレーションなどにおいて長時間積分がなされるが, このとき系のエネルギー収支の観点から, 系の全エネルギーの保存性がしばしば重要視される. そのため, Satoh (2002) (NICAM の非静力学コアの定式化) では, 全エネルギーの式を用いて全エネルギーを時間発展し, 内部エネルギー密度 (理想気体において圧力と等価) は診断により決定する, 完全圧縮系非静力学モデルの保存スキームを提案している. また, Giraldo and Restelli (2008) は, 完全圧縮系非静力学モデルにおける予報変数の選択として, 全エネルギーを用いることを推奨している. その理由として, 彼らは, 非粘性流体のレジームと粘性流体のレジーム両方で物理的に一貫性がある点, 粘性流体おける位置エネルギーおよび運動エネルギーの内部エネルギーへの散逸を自然に考慮できる点を挙げている. したがって, 本研究では, 予報変数の組み合わせとして密度, 運動量, 全エネルギーを用いることにする.

ただし, 大気現象のような低マッハ数圧縮流れに対し, 予報変数として密度, 運動量, 全エネルギーを選択し, 伝統的な圧縮性流れの数値解法をそのまま適用することはできない. なぜならば, 方程式系が stiff になるため, 反復解法において残差が減少せず解がなかなか収束しない問題や数値解の精度が悪化する問題が, CFD 分野ではよく知られているからである. これらを回避するためには, 低マッハ数流れに拡張された圧縮性流れの数値解法を用いらなければならない. その数値解法の詳細は, 第 4 章で述べる.

次に, 状態変数を基本場と擾乱場に分離する. 本論文では, 基本場として静水圧平衡の状態を選択することにし, 状態変数を次のように分離する.

\begin{aligned}ρ(x, t) & = \overline{ρ}^{h}(x) + ρ′ (x, t), \\ p(x, t) & = \overline{p}^{h}(x) + p′ (x, t), \\ ρe_{t}(x, t) & = \overline{(ρe_{t})}^{h}(x) + (ρe_{t})′ (x, t). \end{aligned}  (3.8)

ここで, x = (r, α, β, n_{p}) である. また, 添字\overline{( )}^{h} が静水圧平衡による寄与示し, それらの変数は静水圧平衡,

\frac{∂\overline{p}^{h}}{∂r} = −\overline{ρ}^{h}g (\frac{a}{r})^{2}  (3.9)

を満たしている. ここで, a は惑星の平均半径, g = GM/a^{2} は重力加速度である. ただし, G は万有引力定数, M は惑星の質量である. このとき, 全エネルギー密度の静水圧平衡成分は,

\overline{(ρe_{t})}^{h} = \frac{\overline{p}^{h}}{γ − 1} + \overline{ρ}^{h}ga (1 − \frac{a}{r})  (3.10)

である. このように基本場として静水圧平衡を指定することによって, 偽の鉛直運動量の生成の要因となり得る, 大気の静水圧状態を近似した場合の誤差を取り除くことができる.

3.2 立方球面座標系の導入

前節で示した基礎方程式系に対して, 第 2 章の立方球面座標系を導入する.

立方球面座標系の共変基底を {e_{α}, e_{β}}, 放射方向の基底を e_{r} とすると, 3 次元速度ベクトル u は, 反変成分 (u^{α}, u^{β}, u^{r}) を使って,

u = u^{α}e_{α} + u^{β}e_{β} + u^{r}e_{r}  (3.11)

と書ける. 前節の運動量方程式 (3.2) はベクトル形式で書かれているが, これを各方向の運動量成分 (ρu^{α}, ρu^{β}, ρu^{r}) に対する方程式に分解する. 残りの方程式 (3.1), (3.3) は, スカラー方程式であるのでそのまま用いる.

(3.1)-(3.3) の空間微分演算子に対して, 立方球面座標系における微分演算子 (付録 A を参照) を適用すれば, (ρ, ρu^{α}, ρu^{β}, ρu^{r}, ρe_{t}) に対する方程式系は, アインシュタインの縮約規則を用いて次のようになる.

\frac{∂ρ′}{∂t} + \frac{1}{J} \frac{∂}{∂x^{k}} (Jρu^{k}) = 0,  (3.12)
\frac{∂(ρu^{α})}{∂t} + \frac{1}{J} \frac{∂}{∂x^{k}} (J(ρu^{α}u^{k} + G^{αk}p′ )) = ψ_{H}^{α} + ψ_{M}^{α} + ψ_{C}^{α},  (3.13)
\frac{∂(ρu^{β})}{∂t} + \frac{1}{J} \frac{∂}{∂x^{k}} ((J(ρu^{β}u^{k} + G^{βk}p′ )) = ψ_{H}^{β} + ψ_{M}^{β} + ψ_{C}^{β},  (3.14)
\frac{∂(ρu^{r})}{∂t} + \frac{1}{J} \frac{∂}{∂x^{k}} ((J(ρu^{r}u^{k} + G^{βk}p′ )) = ψ_{M}^{r} + ψ_{C}^{r} + ψ_{G}^{r},  (3.15)
\frac{∂(ρe_{t})′}{∂t} + \frac{1}{J} \frac{∂}{∂x^{k}} (J(ρe_{t} + p)u^{k}) = 0.  (3.16)

ここで, 添字 k は {α, β, r} で展開され, (x^{α}, x^{β}, x^{r}) = (α, β, r) となる. また, G^{ij} は (A.3) で定義される立方球面格子における共変計量テンソル, J は (A.5) で定義される立方球面座標系の変換ヤコビアンである. ψ_{H} = (ψ_{H}^{α}, ψ_{M}^{β}) は, 静水圧平衡にある圧力の水平変化によるソース項であり,

ψ_{H} = −\frac{δ^{2}}{(1 + X^{2})(1 + Y^{2})} (\begin{matrix}(1 + Y^{2})\frac{∂\overline{p}^{h}}{∂α} + XY \frac{∂\overline{p}^{h}}{∂β} \\ XY \frac{∂\overline{p}^{h}}{∂α} + (1 + X^{2})\frac{∂\overline{p}^{h}}{∂β} \end{matrix} )  (3.17)

で与えられる. また, ψ_{M} = (ψ_{M}^{α}, ψ_{M}^{β}, ψ_{M}^{r}) は, 非直交座標系である球面立方座標系を導入したことによって生じるソース項で,

ψ_{M}^{i} = Γ_{jk}^{i}(ρu^{j}u^{k}),  (3.18)

によって与えられる. ただし, Γ_{jk}^{i} は第二種クリストッフェル記号であり, 等角立方球面座標系において次のように得られる.

Γ^{α} = (\begin{matrix}\frac{2XY^{2}}{δ^{2}} & \frac{−Y (1+Y^{2})}{δ^{2}} & \frac{1}{r} \\ \frac{−Y (1+Y^{2})}{δ^{2}} & 0 & 0 \\ \frac{1}{r} & 0 & 0 \end{matrix}) ,  (3.19)

Γ^{β} = (\begin{matrix}0 & \frac{−X(1+X^{2})}{δ^{2}} & 0 \\ \frac{−X(1+X^{2})}{δ^{2}} & \frac{2X^{2}Y}{δ^{2}} & \frac{1}{r} \\ 0 & \frac{1}{r} & 0 \end{matrix}) ,  (3.20)

Γ^{r} = \frac{r(1 + X^{2})(1 + Y^{2})}{δ^{4}} (\begin{matrix}−(1 + X^{2}) & XY & 0 \\ XY & −(1 + Y^{2}) & 0 \\ 0 & 0 & 0 \end{matrix}) .  (3.21)

一方, ψ_{C} = (ψ_{C}^{α}, ψ_{C}^{β}, ψ_{C}^{r}) は, コリオリ効果による強制であり,

ψ_{C} = −2Ω × (ρu)  (3.22)

あるいは,

(ψ_{C})^{i} = −2JG^{ij}ϵ_{jkl}Ω^{k}(ρu^{l})  (3.23)

である. ここで, ϵ_{ijk} はレヴィ=チヴィタ記号, Ω は惑星の自転角速度ベクトルであり, 地理座標系の基底 {e_{λ}, e_{ϕ}, e_{r}} を用いて,

Ω = ω[cos ϕe_{ϕ} + sin ϕe_{r}]  (3.24)

と書かれる. ただし, ω = |Ω| である. また, 等角立方球面座標系において Ω は, 次のように書かれる.

Ω = \begin{cases}ω [\frac{δ}{r(1+Y^{2})}e_{β} + \frac{Y}{δ} e_{r}] & for equatorial patch(n_{p} ≤ 4) \\ sω [−\frac{Xδ}{r(1+X^{2})}e_{α} − \frac{δ}{r(1+Y^{2})}e_{β} + \frac{1}{δ}e_{r}] & for polar patch(n_{p} = 5, 6) \end{cases} (3.25)

によって与えられる. ここで, s はパッチの指示子であり, 北極を含むパッチ (n_{p} = 5) に対して s = 1, 南極を含むパッチ (n_{p} = 6) に対して s = −1 である. 最後に, 放射方向の運動量方程式中の重力項 ψ_{G}^{r} は,

ψ_{G}^{r} = −ρ′g (\frac{a}{r})^{2}  (3.26)

で表される.

今, 保存変数を並べた解ベクトル

q = (ρ, ρu^{α}, ρu^{β}, ρu^{r}, ρe_{t})^{T}  (3.27)

を導入する. このとき, 非静力学方程式系 (3.12)-(3.16) は, この解ベクトル q を用いて簡潔に書ける. すなわち,

\frac{∂}{∂t} q + \frac{1}{J} \frac{∂}{∂x^{k}} F_{k} = Ψ_{H} + Ψ_{M} + Ψ_{C} + Ψ_{G}.  (3.28)

ここで, F_{k} は非粘性流束の各方向成分であり,

F_{k} = (\begin{matrix}ρu^{k} \\ ρu^{α}u^{k} + G^{αk}p′ \\ ρu^{β}u^{k} + G^{βk}p′ \\ ρu^{r}u^{k} + G^{rk}p′ \\ (ρe_{t} + p)u^{k} \end{matrix})  (3.29)

と定義し,

F = F_{α} e_{α} + F_{β} e_{β} + F_{r} e_{r}  (3.30)

と書くことにする. また, ψ_{H}, ψ_{M}, ψ_{C}, ψ_{G} をそれぞれ, Ψ_{H} = (0, ψ_{H}^{α}, ψ_{H}^{β}, 0, 0)^{T} , Ψ_{M} = (0, ψ_{M}^{α}, ψ_{M}^{β}, ψ_{M}^{r}, 0)^{T} , Ψ_{C} = (0, ψ_{C}^{α}, ψ_{C}^{β}, ψ_{C}^{r}, 0)^{T} , Ψ_{G} = (0, 0, 0, ψ_{G}, 0) と再定義した.

3.3 浅い大気近似を施した場合の定式化

地球大気のように惑星半径に対し大気層が十分薄い場合には, 浅い大気近似 (shallow atmosphere approxiamation) がよく行われる. DECAMO では, オプションとして浅い大気近似を施した非静力学方程式を支配方程式系として用いることもできる. この場合に, 計量テンソル G^{ij}, 変換ヤコビアン J に現れる r は惑星半径 a に置き換えられる. また, 浅い大気近似のもとでエネルギー保存則, 角運動量保存則が失われないようにメトリック項とコリオリ項の一部を無視する, 伝統的近似 (traditional approximation) を同時に行う. これらの近似が施されるとき, メトリック項における Γ_{jk}^{i} は,

Γ_{αr}^{α} = Γ_{rα}^{α} = 0, Γ_{βr}^{β} = Γ_{rβ}^{β} = 0, Γ^{r} = 0.  (3.31)

となる. また, 自転角速度ベクトル (3.25) の水平方向の寄与は無視される.

その結果, ソース項は次のように簡単化される. 初めに, メトリック項 (3.18) は,

(ψ_{M})_{shallow} = \frac{2}{δ^{2}} (\begin{matrix}−XY^{2}ρu^{α}u^{α} + Y (1 + Y^{2})ρu^{α}u^{β} \\ X(1 + X^{2})ρu^{α}u^{β} − X^{2}Y ρu^{β}u^{β} \end{matrix})  (3.32)

となる. また, コリオリ力の強制 (3.23) は, 赤道周りのパッチ (n_{p} = 1, 2, 3, 4) に対しては,

(ψ_{C})_{shallow} = \frac{2ω}{δ^{2}} (\begin{matrix}−XY^{2} & Y (1 + Y^{2}) \\ −Y (1 + X^{2}) & XY^{2} \end{matrix}) (\begin{matrix}ρu^{α} \\ ρu^{β} \end{matrix}) ,  (3.33)

極のパッチ (n_{p} = 5, 6) に対しては,

(ψ_{C})_{shallow} = \frac{2sω}{δ^{2}} (\begin{matrix}−XY & (1 + Y^{2}) \\ −(1 + X^{2}) & XY \end{matrix}) (\begin{matrix}ρu^{α} \\ ρu^{β} \end{matrix})  (3.34)

となる. 最後に, 浅い大気近似において重力ポテンシャルは,

(Φ)_{shallow} = gz,  (3.35)

したがって重力項 (3.26) は,

(ψ_{G}^{r})_{shallow} = −ρ′g  (3.36)

と簡単化される.

第4章数値解法

本章では, 非静力学方程式 (3.12)-(3.16) を数値的に解くために用いられる手法を示す. 詳細に触れる前に, 初めにその数値解法の概要を述べておく.

非静力学方程式 (3.12)-(3.16) は, 時間微分と空間微分を含む双曲型の偏微分方程式系 (partial differential equations, PDE) である. この双曲型偏微分方程式を数値的に解くために, 各格子セルに対してこれらの方程式に含まれる空間微分を離散化し, 偏微分方程式 (3.12)-(3.16) の初期値問題を, 各格子セル上で時間方向に対する常微分方程式 (ordinary differential equation: ODE) を解く問題に帰着させる, 線の方法 (method of lines) を用いる.

本モデルの空間離散化は, 有限体積法によって行われる. 有限体積法では, まず保存形式で記述された支配方程式系を任意の検査セルに対し体積分し, フラックスの発散・収束項を検査セル界面におけるフラックスの面積分に変換することから始める. 空間高精度化は, 検査セル内の物理量の空間分布を高次補間する (再構築) ことによって達成される. この再構築の手法には, Central ENO 法が適用される. これにより原理的には任意の空間精度を実現できるが, 本研究では結果的に水平方向に空間 3 次精度, 鉛直方向に空間 2 次精度となる再構築を行った. その後, 検査セル界面におけるフラックスを一意に決定するための数値流束関数を導入する. 本研究ではこの数値流束関数として, AUSM^{+}-up スキームが適用される.

次に, 時間積分の手法について述べられる. 大気大循環モデルでは通常, 水平方向に比べて鉛直方向の格子幅がずっと小さい. そのため, 陽解法の時間積分スキームにおいて, CFL 条件を満たす時間ステップ幅は厳しく制約される. この制約を緩和して効率的に時間積分を行うために, 水平方向の流体運動を陽的に, 鉛直方向の流体運動を陰的にそれぞれ分別的に取り扱う. 分別的に取り扱われた水平方向と鉛直方向の運動は, 時間に関して 2 次精度を保証する Strang-carryover 法によって組み合わされる.

4.1 有限体積法による空間離散化

DECAMO では, 非静力学方程式系 (3.28) の空間離散化に, 高精度有限体積法を用いる. 初めに, 立方球面格子の各パッチに分布する水平格子セルの中心および鉛直格子セルの中心の座標を, 整数の添字を使って

α_{i} = (i − \frac{1}{2})Δ − \frac{π}{4} , β_{j} = (j − \frac{1}{2})Δ − \frac{π}{4} , z_{k} = (k − \frac{1}{2})Δz  (4.1)

と記述することにする. ここで, 水平方向のインデックス (i, j) は [1, N_{c}]^{2} の範囲から, 鉛直方向のインデックス k は [1, N_{r}] の範囲から選ばれる. N_{r} は放射方向の鉛直レベル数である. 今後, 各パッチ内において, (α_{i}, β_{j}, z_{k}) を中心とする検査体積を, インデックス (i, j, k) を添える (例えば, Z_{ijk}) ことによって指定することにする.

次に, 支配方程式系 (3.28) を検査セル Z_{i,j,k}(その体積を |Z| と書く) に渡って体積積分することにより, 検査セル Z_{i,j,k} に対して体積平均された (3.28) の半離散化式を導く. フラックス項の体積積分は, ガウスの発散定理により, 検査体積の境界 ∂Z_{i,j,k} に対する面積分に書き換える. これにより,

\frac{d}{dt} \overline{\overline{q}}_{i,j,k} + \frac{1}{|Z|_{i,j,k}} ∫∫_{∂Z_{i,j,k}} F · n dS = \overline{\overline{Ψ_{H}}}_{i,j,k} + \overline{\overline{Ψ_{M}}}_{i,j,k} + \overline{\overline{Ψ_{C}}}_{i,j,k} + \overline{\overline{Ψ_{G}}}_{i,j,k} (4.2)

を得る. ここで, 非粘性流束の積分は, 法線ベクトル n と微小面素 dS をもつ面に渡って積分される. なお, F · n はベクトル量であり, 検査体積の境界面を垂直に横切るそれぞれの状態変数の外向きフラックスを表している. 二重上線は, 次の形式の三次元体積平均量を表してる.

\overline{\overline{ϕ}}_{i,j,k} = \frac{1}{|Z|_{i,j,k}} ∫_{Z_{i,j,k}} ϕ dV  (4.3)

ここで, dV = Jdαdβdr は検査体積の微小体積素である. したがって, \overline{\overline{q}}_{i,j,k} は検査セル Z_{i,j,k} における状態ベクトルの体積積分値を意味する. 最後に, \overline{\overline{Ψ_{H}}}, \overline{\overline{Ψ_{M}}}, \overline{\overline{Ψ_{C}}}, \overline{\overline{Ψ_{G}}} は, それぞれ, 圧力の静水圧部分の水平変化による寄与, 球面立方格子の幾何的な寄与, コリオリ効果の寄与, 重力の寄与によるソース項の体積積分値を表している.

今, 検査セルの界面に渡る面積分による平均量を一本上線で表すことにする. 例えば, 検査セル Z_{ijk} の界面 (i + \frac{1}{2}, j, k) を横切るフラックスの平均値を, \overline{F}_{i+ \frac{1}{2} ,j,k} と書き,

\overline{F}_{i+ \frac{1}{2} ,j,k} = \frac{1}{|∂Z|_{i+ \frac{1}{2} ,j,k}} ∫∫_{∂Z_{i+ \frac{1}{2} ,j,k}} F · n dS  (4.4)

のように定義する. したがって, 界面 (i + \frac{1}{2}, j, k) を横切る全フラックス (F_{i+ \frac{1}{2} ,j,k} と書く) は, フラックスの面平均値と界面の面積を使って単に,

F_{i+ \frac{1}{2} ,j,k} = \overline{F}_{i+ \frac{1}{2} ,j,k} |∂Z|_{i+ \frac{1}{2} ,j,k}  (4.5)

となる.

(4.5) のように各方向の検査セル界面を横切る全フラックスを書くことにすれば, 各検査セルに渡って積分された支配方程式系 (4.2) は次のように書き直すことができる.

\frac{d}{dt} \overline{\overline{q}}_{i,j,k} = H_{i,j,k}(q) + V_{i,j,k}(q).  (4.6)

ここで, 水平方向の運動と鉛直方向の運動を分離し, それらを表す項を以下のように, H, V とおいた.

H_{i,j,k}(q) = −\frac{1}{|Z|_{i,j,k}} [F_{i+ \frac{1}{2} ,j,k} − F_{i− \frac{1}{2} ,j,k} + F_{i,j+ \frac{1}{2} ,k} − F_{i,j− \frac{1}{2} ,k}]
+\overline{\overline{Ψ_{H}}} + \overline{\overline{Ψ_{M}}} + \overline{\overline{Ψ_{C}}},  (4.7)
V_{i,j,k}(q) = −\frac{1}{|Z|_{i,j,k}} [F_{i,j,k+ \frac{1}{2}} − F_{i,j,k− \frac{1}{2}} ] + \overline{\overline{Ψ_{G}}}.  (4.8)

大気循環シミュレーションにおいて, 多くの場合格子セルの鉛直方向の幅は水平方向に比べ相対的に小さい (10^{−2}〜10^{−3} 倍) ため, CFL 条件を介して時間ステップ幅は厳しく制限される. 故に, 数値的に安定な時間ステップ幅は鉛直格子セルの幅に非依存にするために, 時間積分に関して水平方向の運動による寄与 H_{i,j,k} を陽的, 鉛直方向の運動による寄与 V_{i,j,k} を陰的に取り扱う.

(4.6) に現れるフラックス項の面積分やソース項の体積積分は, ガウス求積法 (Gaussian quadrature method) を用いて, 数値的に求められる. ガウス求積法は, 被積分関数の評価点 x_{i}(i = 1, ⋯, n) とそれに対応する重み w_{i} を適切に選ぶことで, 定積分を 2n − 1 次 (あるいはそれ以下) の多項式の積分結果と一致させる. ガウス求積法の公式は, 習慣的に [−1, 1] の積分区間に対し, 次のように書かれる.

 ∫_{−1}^{1} f(x)dx ≈ ∑_{i=1}^{n} w_{i}f(x_{i})  (4.9)

空間精度に関して K 次精度を保証するためには, K − 1 次以上の多項式によるガウス求積法を用いる必要がある. 本研究では鉛直方向の運動に関しては, 空間精度は 2 次精度までしか求めないので, 鉛直方向のセル界面における面積分の求積点は一点でよい. 一方, 水平方向の運動に関しては, 空間精度 3 次精度を要求するので, 水平方向のセル界面の面積分には複数の求積点を必要とする. 結果的に, 検査セルの水平方向の界面 Z_{i± \frac{1}{2} ,j± \frac{1}{2} ,k} におけるフラックス項の面積分は, 各々の面で 2 個の求積点を用いる. 例えば, 界面 Z_{i+ \frac{1}{2} ,j,k} の面積分に対するこの 2 個の求積点の位置は, \{(α_{i+ \frac{1}{2}} , β_{j} − Δ\sqrt{3}/6, z_{k}), (α_{i+ \frac{1}{2}} , β_{j} + Δ\sqrt{3}/6, z_{k})\} である. 一方, 鉛直方向の界面 Z_{i,j,k± \frac{1}{2}} におけるフラックス項の面積分の求積点は, 点 (α_{i}, β_{j}, z_{k± \frac{1}{2}} ) の 1 点で良い. また, ソース項の体積積分では, 計 4 点 (α_{i} ± Δ\sqrt{3}/6, β_{j} ± Δ\sqrt{3}/6, z_{k}) の求積点を用いる. なお, 各評価点に対する重みは同じで良い. 図 (4.1) は, 検査セル Z_{ijk} に対する面積分や体積分を上のガウス求積法で求めるとき, 使われる求積点の位置を示している.

図 4.1: 検査セル Z_{i,j,k} に対するフラックス項の面積分およびソース項の体積分をガウス求積法で求めるときに用いる求積点の位置. 面積分の各求積点は記号 S, 体積分の各求積点は記号 V を用いて表している. 曲線 l_{α±} は, それぞれ, z = z_{k} の球面と α = α_{i} ± Δ/(2\sqrt{3}) の曲面の交点である. 同様に, 曲線 l_{β±} は, それぞれ, z = z_{k} の球面と β = β_{i} ± Δ/(2\sqrt{3}) の曲面の交点である. z = z_{k} の球面上に存在する各求積点 S は曲線 l_{α±} または曲線 l_{β±} 上に存在する. また, 各求積点 V は l_{α±} と曲線 l_{β±} の交点である.

4.2 再構築

検査セルに対するフラックス項の面積分やソース項の体積積分を高精度に行うためには, 検査セル内における状態変数の空間的な変化を考慮する必要がある. しかし, 有限体積法を用いる場合, 各検査セル内では体積平均量しか知り得ないので, 隣接する検査セルの情報を用いて多項式を生成し, 検査セル内の状態変数の空間的な変化を近似する. 一般に, K 次多項式による再構築は, なめらかな問題に対して K + 1 次精度の空間離散化を与える.

本研究では, この再構築の手法として, 水平方向に対しては Ivan (2010); Ivan et al. (2012) が提案した CENO(Central ENO) 法を用いる.

4.2.1 Central ENO 法による水平方向の再構築

CENO 法は, なめらかな領域では高精度を得るために, k-exact 最小二乗法に基いて各検査セルに対し区分的な再構築を行う. また, なめらかでない解が予想される領域において解の単調性を制御するために, 制限的に線形再構築を併用する.

検査セル I(i, j, k)(= Z_{i,j,k}) 内の任意の位置における解変数 q_{I}(α, β) の変化は, 次の多項式で近似される.

q_{I}^{K}(α, β) = ∑_{\begin{matrix}p_{1},p_{2}=0 \\ (p_{1}+p_{2}≤K) \end{matrix}}^{K}(α − α_{I})^{p_{1}}(β − β_{I})^{p_{2}}D_{p_{1}p_{2}}^{K}.  (4.10)

ここで, K は多項式の次数, (α_{I}, β_{I}) は (等角座標空間における) 検査セル Z_{i,j} の中心位置で α_{I} = α_{i}, β_{I} = β_{j} である. また, D_{p_{1}p_{2}}^{K} は高次多項式の係数であり, Z_{i,j,k} とその近傍の検査セルにおける体積平均値を用いて決定される. 本研究では空間 3 次精度スキームを目的とするため K = 2 の場合のみ扱うが, (4.10) は任意精度の再構築を可能にする. なお, 単調性の維持のために, 解像限界または不連続が予想される流れの領域において, 多項式の次数 K は 1 まで落とされ, リミッタが適用される. その基準については, 後ほど述べる.

決定されるべき係数 D_{p_{1}p_{2}}^{K} の数は, N_{D} = (K + 1)(K + 2)/2 によって与えられる. したがって, K = 2 の場合には係数の数は 6 個である. 係数 D_{p_{1}p_{2}}^{K} の決定は, (4.10) から得られる値がセル I(i, j, k) のステンシルに含まれる各検査セル内の平均値と一致することを課すことで導かれる, D_{p_{1}p_{2}}^{K} に対する線形過剰決定系を解くことによってなされる. この D_{p_{1}p_{2}}^{K} に対する線形過剰決定系を導くにあたって, 再構築の手法に課す条件は, 具体的には次のようになる.

1. 解の再構築は, 正確に N ≤ K 次の多項式を再現しなければならない.

2. 解の再構築は, 検査セル内の平均値を保存しなければならない.

3. 解の再構築に用いるステンシルは, 検査セル I(i, j, k) にトポロジー的に近傍でかつその中心に物理的に近い検査セルによって構成されなければならない (コンパクトサポート).

検査セル I(i, j, k) に対して, 一番目の条件は「K-完全性 (K-exactness)」の性質を表している. すなわち, 解が次数 K あるいはそれ以下の多項式である場合には, 解が正確に再構築されることを意味する. 2 番目の条件は, 再構築された多項式関数が検査セル I(i, j, k) の体積平均値を正確に与えることを課す. すなわち,

\overline{\overline{q}}_{i,j,k} = \frac{1}{|Z_{i,j,k}|} ∫_{Z_{i,j,k}} q_{I}^{K}(α, β) dV.  (4.11)

ただし, dV = Jdαdβ である.

コンパクト・ステンシルの最小サイズは, 必要とされる係数の数 N_{D} によって決まる. しかし実際には, メッシュの不規則や格子と相対的な方向による解の補間精度の不均一に対応するために, さらに近傍のセルがステンシルに付け加えられる. 本研究では, K = 2 の k-exact 法に用いるステンシルを, 検査セル I(i, j, k) およびその近接の検査セルの計 9 個で構成する. ただし, パッチの角にある検査セルのステンシルは, 計 8 個で構成する.

検査セル I(i, j, k) に対する再構築において, 厳格に満たすべき (4.11) の制約と伴に, 線形過剰決定系 𝔸D − B = 0 が最小二乗問題的に解かれる. ここで, D は多項式の係数 D_{p_{1}p_{2}}^{K} を成分とするベクトルである. 方程式 𝔸D − B = 0 は, 検査セル I(i, j, k) のステンシルに含まれる検査セル J(γ, δ, k) に対して,

(𝔸D − B)_{γ,δ,k} = (\frac{1}{|Z_{γ,δ,k}|} ∫_{Z_{γ,δ,k}} q_{I}^{K}(α, β) dV − \overline{\overline{q}}_{γ,δ,k}) = 0  (4.12)

によって与えられる. これらの方程式は, 実際の検査セル J(γ, δ, k) の体積平均値 \overline{\overline{q}}_{J} と, 再構築される多項式 q_{I}^{K}(α, β) の検査セル J(γ, δ, k) に渡る体積平均値を一致させる. (4.11) を解析的に満たさせるために, (4.11) の q_{I}^{K} を (4.10) で置き換え, 最初の係数 D_{00}^{K} を他の M = N_{D} − 1 個の未知の項で表せば,

D_{00}^{K} = \overline{\overline{q}}_{I} − ∑_{\begin{matrix}p_{1},p_{2}=0 \\ (p_{1}+p_{2}≤K) \\ (p_{1}+p_{2}≠0) \end{matrix}}^{K} D_{p_{1}p_{2}}^{K}\overline{(α^{p_{1}}β^{p_{2}})}_{I}  (4.13)

を得る. ここで, 冪乗 (p_{1}, p_{2}) の幾何モーメント\overline{(α^{p_{1}}β^{p_{2}})}_{I} を,

\overline{(α^{p_{1}}β^{p_{2}})}_{I} = \frac{1}{|Z_{I}|} ∫_{Z_{I}}(α − α_{I})^{p_{1}}(β − β_{I})^{p_{2}} dV.  (4.14)

と定義する.

(4.12) の q_{I}^{K} に (4.10) を代入し, D_{00}^{K} に対し (4.13) を用いれば, 次の M 個の未知数に対する線形過剰決定系を得る.

[ L_{1} ; L_{2} ; ⋮ ; L_{J} ; ⋮ ; L_{N_{n}} ]_{N_{n}×M} ( D_{01}^{K} ; D_{02}^{K} ; ⋮ ; D_{p_{1}p_{2}}^{K} ; ⋮ ; D_{K0}^{K})_{M×1} − ( w_{1}(\overline{\overline{q}}_{1} − \overline{\overline{q}}_{I}) ; w_{2}(\overline{\overline{q}}_{2} − \overline{\overline{q}}_{I}) ; ⋮ ; w_{J}(\overline{\overline{q}}_{J} − \overline{\overline{q}}_{I}) ; ⋮ ; w_{N_{n}}(\overline{\overline{q}}_{N_{n}} − \overline{\overline{q}}_{I}))_{N_{n}×1} = (0 0 ⋮ 0 ⋮ 0)_{N_{n}×1}  (4.15)

ここで, N_{n} はステンシル内の近傍検査セルの数である. また, 近傍検査セルJ(γ, δ, k) の方程式と対応する係数行列 𝕃 の第 J 行目の列 L_{J} は,

L_{J} = (w_{J}({ \atop α^{0}β^{1}})_{IJ} w_{J}({ \atop α^{0}β^{2}})_{IJ} ⋯ w_{J}({ \atop α^{p_{1}}β^{p_{2}}})_{IJ} ⋯ w_{J}({ \atop α^{K}β^{0}})_{IJ}) (4.16)

によって与えられる. w_{J} はそれぞれの検査セル J に特有の幾何的な重みであり,

w_{J} = \frac{1}{|x_{J} − x_{I}|^{θ}}  (4.17)

と定義される. ここで, x_{J}, x_{I} はそれぞれ近傍検査セル J と検査セル I の位置ベクトルを表す. また, θ = 0, 1, 2 である. この幾何的な重みは, 再構築の局所性を改善するために導入される. 検査セル I と検査セル J の組に対する行列要素の係数 ({ \atop α^{p_{1}}β^{p_{2}}})_{IJ} は, 次の形式で表される.

({ \atop α^{p_{1}}β^{p_{2}}})_{IJ} = (\frac{1}{|Z_{J}|} ∫_{Z_{J}} (α − α_{I})^{p_{1}}(β − β_{I})^{p_{2}} dV ) − \overline{(α^{p_{1}}β^{p_{2}})}_{I}  (4.18)

全ての ({ \atop α^{p_{1}}β^{p_{2}}})_{IJ} は, 付録 B に示すように, \overline{(α^{p_{1}}β^{p_{2}})} のみを使って効率的に求まる.

線形過剰決定系 (4.15) の解は, QR 分解または特異値分解によって求めた 𝕃 の擬逆行列の乗算によって得られる. 各検査セル, 各状態変数 (あるいはプリミティブ変数), この重み付き最小二乗問題が解かれる. 行列 𝕃 は完全に幾何的な情報だけに依存するので, ある検査セル (i, j, k) における全ての最小二乗問題に対して, いつでも同じである. したがって, 行列 𝕃 の逆行列は前もって計算可能であり, 再利用できる. もし擬逆行列を保持するためのメインメモリ容量に問題がないならば, 擬逆行列を再利用することで大幅に計算速度を向上させることができる. 擬逆行列のサイズは, 二次多項式補間 (k = 2) で N_{n} = 8 の場合に 5 × 8, 三次多項式補間 (k = 3) で N_{n} = 12 の場合に 9 × 12 である. したがって, 行列要素を倍精度で保持する場合に必要なメモリ量は各検査セルに対し, それぞれ 320 Byte, 864 Byte 必要である. 地形を含まない場合, 立方球面格子の各パッチ, 各層の同じ水平位置においてこの擬似逆行列は共通なので, 例えば, c90 の水平解像度 (水平格子間隔約 1 度に対応) に対して擬似逆行列を保持するために必要な全メモリ量は, k = 2 で約 3 MByte, k = 3 で約 7 MByte であり, 現在の計算機のメインメモリ容量では問題とならない. したがって, 本研究ではあらかじめ擬逆行列を計算しておき, 再利用する.

CENO 法における単調性の保証は, 解像限界または不連続が予想される流れの領域において, 非制限の k-exact 法による再構築を k = 1 の制限付き線形再構築まで落とすことによってなされる. 再構築の次数を落とすことよる解の単調性の維持は, 解像領域, 非解像領域ともに非一様な精度の再構築をもたらす. しかし, 解がどこでも滑らかである限り, 滑らかな極値が存在する場合でも一様な精度を維持できることが知られている. 解が計算セルに解像されているかを判別するためには, 次に定義される smoothness indicator 𝕊 を用いる.

𝕊 = \frac{α}{max(1 − α, ε)} \frac{N_{n} − N_{D}}{N_{D} − 1} .  (4.19)

ここで, a は,

α = 1 − \begin{matrix}∑ & J & (q_{J}^{k}(X_{J}) − q_{I}^{k}(X_{J})) & 2 \\ ∑ & J & (q_{J}^{k}(X_{J}) − \overline{q}_{I}) & 2 \end{matrix}  (4.20)

のように決定される. 閾値 ε はゼロによる除算を回避するために導入される. ε の適切な値は, 10^{−8} であることが知られている. この smoothness indicator を評価した後, その値はカットオフの基準 𝕊_{c} と比較される. 𝕊 < 𝕊_{c} の場合, 解は解像限界あるいは不連続であると判断され, 高次の k-exact 再構築は制限付きの線形再構築に切り換えられる. 一方, 𝕊 > 𝕊_{c} の場合には, 非制限の高次再構築が適用される. カットオフ 𝕊_{c} の値は, 1000 から 5000 の範囲で選ばれる.

通常, 上述した単調性を保証する手法は, 密度変数や圧力変数に対し適用される. しかし, 物質 (水蒸気, 化学物質など) の輸送などを除き, 大気の大規模な流れに対し単調性の保証は必要でない場合が多い. 特に, 今回の非静力学コアの定式化では, 密度や圧力に対する再構築は静水圧成分からのずれに対し行うため, 密度や圧力が高次再構築により負となる問題は生じないと考えた. また計算効率の観点からも, 最終的に本研究では単調性の保証は考慮しないことにした.

4.2.2 鉛直方向の再構築

鉛直方向は, 2 次多項式による再構築を行う. すなわち, 検査セル I = (Z_{i,j,k}) 内の鉛直方向の変化は,

q_{I}(z) = ∑_{p=0}^{2} D_{p}(z − z_{I})^{p}  (4.21)

で近似される. ここで, z_{I} = z_{k} である. 係数 D_{p} は, q_{I}(z) を検査セル I および近接する検査セル Z_{i,j,k± \frac{1}{2}} における体積平均値と単に一致させることにより決定する. その結果, 係数 D_{p} は次のように決定される.

\begin{aligned}D_{0} & = q_{k}, \\ D_{1} & = −\frac{δ_{+}}{δ_{−}(δ_{−} + δ_{+})} q_{k−1} + \frac{δ_{+} − δ_{−}}{δ_{−}δ_{+}} q_{k} + \frac{δ_{−}}{δ_{+}(δ_{−} + δ_{+})} q_{k+1}, \\ D_{2} & = \frac{1}{δ_{−}(δ_{−} + δ_{+})} q_{k−1} − \frac{1}{δ_{−}δ_{+}} q_{k} + \frac{1}{δ_{+}(δ_{−} + δ_{+})} q_{k+1}. \end{aligned}  (4.22)

ここで, δ_{+} = z_{k+1} − z_{k}, δ_{−} = z_{k} − z_{k−1} である.

\overline{main.tex 2013/04/02(河合佑太)}

ただし, 上の方法において, 境界下端および上端に接する検査セルに対してはステンシルの一部が領域外にはみ出してしまう. 本研究では, これらの検査セルに対しても同様の再構築を適用できるように, ゴーストセルを導入する. ゴーストセル Z_{i,j,0}, Z_{i,j,N_{r}+1} の状態変数の値は, 境界端において境界条件が満たされるように選ばれる. ゴーストセルを導入することで, 境界近傍で再構築が線形になることを回避でき, また境界条件を精度良く取り扱える.

4.3 検査セル境界における数値流束の評価

風上有限体積法において, 検査セル界面におけるフラックス項は, その界面に接する左右のそれぞれのセルの再構築から得られた値を用いて評価される. しかし, 再構築はそれぞれの検査セルで区分的に行われるため, 界面の左側と右側で再構築される値は一般的には一致しない. この検査セル界面のフラックスを一意に決定する関数は数値流束関数と呼ばれ, それによって求められたフラックスは数値流束 (numerical flux) と言われる. この数値流束関数を与え, 同時にフラックスの風上化を行う手法として, Roe の流束差分離 (Flux difference splitting: FDS) 法 (Roe, 1981), Van Leer の流束ベクトル分離 (Flux vector splitting: FVS) 法が有名である. また近年では, 風上分離法 (Advective Upstream Splitting Method: AUSM)(Liou, 1996, 2006) が考案されている. AUSM 族のスキームは, 流束ヤコビアンの固有値の構造を陽には必要としないため, FDS 法や FVS 法に比べ, 計算効率が高く, 実装が容易である.

大気大循環のような低マッハ数流れに対する各数値流束関数の適用性は, Ullrich et al. (2010); Ullrich and Jablonowski (2012a,b) によって調べられている. Ullrich et al. (2010) では, 全球浅水モデルの標準テストケース (Williamson et al., 1992) を用いて, Rusanov スキーム, Roe スキーム, AUSM^{+}-up スキームの性能を比較し, AUSM^{+}-up スキームを用いた場合に全体的に最も精度の良い結果が得られることを確かめた. また, 彼らが行った浅水系のテストケースでは, Roe スキームと AUSM^{+}-up スキームの精度の差はあまり大きくないという結果が出た. しかし, Roe スキームなどの近似リーマン解法を低マッハ数流れに対してそのまま用いると, 解の精度や反復計算における収束性が悪化することが知られている. これを回避するためには, Weiss and Smith. (1995) などの前処理法を近似リーマン解法に適用する必要がある. 一方, Liou (2006) が開発した AUSM^{+}-up スキームは, 低マッハ数流れに対しても適用できるように拡張された AUSM 族スキームである. したがって, 本研究では, 計算効率と低マッハ数流れへの適用性の観点から, 数値流束の評価方法として AUSM^{+}-up スキームを用いることにする.

4.3.1 直交化

一般に数値流束の評価において, 数値流束関数に渡される速度ベクトルの成分は直交座標系でのベクトル成分である. したがって, 立方球面座標系における速度ベクトル成分 (u^{α}, u^{β}) は, 各検査セル界面において直交座標系での速度成分に変換されなければならない. 各検査セル界面における 3 次元直交座標系は, 対象とする界面に垂直な基底 e_{⊥} と, 界面に平行な 2 つの基底 e_{1}, e_{2} によって構成する. 直交基底が定義された点において, 任意のベクトル v は, 立方球面座標系における自然基底 {g_{α}, g_{β}, g_{r}} を使って,

v = v^{α}g_{α} + v^{β}g_{β} + v^{r}g_{r},  (4.23)

あるいは直交基底 {e_{⊥}, e_{1}, e_{2}} を使って,

v = v^{⊥}e_{⊥} + v^{1}e_{1} + v^{2}e_{2}  (4.24)

と書ける. これらの座標系間の変換は, 基底直交化行列 O_{d}(α, β, r) を介してなされる. すなわち,

(\begin{matrix}v^{⊥} \\ v^{1} \\ v^{2} \end{matrix}) = O_{d}(α, β, r) (\begin{matrix}v^{α} \\ v^{β} \\ v^{r} \end{matrix}) .  (4.25)

ここで, d ∈ α, β, r であり, 定数に固定される座標を表す. 同様に, 直交基底から自然基底への変換は, 逆行列 O_{d}^{−1}(α, β, r) を用いて,

(\begin{matrix}v^{α} \\ v^{β} \\ v^{r} \end{matrix}) = O_{d}^{−1}(α, β, r) (\begin{matrix}v^{⊥} \\ v^{1} \\ v^{2} \end{matrix}) .  (4.26)

球面立方座標系の自然基底 {g_{α}, g_{β}} は, これらの基底に対して垂直でかつ正規化された放射方向の基底を構成するので, 単に e_{2} = g_{r} とすればよい. したがって, 基底変換行列は 2 × 2 行列となる. 結果的に, 一定の α または β の検査セル界面における変換行列はそれぞれ,

O_{α} = (\begin{matrix}\frac{\sqrt{1 + X^{2}}}{δ} & 0 \\ \frac{−XY^{2}\sqrt{1 + X^{2}}}{δ^{2}} & \frac{(1 + Y^{2})\sqrt{1 + X^{2}}}{δ^{2}} \end{matrix} ) ,  (4.27)

O_{β} = (\begin{matrix}0 & \frac{\sqrt{1 + Y^{2}}}{δ} \\ \frac{(1 + X^{2})\sqrt{1 + X^{2}}}{δ^{2}} & \frac{−XY \sqrt{1 + Y^{2}}}{δ^{2}} \end{matrix} )  (4.28)

によって与えられる. 一方, 逆変換行列は,

\overline{main.tex 2013/04/02(河合佑太)}

O_{α}^{−1} = (\begin{matrix}\frac{δ}{\sqrt{1 + X^{2}}} & 0 \\ \frac{XY δ}{(1 + Y^{2})\sqrt{1 + X^{2}}} & \frac{δ^{2}}{(1 + Y^{2})\sqrt{1 + X^{2}}} \end{matrix} ) ,  (4.29)

O_{β}^{−1} = (\begin{matrix}\frac{XY δ}{(1 + X^{2})\sqrt{1 + Y^{2}}} & \frac{δ^{2}}{(1 + X^{2})\sqrt{1 + Y^{2}}} \\ \frac{δ}{\sqrt{1 + Y^{2}}} & 0 \end{matrix} )  (4.30)

によって与えられる. 上で定義された行列は, ベクトル v の水平成分に対してのみ作用する.

4.3.2 AUSM^{+}-up スキーム

最初に, 対象とする検査セル界面において, 界面に垂直な基底 e_{⊥} と, 界面に平行な 2 つの基底 e_{1}, e_{2} によって構成される, 3 次元直交座標系を導入する. この直交座標系において速度ベクトルを v = v^{⊥}e_{⊥} + v^{1}e_{1} + v^{2}e_{2} と書くことにする. この各方向の速度成分 (v^{⊥}, v^{1}, v^{2}) を使った状態ベクトルを,

Q = (ρ, ρv^{⊥}, v^{1}, v^{2}, ρe_{t})^{T}  (4.31)

また e_{⊥} 方向の非粘性流束を

F_{inv} = (ρv^{⊥}, ρv^{⊥}v^{⊥} + p, ρv^{1}v^{⊥}, ρv^{2}v^{⊥}, ρHv^{⊥})^{T}  (4.32)

と定義する. ここで, H = e_{t} + p/ρ は全エントロピーである.

AUSM 族のスキームは共通して, 第一段階として非粘性流束を対流フラックスと圧力フラックスに陽に分割する. すなわち,

F_{inv} = F^{(c)} + P = ṁψ⃗ + P .  (4.33)

ここで, ṁ = ρv^{⊥}, ψ⃗ = (1. v^{⊥}, v^{1}, v^{2}, H)^{T} である. 対流フラックスは保存則のシステム全体でスカラーの質量フラックス ṁ を共有し, ベクトル量 ψ⃗ は ṁ によって移流される点に注意が必要である. 圧力フラックスは圧力項のみを含む.

P = (0, p, 0, 0, 0)^{T} .  (4.34)

数値流束 f_{1/2} は, 離散化前の (4.33) の表式を保ちながら, 共通の質量フラックスを用いて,

f_{1/2} = ṁ_{1/2}ψ⃗_{L/R} + p_{1/2}  (4.35)

と書ける. なお, セル界面における量であることを添字 1/2 で表している. さらに, ψ⃗_{L/R} は簡単な風上化の原理に基づき,

ψ⃗_{L/R} = \begin{cases}ψ⃗_{L} & (ṁ_{1/2} > 0) \\ ψ⃗_{R} & otherwise \end{cases}  (4.36)

によって決定される. ψ⃗_{L}, ψ⃗_{R} は, 左右の検査セルにおける再構築によって得られるセル界面での値である.

セル界面の質量フラックス ṁ_{1/2} は,

ṁ_{1/2} = v_{1/2}^{⊥} ρ_{L/R} = a_{1/2}M_{1/2} \begin{cases}ρ⃗_{L} & (v_{1/2}^{⊥} > 0) \\ ρ⃗_{R} & otherwise. \end{cases}  (4.37)

の形式で書かれる. ここで, v_{1/2}^{⊥}, M_{1/2}, a_{1/2} はそれぞれ, 界面における対流速度, マッハ数, 音速である. ρ_{L/R} は, v_{1/2}^{⊥} によって移流される密度であり, 風上化の概念により v_{1/2}^{⊥} の方向によって規定される. 作業変数として導入されたマッハ数 M_{1/2} の表し方は, AUSM 族の各スキームにより異なる. その定式化にあたる理論的な背景や低マッハ数流れのための全速度拡張の詳細については, Liou (2006) を参照されたい. ここでは, AUSM^{+}-up スキームにおける M_{1/2} や p_{1/2} の表式を単に示すだけにする. 初めに, セル界面左右のマッハ数を,

M_{L/R} = \frac{v_{L/R}^{⊥}}{a_{1/2}} ,  (4.38)

局所的な平均マッハ数を

\overline{M}^{2} = \frac{M_{L}^{2} + M_{R}^{2}}{2}  (4.39)

と定義する. また, セル界面における音速 a_{1/2} は, ここでは簡単に左右の平均として a_{1/2} = (a_{L} + a_{R})/2 によって定義することにする. このとき全速度拡張された AUSM^{+}-up スキームにおけるセル界面のマッハ数は,

M_{1/2} = M_{(4)}^{+}(M_{L}) + M_{(4)}^{−}(M_{R}) − \frac{K_{p}}{f_{a}} max(1 − σ\overline{M}^{2}, 0)\frac{p_{R} − p_{L}}{ρ_{1/2}a_{1/2}^{2}}  (4.40)

と定義される. ただし,

M_{(2)}^{±}(M) = \frac{1}{4}(M ± 1)^{2}, M_{(2)}^{±}(M) = \begin{cases}\frac{1}{2}(M ± |M|) & (|M| ≥ 1) \\ M_{(2)}^{±}(M)(1 ∓ 16β̃M_{(2)}^{∓}(M)) & (otherwise) \end{cases} (4.41)

一方, 全速度拡張された AUSM^{+}-up スキームにおけるセル界面の圧力は,

p_{1/2} = P_{(5)}^{+}(M_{L})p_{L} + P_{(5)}^{−}(M_{R})p_{R} − K_{u}P_{(5)}^{+}P_{(5)}^{−}(ρ_{L} + ρ_{R})(f_{a}a_{1/2})(v_{R}^{⊥} − v_{L}^{⊥}) (4.42)

と定義される. ただし,

P_{(5)}^{±}(M) = \begin{cases}\frac{1}{2}(1 ± sign(M)) & (|M| ≥ 1) \\ M_{(2)}^{±}(M) [(±2 − M) ∓ 16αMM_{(2)}^{∓}(M)] & (otherwise). \end{cases}  (4.43)

上で与えられた M_{1/2}, p_{1/2} には, 低マッハ数流れに対して拡散項をスケーリングするための関数

f_{a}(M_{o}) = M_{o}(2 − M_{o}) ∈ [0, 1], M_{o}^{2} = min(1, max(\overline{M}^{2}, M_{∞}^{2})) ∈ [0, 1] (4.44)

が含まれる. このスケーリング関数は, 低マッハ数において圧力の小さな振動が最低次の速度場に悪影響を及ぼさないために導入されたここで, M_{∞} は無限大マッハ数である. また,

α = 3(−4 + 5f_{a}^{2})/16, β = 1/8.  (4.45)

最後に, M_{1/2}, p_{1/2} には, フリーパラメータ K_{p}, K_{u}, σ ∈ [0, 1] が存在するが, Liou (2006) では,

K_{p} = 0.25, K_{u} = 0.75, σ = 1.0

と選ばれている. 本研究では上の値を基準値とし, それらを (K_{p})_{ref}, (K_{u})_{ref} と表すことにする.

しかし, DECAMO の定式化において, 圧力勾配項は摂動場の圧力 p′(静水圧平衡からのずれ) を使って書かれるため, スケーリング関数 f_{a} は必要ない. f_{a} = 1 とすれば AUSM^{+}-up スキームは全速度拡張から基本スキームに戻るので, (4.40), (4.42) をそれぞれ,

M_{1/2} = M_{(4)}^{+}(M_{L}) + M_{(4)}^{−}(M_{R}) − K_{p} max(1 − σ\overline{M}^{2}, 0)\frac{p′_{R} − p′_{L}}{ρ_{1/2}a_{1/2}^{2}} ,  (4.46)
p′_{1/2} = P_{(5)}^{+}(M_{L})p′_{L}+P_{(5)}^{−}(M_{R})p′_{R}−K_{u}P_{(5)}^{+}P_{(5)}^{−}(ρ_{L}+ρ_{R})(a_{1/2})(v_{R}^{⊥}−v_{L}^{⊥})   (4.47)

と書き換えることにする.

4.3.3 鉛直方向の数値流束のための AUSM^{+}-up スキームの凖線形化

総観スケールの大気の運動において鉛直速度は非常に小さい (cm/s オーダ) ので, AUSM^{+}-up スキームを線形化することによって, 数値流束関数はより簡単化される. この手順は, Ullrich and Jablonowski (2012a) において説明されている. この線形化により鉛直方向のフラックス項の離散化形式がずっと簡単になるため, 陰的な鉛直方向の時間積分で用いられる非線形方程式の解析的なヤコビアンを, 比較的簡単な形式で得ることができる.

凖線形化された AUSM^{+}-up スキームでは, 鉛直方向の検査セル界面における線形化の基準音速を, 静水圧大気の音速\overline{a}_{1/2}^{h} = \sqrt{γ\overline{p}_{1/2}^{h}/\overline{ρ}_{1/2}^{h}} によって与える. この凖線形化の結果, 鉛直方向の検査セル界面における界面に垂直な速度 v_{1/2}^{⊥} は,

w_{1/2} = v_{1/2}^{⊥} = \frac{(ρv^{⊥})_{L} + (ρv^{⊥})_{R}}{ρ_{L} + ρ_{R}} − (K_{p})_{V} \frac{p′_{R} − p′_{L}}{\overline{ρ}_{1/2}^{h}\overline{a}_{1/2}^{h}} ,  (4.48)

検査セル界面における圧力は,

p′_{1/2} = \frac{p′_{L} + p′_{R}}{2} − (K_{u})_{V} \frac{\overline{ρ}_{1/2}^{h}\overline{a}_{1/2}^{h}}{2} ((v^{⊥})_{R} − (v^{⊥})_{L})  (4.49)

と簡単化される. ここで, (K_{p})_{V} , (K_{u})_{V} の添字 V は, 鉛直方向のセル界面に対して適用される AUSM^{+}-up スキームのフリーパラメータであることを示す. 一方, 水平方向のセル界面に対する AUSM^{+}-up スキームのフリーパラメータは添字 H を用いて, (K_{p})_{H}, (K_{u})_{H} と表すことにする.

4.4 時間積分法

水平方向の陽的な時間積分と鉛直方向の陰的な時間積分を, 時間精度を十分に維持しながら組み合わせるために, Strange-carryover 法を適用する. 水平方向の陽的な時間積分には 4 次の Runge-Kutta 法を用いる. Strange-carryover 法によって, 陽的な時間積分と陰的な時間積分の結果が組み合わされるとき, 全体として, 時間について 2 次精度が達成される.

一番最初の時間ステップにおいて, まず Δt/2 だけ陰的な時間積分を行い, q^{(1)} が求められる.

\frac{q^{(1)} − q^{n}}{(Δt/2)} − V (q^{(1)}) = 0.  (4.50)

非線形方程式 (4.50) は, ニュートン法を用いて数値的に解かれる. ニュートン法の各復計算においては, 次の形式の線形方程式を解くことになる.

[\frac{𝕀}{2Δt} − \frac{∂V}{∂q} (q^{(1),m})] (Δq^{(1)})^{m} = R(q^{(1),m}) ≡ \frac{q^{(1),m} − q^{n}}{Δt/2} + V (q^{(1),m}) (4.51)

ここで, 𝕀はランク5N_{r} の単位行列, mはニュートン反復の回数を表す添字, (Δq^{(1)})^{m} = q^{(1),m+1}−q^{(1),m} である. ニュートン法におけるヤコビアン\frac{∂V}{∂q} は, 本研究では解析的に計算される. 今の場合, 鉛直方向の数値流束の評価に凖線形化された AUSM^{+}-up スキームを用いるため, 比較的容易に\frac{∂V}{∂q} を求めることができる. より高次の再構築や複雑な境界条件が与えられる場合には, このヤコビアンを解析的に与えることは困難であるため, 数値的に与える必要がある. DECAMO では線形方程式の解法として, 直接解法あるいは Krylov 部分空間反復法の一つである一般化最小残差法 (Generalized Minimal RESidual method : GMRES 法) を用いる. なお, 今回のテスト計算では鉛直格子分割数が小さいので, 両者の計算効率の差は小さい. したがって, 前者の方法を使って第 5 章の数値実験を行った. 陽解法と陰解法の時間積分の結果を組み合わせるために, 線形の陰的 Rosenbrock 法が用いられるが, この手法においてニュートン反復に非線形収束を求めない. したがって本研究では, ニュートン反復は基本的に一度だけ行う. なお, 完全な収束解は, 反復計算をさらに数度だけ行えば得られることを確認している.

陰的な時間積分の後, その結果 q^{(1)} を使って陽的な時間積分を実行する. 陽的な時間積分として 4 次の Runge-Kutta 法用いる場合, 陽的なステップは,

\begin{aligned}q^{(2)} & = q^{(1)} + \frac{Δt}{2}H(q^{(1)}), \\ q^{(3)} & = q^{(1)} + \frac{Δt}{2}H(q^{(2)}), \\ q^{(4)} & = q^{(1)} + \frac{Δt}{2}H(q^{(3)}), \\ q^{∗} & = −\frac{1}{3}q^{(1)} + \frac{1}{3}q^{(2)} + \frac{2}{3}q^{(3)} + \frac{1}{3}q^{(4)} + \frac{Δt}{6}q^{(4)} \end{aligned}  (4.52)

となる. 最後に, 陽的な時間積分で求められた q^{∗} を用いて, 陰的な時間積分を Δt/2 だけ進める.

\frac{q^{n+1} − q^{∗}}{(Δt/2)} − V (q^{n+1}) = 0.  (4.53)

この段階で, 陰的な時間ステップによる強制として,

G = q^{n+1} − q^{∗}  (4.54)

を保持しておく. そして, 次の時間ステップの陽的な時間積分のインプットにおいて, この強制を加える. すなわち,

q^{(1)} = q^{n} + G.  (4.55)

したがって, Strange-carryover 法は, 一番最初の時間ステップ以外は, 時間ステップ一回に対し, 陰的な時間積分を一回だけ行えば良いように見える.

4.5 パッチ境界の取り扱い

パッチ境界に隣接する検査セルにおいて再構築を行うとき, そのステンシルは隣のパッチにはみ出すため特別な取り扱いが必要である. 本研究では, 鉛直方向の再構築と同様, パッチ境界の外側にゴーストセルの層を導入し, このゴーストセルは対応する隣のパッチ内の境界近傍セルと完全に一致させることにする. したがって, ゴーストセルの体積平均値には, 対応する隣のパッチ内のセルの値が与えられる. 図 4.2 は, 通常のステンシルと隣接するパッチにステンシルがはみ出す場合を図示したものである. この図において, 対象とする検査セルは赤点で, そのステンシルが含む赤色の点以外の検査セルは紫色の点で示される. ただし, ステンシルが異なるパッチにはみ出す場合には, パッチ外に対応する検査セル (ゴーストセル) は濃い青色の点で指示される.

Ullrich and Jablonowski (2012b) とは異なり, 再構築で用いられる CENO 法は非構造格子に対するものである. したがって, ゴーストセルの位置を対象パッチの格子線の延長上に配置し, 隣のパッチ上の複数の検査体積値からゴーストセル上の検査体積値にリマッピングする, 複雑な操作は必要としない. また, CENO 法ではステンシルは静的に決定されるので, ステンシルが異なるパッチにはみ出す場合でも, そのステンシルが含む異なるパッチ上の検査セルは既知である. このことは, 袖領域のデータ交換の処理を簡単にするだけでなく, パッチ境界における数値流束の一意性を恣意的な操作なしに保証することを可能にする. したがって, 立方球面格子のパッチ間でリマッピングが用いられる場合に, 有限体積法の解変数の保存性が失われる問題は, 本研究の方法では発生しない. ただし, それぞれのパッチにおいて局所座標系が導入されているので, ベクトル量は, 対象となるパッチの座標系での量に座標変換されなければならない.

図 4.2: 検査セルとそれに対応するステンシルの関係. 境界に隣合う検査セルにおいて, ステンシルは異なるパッチにはみ出す. 濃い青色の点は, 青色のパッチのゴーストセルと関連付けられる, 近隣のパッチ内の検査セルを示している. パッチの角の検査セルに対するステンシル (Corner Stencil) において, 保持される検査セルの数は他のステンシルよりも一つ少ない.

第5章数値実験

DECAMO の支配方程式系は完全圧縮系のため, 乾燥大気における全ての種類の波動が解として含まれる. したがって, DECAMO の大気循環モデルとしての数値精度や特性を検証する上で, まず初めに, 数値モデルにおける基本的な大気波動を振る舞いを調べることが必要である. 故に本章では, DECAMO の性能を検証するために, 音波のテストケース (Tomita and Satoh, 2004), 内部重力波のテストケース (Tomita and Satoh, 2004; Jablonowski et al., 2008) 傾圧不安定のテストケース (Jablonowski and Williamson, 2006) を行う.

5.1 テストケース 1: 音波

DECAMO の支配方程式系は完全圧縮系のため, モデルは音波を含む. しかし, 地球流体力学的な興味において音波が重要な場合は少ない一方で, 高周波の波は計算不安定の要因になりやすい. そのため, 多くの非静力学モデルでは divergence damping(Skamarock and Klemp, 1992) を用いることによって, 音波を選択的に減衰させて制御する. DECAMO における数値安定性の仕組みは, AUSM^{+}-up スキームに備え付けられている拡散効果によってもたらされる. この効果の強さは, フリーパラメータ (K_{p})_{H}, (K_{u})_{H}, (K_{p})_{V} , (K_{u})_{V} に依存する. 本テストケースの目的は, DECAMO における音波伝播の振る舞いを調べることと, 総観規模の流れに対する数値安定性のためのフリーパラメータの適切な値を検証することである.

この音波伝播テストにおける, 基本場と初期場の設定を示す. 基本場は T_{0} = 300K の等温大気を与える. このとき, 基本場の圧力は,

\overline{p}^{h}(z) = P_{0} exp (−\frac{z}{H} ).  (5.1)

ここで, P_{0}(= 1000hPa) は地上圧力, H = RT_{0}/g は等温大気のスケールハイトである. また, 初期条件として, 次の圧力擾乱を圧力の基本場に重ねる.

p′ = Δp f(λ, θ)g(z),  (5.2)

ここで, Δp は圧力擾乱の振幅, f(λ, θ), g(z) はそれぞれ圧力擾乱の水平方向および鉛直方向の分布関数であり, 以下のように定義する.

f(λ, θ) = \begin{cases}\frac{1}{2}(1 + cos(πl/R)) & (r < R), \\ 0 & (r > R), \end{cases}  (5.3)

g(z) = sin (\frac{n_{v}πz}{z_{T}} ).  (5.4)

ただし, R はコサインベルの半径, n_{v} は鉛直モード数, z_{T} はモデル上端高度, l はコサインベル中心 (λ_{0}, ϕ_{0}) から (λ, ϕ) までの大円距離を表す. ここでは, 圧力擾乱に関するパラメータを a = 6371 km, Δp = 100 hPa, R = a/3 km, (λ_{0}, ϕ_{0}) = (0, 0), n_{v} = 1, z_{T} = 10 km に設定した.

本ケースでは, 数値モデルにおける音波伝播の振る舞いを調べる中で, AUSM^{+}-up スキームの数値安定性と関連するフリーパラメータ (K_{p})_{H}, (K_{u})_{H}, (K_{p})_{V} , (K_{u})_{V} 依存性を考察する. そのために, 表 5.1 に示すような 4 種類のサブテストを行う.

ケース (K_{p})_{H}/(K_{p})_{ref} (K_{u})_{H}/(K_{u})_{ref} (K_{p})_{V} /(K_{p})_{ref} (K_{u})_{V} /(K_{u})_{ref}
1-1  1, 10^{−1}, 10^{−2}  1  1  1
1-2  1  1, 10^{−1}, 10^{−2}  1  1
1-3  1  1  5, 1, 10^{−1}  1
1-4  1  1  1  5, 1, 10^{−1}

表 5.1: 各サブテストにおける AUSM^{+}-up スキームの数値安定性のためのフリーパラメータの設定

各サブテスト間で共通の設定を説明する. 解像度は, 立方球面格子の水平解像度は c46(格子間隔は約 200km), 鉛直解像度は等間隔に 20 層 (Δz = 500m) とする. また, 本テストケースでは自転角速度はゼロ (f = 0) とする. 領域の上端と下端における境界条件にはスリップ条件を課した. 時間積分は 48 時間行い, 時間ステップ幅を 450 秒, 鉛直方向の陰解法におけるニュートン反復の最大反復回数を 1 回に設定した.

以下テストケース 1 の結果を示す. 図 5.1 は, テストケース 1-1 の(K_{p})_{H}/(K_{p})_{ref} = 10^{−2} に対する t = 0, 4, 8 時間後の圧力擾乱の鉛直分布と水平分布を示している. このテストにおいて圧力擾乱は, 球面上を大円として音波の速度で伝搬する. この波のモードは, 静力学的な流れのレジームにおける圧縮モードすなわちラブ波 (lamb wave) に対応している.

図 5.2 は, テストケース 1-1 における運動エネルギー, 全ポテンシャルエネルギー (ポテンシャルエネルギーと内部エネルギーの和), 全エネルギーの時間的な変化を示している. 約 16 時間後に存在する運動エネルギーと全ポテンシャルエネルギーの最初のピークは, 初期擾乱の位置の対極に擾乱が到達した時刻に対応している. また, 約 33 時間後に存在する 2 回目のピークは, 擾乱が初期擾乱の位置に戻ってきた時刻に対応している. これより, モデルにおけるラブ波の伝播速度はおおまかには 338 ms^{−1} であることが推定される. 一方, 基本場は 300K の等温大気であるので, 音速の理論値は 347 ms^{−1} と計算される. したがって, モデルの音波に対する位相誤差は最大でも数パーセントであることが分かる.

図 5.2 および図 5.3 はそれぞれ, ケース 1-1, 1-2 における各エネルギーの時間変化に加えて, それらの (K_{p})_{H}, (K_{u})_{H} 依存性を示している. これらの図より, (K_{p})_{H}, (K_{u})_{H} を Liou (2006) の値に設定した場合に, いくらか運動エネルギーが減衰することが分かる. 一方で, サブケース 1-3, 1-4 の結果から, (K_{p})_{V} , (K_{u})_{V} の値に依る各エネルギーの時間変化の違いは, ケース 1-1, 1-2 の場合に比べて非常に小さかった. 流れの場に関しては, サブケース 1-1, 1-2 において, 各フリーパラメータ間であるいはサブケース間で, 高波数のノイズの減衰の程度に若干の違いはあったもの, 卓越する構造に大きな違いは現れなかった. しかし, サブケース 1-3, 1-4 では,

図 5.1: テストケース 1-1 の (K_{p})_{H}/(K_{p})_{ref} = 10^{−1} に対する音波伝播の結果. 上段は経度 0 度における圧力擾乱の鉛直分布, 下段は z = 250 m における圧力擾乱の水平分布を表す. また, (a), (b), (c) はそれぞれ t = 0, 4, 8 時間後の圧力擾乱の鉛直・水平分布であることを示している.

図 5.2: ケース 1-1 における運動エネルギー, 全ポテンシャルエネルギー, 全エネルギーの時間的な変化およびそれらの (K_{p})_{H} 依存性. 黒の実線, 赤の破線, 黄緑の一点破線は, (K_{p})_{H} の値をそれぞれ (K_{p})_{ref}, (K_{p})_{ref} × 10^{−1}, (K_{p})_{ref} × 10^{−2} に設定したときの結果を示している. 全ての場合に対して, (K_{u})_{H}, (K_{p})_{V} , (K_{u})_{V} の値は Liou(2006) の基準値に固定している.

図 5.3: ケース 1-2 における運動エネルギー, 全ポテンシャルエネルギー, 全エネルギーの時間的な変化およびそれらの (K_{u})_{H} 依存性. 黒の実線, 赤の破線, 黄緑の一点破線は, (K_{u})_{H} の値をそれぞれ (K_{u})_{ref}, (K_{u})_{ref} × 10^{−1}, (K_{u})_{ref} × 10^{−2} に設定したときの結果を示している. 全ての場合に対して, (K_{p})_{H}, (K_{p})_{V} , (K_{u})_{V} の値は Liou(2006) の基準値に固定している.

図 5.4: サブケース 1-3, 1-4 における 48 時間後の赤道上での鉛直速度の鉛直分布. 単位は m/s.

(a) は (K_{p})_{V} /(K_{p})_{ref} = 1, (K_{u})_{V} /(K_{u})_{ref} = 1, (b) は (K_{p})_{V} /(K_{p})_{ref} = 5, (K_{u})_{V} /(K_{u})_{ref} = 1 (c) は (K_{p})_{V} /(K_{p})_{ref} = 1, (K_{u})_{V} /(K_{u})_{ref} = 10^{−1} (d) は (K_{p})_{V} /(K_{p})_{ref} = 1, (K_{u})_{V} /(K_{u})_{ref} = 5 の場合である.

(K_{p})_{V} , (K_{u})_{V} の大きさによって, 鉛直速度の分布に違いが現れた. この鉛直速度の鉛直分布に対する (K_{p})_{V} , (K_{u})_{V} 依存性は, 図 5.4 に示される. サブケース 1-3 において, 圧力に依存する拡散が大きすぎる ((K_{u})_{V} /(K_{u})_{ref} > O(1)) 場合に, 図 5.4(b) に示すように, 鉛直速度場において鉛直方向に偽の振動が見られた. これは, 圧力に関する数値安定性は非線形的に働くため, ニュートン反復が一回では拡散の効果が正しく得られないためだと考えられる. 実際, 反復回数を二回に増やすことによってこの偽の振動は消滅する. また, サブケース 1-4 では, 鉛直速度に依存する拡散が十分でない ((K_{u})_{V} /(K_{u})_{ref} < O(1)) 場合に, 図 5.4(c) に示すように, 鉛直速度の鉛直分布に激しい偽の振動が現れた. なお領域のアスペクト比が大きい場合に, 偽の振動を消去するために比較的大きな (K_{u})_{V} (すなわち, (K_{u})_{V} /(K_{u})_{ref} ≥ O(1)) が必要であることは, Urllich et al(2012) においても述べられている. ただし, 彼らが採用している支配方程式系は, 全エネルギーの時間発展式の代わりにエントロピー密度を使っている点が, DECAMO とは異なっている.

上の結果から, 以後のテストケースでは, 水平解像度 c48, 鉛直格子間隔 500 m に対しては, 数値安定性のためのパラメータを

(K_{p})_{H}/(K_{p})_{ref} = 10^{−1}, (K_{u})_{H}/(K_{u})_{ref} = 10^{−1},
(K_{p})_{V} /(K_{p})_{ref} = 1, (K_{u})_{V} /(K_{u})_{ref} = 5 (5.5)

と設定する. 異なる格子解像度に対しては, 上記のパラメータと対応する拡散時間スケールと同程度になるように, (K_{p})_{H}, (K_{u})_{H}, (K_{p})_{V} , (K_{u})_{V} の値を調節することにする.

最後に, 系の全質量・全エネルギーの保存特性について述べる. DECAMO で用いられる非静力学方程式系は完全圧縮系であり, フラックス形式で記述されているため, 数値計算法が適切に選ばれるとき, 系の質量・運動量・全エネルギーは本質的には保存される. 本テストケースでは, 図 5.2 に示すように, 系の全エネルギーが正しく保存していることが分かる. 定量的には, 系の全エネルギーは初期の値に対して 10^{−12} % の範囲で変動する一方, 図には示していない系の全質量は初期の値に対して 10^{−14} % の範囲で変動する. したがって, 系の全質量や全エネルギーは, 丸め誤差の範囲で保存することが確認された. なお, 他のサブテストでもこれらの変動量はほとんど変わらなかった.

5.2 テストケース 2: 内部重力波

本テストケースでは, 非回転系における内部重力波の伝播を調べる. このテストは, Skamarock and Klemp (1994) の 2 次元領域非静力学モデルにおけるテストを全球領域に拡張したものである. 前ケースの音波とは異なり, 内部重力波は温位の初期擾乱によって引き起こされる.

基本場は, 一定の Brunt-Väisälä 振動数 N をもつ静止大気によって構成する. このとき, 圧力の基本場は次のように与えられる.

p(z) = p_{0} [(1 − \frac{S}{T_{0}} ) + \frac{S}{T_{0}} exp (−\frac{N^{2}z}{g} )]^{c_{p}/R_{d}} .  (5.6)

ここで, p_{0} = 1000 hPa, T_{0} = 300 K, S = g^{2}/(c_{p}N^{2}) は静的安定度である. 温位の場 θ は, 静水圧バランスにある基本場\overline{θ}^{h} と擾乱場 θ′ に分割される.

θ(λ, ϕ, z) = T_{0} exp (\frac{N^{2}z}{g} ) + Δθf(λ, ϕ)g(z)  (5.7)

ここで, Δθ = 10 K は温位擾乱の最大値である. f(λ, ϕ), g(z) はそれぞれ擾乱の水平, 鉛直分布を決める関数であり, (5.3), (5.4) によって与える. ただし, コサインベルの半径は R = a/3, 鉛直モード数は n_{z} = 1, 2 とする.

本テストケースでは, 全部で 4 つのサブテストを行う. サブケース 4-1, 4-2, 4-3 では, 浮力振動数 N や鉛直モード数 n_{v} を変化させて, モデルが内部重力波の伝播や場の分布を正しく表現できるかを検証する (Tomita and Satoh (2004) のテストケース 2 に対応する). 具体的には, サブケース 4-1, 4-2, 4-3 においてそれぞれ, (N , n_{v}) = (0.01, 1), (0.01, 2), (0.01, 3) と設定する. ただし, N の単位は秒である. モデルの解像度や時間ステップはテストケース 1 と同様に, 水平解像度 c48(水平格子間隔約 200 km), 鉛直層数全 20 層, 時間ステップを 450 s とした. 一方, サブケース 4-4 では, サブケース 4-1 の水平解像度を c90(水平格子間隔約 100km), 時間ステップを 225 s に設定して計算を行う (Jablonowski et al. (2008) のテストケース 6-0-0 に対応する). そして, Jablonowski et al. (2008) のテストケース 6-0-0 に対する他の全球非静力学モデルの計算結果と比較する.

サブケース 2-1, 2-2, 2-3 における 2 日後の赤道上の温位擾乱の鉛直分布を図 5.5 に示す. 内部重力波の線形理論によれば, 一定の Brunt-Väisälä 振動数 N で成層している基本場において, 鉛直モード数 n_{v} の内部重力波の位相速度は,

c_{g} = \frac{Nz_{T}}{πn_{v}}  (5.8)

によって与えられる. ただし, 水平波数は鉛直波数に対して十分小さいことを仮定した. したがって, サブケース 2-1, 2-2, 2-3 において推定される内部重力波の位相速度は, 31.8, 63.7, 15.9 ms^{−1} と求められる. 一方, 図 5.5 から読み取れる内部重力波の位相速度は, 線形理論による推定値と比較して, サブケース 2-1, 2-2 では 10%, サブケース 2-3 では 20%ほど速い. この要因は温位擾乱の振幅が大きい (Δθ = 10K) ためであり, 振幅を小さくすれば内部重力波の位相速度は線形論の値に近づく.

サブケース 2-4 における 4 日後の赤道上の温位擾乱の鉛直分布を図 5.6 (A) の (a) に示す. また, 図 5.6 には, 地球流体電脳倶楽部で開発されている DIGMO(Dennou Icosahedral Grid MOdel)^{*1} および NCAR で開発されているコミュニティ大気モデル群^{*2} CAM-EUL, CAM-FV, CAM-SE の各力学コアを用いた場合のテスト結果 (Ullrich and Jablonowski, 2012b) を一緒に示す. DIGMO は, 筆者が Tomita and Satoh (2004) に基づいて実装した有限体積法による正二十面格子大気モデルである. また, CAM-EUL はスペクトル変換法に基づく力学コア, CAM-FV は有限体積法, セミラグラジアン法に基づく力学コア, CAM-SE はスペクトル要素法による力学コアである. CAM-EUL に関しては, 陽な拡散を加えない場合と 4 次の超粘性 (K_{4} = 5 × 10^{14}m^{4}s^{−1}) を加えた場合の両方の結果を示している. また, CAM-FV に関しては divergence damping を無効にしてあり, CAM-SE の結果は 4 次の超粘性 (K_{4} = 9.6 × 10^{14}m^{4}s^{−1}) を加えた結果を示している. DIGMO, CAM-FV, CAM-SE の解像度は約 1^{∘}× 約 1^{∘} L20, CAM-EUL はそれと同等な T106 L20 に設定されている. 図 5.6 より, DECAMO の結果は他のモデルの結果と競合的であることが分かる. DECAMO の結果は, DIGMO(図 5.6(A) の (b)), 超粘性ありの CAM-EUL(図 5.6(B) の (b)) や CAM-FV(図 5.6(B) の (c)) の結果と比較して, 拡散がより小さいこ

*1http://www.gfd-dennou.org/member/ykawai/work/IGModel.htm

*2http://www.cesm.ucar.edu/models/atm-cam/ \overline{main.tex 2013/04/02(河合佑太)}

図 5.5: 48 時間後の赤道上での温度摂動の鉛直分布. 単位は K. (a) はサブケース 2-1(N = 0.01 s^{−1}, n_{v} = 1), (b) はサブケース 2-2(N = 0.02 s^{−1}, n_{v} = 1) (c) はサブケース 2-3(N = 0.01 s^{−1}, n_{v} = 2) の結果を示している.

とが分かる. また, 拡散なしの CAM-EUL(図 5.6(B) の (a)) や CAM-SE(図 5.6(B) の (d)) に見られるように, 重力波の先端部に急勾配を維持することに成功している. 加えて, DECAMO の結果から読み取れる重力波の伝播速度は, 他のモデルと一貫性があることも確認できた. なお DECAMO の結果において, 重力波が伝播する前の領域で見られるノイズは立方格子の不規則性によって生じるものである. しかしながら, そのノイズの大きさは非常に小さく, 卓越的な構造が一度現れれば見えなくなる程度のものである.

図 5.6: 96 時間後の赤道上における温位摂動の鉛直分布と各数値モデル間での結果の比較. 単位は K. (A)-(a) は DECAMO, (A)-(b) は DIGMO(Tomita and Satoh (2004) と同様の定式化と数値手法を用いている) の結果を示している. また, (B)-(a) は CAM-EUL(拡散なし), (B)-(a) は CAM-EUL(拡散あり), (B)-(c) は CAM-FV(拡散なし), (B)-(d) は CAM-SE(拡散あり) の結果を示している. なお, (B) の図は Ullrich and Jablonowski (2012b) の図 10 を引用している.

5.3 テストケース 3: 傾圧不安定

本節では, Jablonowski and Williamson (2006); Jablonowski et al. (2008) で提案されている, 傾圧不安定のテストを行う. 本テストケースは非常に理想化されているものの, 背景場は, 実際の大気のそれとよく似た構造が選ばれる. したがって, このテストは, 大気モデルによって現実的な大気の流れを模倣するにあたっての良いベンチマークを与えるだろう.

Jablonowski and Williamson (2006) のテストケースは, 二つの部分に分けられる. 最初に, 力学コアは, 定常状態で傾度風バランスしている初期条件を使って初期化される. この初期条件は, 静力学プリミティブ方程式の解析解によって与えられる. したがって, この初期化は, 数値丸め誤差や離散化誤差, 重力波によって定常状態が破壊されることなく, 力学コアがどの程度定常状態を維持できるかを明らかにする. 次に, 相対的に大規模であるが局所的なガウシアン型の摂動を北半球中緯度の帯状流に重ね合わせる. この小振幅の擾乱は, 数日間に渡る傾圧波の成長のきっかけとなる.

Jablonowski and Williamson (2006) のテストケースはもともと, η = p/p_{s} によって定義される圧力に基づく座標系で定式化されている. そのため, η 系で記述された初期条件を高度座標系に変換する必要があるが, この座標変換は Jablonowski and Williamson (2006) の付録で説明される数値的な手法によって行った. なお, このテストケースの大気下端は等ジオポテンシャル面と一致しないため, 高度座標系を採用しているモデルに対してより厳密に同じ設定を行うためには, 地形の導入が必要となる. しかし, 本研究では, モデルの定式化において地形が考慮されていない. 一方で, Jablonowski and Williamson (2006) は鉛直座標系の違いに伴う境界条件などの差異は, 10 日後までの傾圧波の時間進展に重要でないと述べている. したがって, 本研究では, このテストケースにおける大気下端を, ジオポテンシャルがゼロの等値面とする. この大気下端の取り扱いの違いのために, Jablonowski and Williamson (2006) などで示される他のモデルの結果と細部を比較する際には注意が必要である.

次に, 傾圧不安定のテストケースにおいて, 本モデルに対して行った設定を述べる. 背景場\overline{p}^{h}, \overline{ρ}^{h} は, 擾乱を除いた初期条件によって与える. また, 定常帯状流に重ねられる摂動は (λ_{c}, ϕ_{c}) = (π/9, 2π/9) に中心を持ち, その分布は次のように与えられる.

u′(λ, ϕ, η) = u_{p} exp (− (\frac{r_{c}}{R})^{2}) .  (5.9)

ここで, R = a/10 は摂動の半径, u_{p} = 10 ms^{−1} は最大振幅, r_{c} は摂動中心からの大円距離である. 解像度は c90(水平格子間隔約 111 km), 鉛直層数 26 層, また時間刻み幅は 225 s に設定し, 10 日間の時間積分を行った. また, このテストケースは浅い大気近似を仮定にしているので, DECAMO の浅い大気近似のオプションを有効にした. なお, 大気下端に現れる境界層を少しでも表現できるように, 大気下層で解像度を高くする. このために, 鉛直格子は次のように不等間隔に配置される.

z_{k} = \frac{z_{T}}{\sqrt{b + 1} − 1} [\sqrt{bξ_{k}^{2} + 1} − 1]  (5.10)

ここで, ξ_{k} = (k − \frac{1}{2})/N_{z} (k = 1, …, N_{z}) である. また, b は上層に向かってどれくらい速く格子間隔を一定にするかを決定するパラメータである. ここでは b = 10, 大気上端高度 z_{T} を 30 km とした.

以下にテストケースの結果を示す. 初めに, 定常状態のテストケースの結果を, 図 5.7 に示す. ここでは, 初期からバランス場をどの程度維持できているかを示す指標として, 以下で定義される誤差ノルム l_{2}(p_{s}(λ, ϕ, t) − p_{0})) を用いる.

l_{2}(p_{s}(t) − p_{0}) = \frac{∑_{i,j}(p_{s}(λ_{i}, ϕ_{j}, t) − p_{0})w_{j}(ϕ_{j})}{∑_{i,j} w_{j}(ϕ_{j})}.  (5.11)

ここで, p_{0}(= 1000 hPa) は Jablonowski and Williamson (2006) の η 座標系 (圧力座標系とシグマ座標系のハイブリッド) を用いる定式化において定義される, 初期の地表面圧力である. また, w_{j} = | sin ϕ_{j+1/2} −sin ϕ_{j−1/2}| である. この指標は, Lauritzen et al. (2010) の定常状態のテストケースで用いられるものと同じである. 本研究では 5.11 を計算するために, 立方球面格子上の圧力の検査体積平均値を再構築を用いて緯度経度格子上に補間し, さらに幾何的な鉛直座標から η 座標系の地表面への補間には, 線形内挿あるいは外挿を用いる. Lauritzen et al. (2010) は, 誤差ノルム l_{2}(p_{s}(t) − p_{0}) が 2 × 10^{−1} ∼ 4 × 10^{−1} hPa に達すると偽の波が指数関数的に成長し始めるため, バランス場の維持に対する評価の閾値を l_{2}(p_{s}(t) − p_{0}) = 0.5 hPa としている. 図 5.7 において, 閾値 l_{2}(p_{s}(t) − p_{0}) = 0.5 hPa の直線を赤線で示す. この閾値に達するまでの期間は, 本モデルでは, 水平解像度 c46 に対して 8 日間, 水平解像度 c90 に対して 10 日間強ほどであることが分かる. このように, 一般に正二十面格子や立方球面格子を用いた力学コアは, 緯度経度格子によるスペクトルモデルに比べて, ずっと早く初期場の東西対称性が失われてしまう. 東西対称性の維持時間に関して, 他の立方球面格子モデルの結果と比較するためには, Jablonowski and Williamson (2006, 図 4) が参考になる. これより, 本モデルにおける東西対称性の維持時間は, 立方球面格子を用いた有限体積法による GEOS FV CUBED(水平方向に 2 次精度の再構築) に比べ同程度かやや長く, 立方球面格子を用いたスペクトル要素法による HOMEE (水平方向に 4 次精度の再構築) に比べ 1 週間ほど短いことが分かる. 本モデルの水平方向に関して再構築が 3 次精度であることを考慮すれば, この結果は整合的である.

傾圧不安定のテストケースの結果として, 4, 6, 8, 10 日目の 850 hPa 面の温度分布と地表面圧力を, それぞれ図 5.8, 図 5.9 に示す. 他の力学コアによる同テストの結果は, Jablonowski and Williamson (2006); Lauritzen et al. (2010) より得られる. これらの比較により, DECAMO の結果は地表面圧力の極値の位置を正しく捉えていることが分かる. また, 8 日目以降の温度場において, 前線付近の急な水平温度勾配についても他のモデルと同じように維持されていることが分かる. 一方で, 地表面圧力の極値の強さは 10hPa ほど他のモデルの結果よりも大きかった.

地表面圧力の極値がやや大きい傾向は, Lauritzen et al. (2010, 図16)の正二十面格子モデル ICON の結果にも見られる. しかし, 本研究での大気下端は, Jablonowski and Williamson (2006) の大気下端と同じに設定できていないことに注意が必要である. すなわち, 緯度 40 度より高緯度では, 本来の正しい大気下端が計算領域外と

図 5.7: 水平解像度 c46, c90 に対する定常状態テストにおける l_{2}(p_{s}(t) − p_{0}) の 10 日間の時間変化. バランス場の維持状態を評価するための閾値 l_{2}(p_{s}(t) − p_{0}) = 0.5 hPa の直線は, 赤線で示される.

なる (Jablonowski and Williamson, 2006, 図 1) ことが, 傾圧波の発達の予想に多少なりとも影響を与える可能性が考えられる. したがって, この問題についてより議論を進めるためには, まず地形を導入する必要であろう.

また, 8 日目あたりから, 立方球面格子の格子インプリンティングに伴う偽の波が, 南半球の傾圧帯で若干確認され始める. この特徴は, 定常状態のテストケースの結果からも推測される. 一方, 水平方向に関して 4 次精度の再構築を行った Ullrich and Jablonowski (2012b) の Mcore の同テストケースの結果において, 少なくとも 10 日目までの温度や地表面圧力の水平分布の中に格子インプリンティングの影響は見られなかった. したがって, 本モデルにおいても水平方向の空間離散化を 4 次精度まで高精度化すれば, この格子インプリンティングの問題は改善されるだろう. その可能性は, Lauritzen et al. (2010, 図 9) の水平方向に 2 次精度の再構築を行う GEOS FV CUBED の結果と, 水平方向に 4 次精度の再構築を行う HOMEE の結果の比較からも示唆される.

図 5.8: 傾圧不安定のテストケースにおける (a)4 日後 (b)6 日後 (c)8 日後 (d)10 日後の地表面圧力. ただし, この地表面は幾何的な地表面ではなく, Jablonowski and Williamson (2006) で規定される σ 座標系の σ = 1 の等値面である.

\underline{47}

{図 5.9: 傾圧不安定のテストケースにおける (a)4 日後 (b)6 日後 (c)8 日後 (d)10 日後の 850hPa \atop 面の温度分布.} \overline{main.tex 2013/04/02(河合佑太)}

第6章結論

本論文では, 空間高精度で計算効率が良く, 系の保存特性を完全に満たす, 高精度有限体積法による全球完全圧縮非静力学モデル DECAMO を作成した. 支配方程式系は, 保存形式の密度, 運動量, 全エネルギーの時間発展式によって記述される. モデルの格子には立方球面格子を採用することによって, 極の特異性を回避する. 高精度風上有限体積法で重要となる再構築には, Central ENO 法 (Ivan et al., 2012) を用いた. この方法は原理的には任意空間精度を可能にするが, 今回は水平方向に 3 次精度を保証する 2 次多項式による再構築を行った. また, この手法では立方球面格子の各パッチ間での格子線上への場のリマッピングを必要としないため, 系の保存量を正確に保存することができる. 検査セルの水平界面における数値流束の評価には, Ullrich and Jablonowski (2012b) と同様に, 風上分離法の一つである AUSM^{+}-up スキーム (Liou, 2006) を適用した. 時間積分スキームについては, 鉛直方向に伝播する位相速度の速い波による時間刻み幅の制約を回避するために, 鉛直方向の運動に関係する項は陰解法を用いた. この際, 陰解法の計算負荷を軽減するために, 鉛直方向の数値流束の評価に対しては, 準線形化した AUSM^{+}-up スキーム (Ullrich and Jablonowski, 2012a) を使った. 水平方向の運動の項は陽的に時間積分し, 両者の寄与は演算子分割 Runge-Kutta-Rosenbrock スキームによって組み合わされる. 結果的に, 全体では空間 2 次精度, 時間 2 次精度の数値モデルを構築した.

開発した数値モデルの性能を検証するために, 三種類の非静力学コアに対する標準テストケース (音波, 内部重力波, 傾圧不安定のテストケース) を実施した. テストケースから得られた結果は, 他の高精度非静力学モデルから得られた結果と一貫性があることが示された. また, 内部重力波や傾圧不安定のテストケースの結果を通して, 拡散が少なく, 場の急勾配を維持できることが確認できた. これは, 風上分離法や高精度再構築によるものだと思われる.

すぐに行うべき課題は, 水平方向の空間精度を 4 次精度まで拡張することと, その収束性を定量的に確認することである. また, より高解像度, 長時間の積分を必要とする大気の数値シミュレーションを行うためには, 現在実装済みの OpenMP によるスレッド並列化に加え, MPI(Message Passing Interface) によるプロセス並列化が必須である. この段階において, 初めて本モデルの並列計算環境におけるスケーラビリティの議論が可能となる. また, 本モデルは, プログラミング言語として C++ を用いて, オブジェクト指向に基づき設計されている. この設計は単にクラス化によるコードの再利用性や可読性を向上するだけでなく, 格子の局所的な細分化や並列プロセスの動的な割り当てによる並列処理の最適化といった, 複雑な処理を含む拡張を容易にするだろう. 最後に, 本モデルの力学コアは深い大気の運動を取り扱えるため, 今後は木星表層あるいは内部の対流活動を数値シミュレーションしたいと考えている.

謝辞

本論文の作成に当たり, 多くの方にご協力頂きました. 指導教官である林祥介教授には, 研究グループの研究計画もある中, 自分が興味を持った内容を自由に研究することを許容して下さったことに大変感謝しています. そして大変ご多忙にも関わらず, 研究に対するコメントだけでなく, 健康面にもご配慮下さいました. 基礎理論講義では, どのような疑問に対しても, 時間を割いて丁寧に説明して下さいました. また, 岩山隆寛准教授には, 普段の研究の相談はもちろん, 物理数学や共通の趣味の六甲ハイキングなど楽しいお話をして下さり, また挽きたてのコーヒーをたくさんご馳走になりました. 高橋芳幸特任助教には, 毎週の読書会において抱いた地球流体力学の疑問に対し, 一緒に向き合い答えを下さりました. そして, 日頃の研究室の活動, 進学に対する不安など様々な面でお世話になりました. かつて本研究室に所属されていた西澤誠也助教, 佐々木洋平氏には, 数値モデル開発において, 多くのきっかけや知識を与えて下さりました. NICAM 開発者の一人である富田浩文様は, 大変お忙しいにも関わらず, DECAMO の開発のきっかけとなった正二十面格子モデル DIGMO に関する相談に乗って下さいました. また, 地球および惑星大気研究室の関係者の皆さまにも大変お世話になりました. この場を借りて感謝申し上げます.

付録A {球面立方座標系の幾何的な \atop 性質}

この付録では, 深い大気の非静力学方程式の定式化において, 立方球面座標系を導入する際に必要な, メトリック項や座標変換などについてまとめる. 本付録全般に渡り, アインシュタインの縮約規則を適用する. また, 立方球面座標において, インデックスは, 各座標軸に対応して α, β, r の値をとるとする.

A.1 計量

計量は, 立方球面格子の各パッチで共通であるが, パッチ内で座標に応じて変化する. 一定の r の球面における 2 次元の計量は, 共変インデックスあるいは反変インデックスのいずれかを使って, 次のように定義される (G^{ij} = (G_{ij})^{−1} であることに注意).

 Ĝ_{ij} = \frac{r^{2}(1 + Y^{2})(1 + X^{2})}{δ^{4}} (\begin{matrix}1 + X^{2} & −XY \\ −XY & 1 + Y^{2} \end{matrix}) ,  (A.1)

 Ĝ^{ij} = \frac{δ^{2}}{r^{2}(1 + Y^{2})(1 + X^{2})} (\begin{matrix}1 + Y^{2} & XY \\ XY & 1 + X^{2} \end{matrix}) .  (A.2)

立方球面座標系の放射方向の基底ベクトルは, 球座標系と同様に, 一定の r 面に対しどこでも直交し, また正規化されている. 故に, 3 次元問題に対する完全な計量は, 放射方向の単位成分に沿って 2 次元成分に分解される. したがって,

 G_{ij} = (\begin{matrix}Ĝ_{ij} & 0 \\ 0 & 1 \end{matrix}) ,  G^{ij} = (\begin{matrix}Ĝ^{ij} & 0 \\ 0 & 1 \end{matrix}) .  (A.3)

曲線座標系において, 計量は座標系の直交性に加え基底の長さを決定する役割がある. 2 つのベクトルの内積は次のように定義される.

< v, w >= G_{ij}v^{i}w^{j},  (A.4)

ここで, v^{i}, w^{j} はそれぞれ立方球面座標系の各基底に対応する v, w の成分である. 立方球面座標系の基底を使って書かれたベクトルに対し, その大きさは上で定義された内積を介して計算される. すなわち, |v| =< v, v >.

共変計量テンソルの行列式, すなわち変換ヤコビアン,

 J = \sqrt{detG_{ij}} = \frac{r^{2}(1 + X^{2})(1 + Y^{2})}{δ^{3}},  (A.5)

は, 無限小の体積素 dV = Jdαdβdr を決定する. 同様に, G_{ij} の余因子行列の行列式 J_{α}, J_{β}, J_{r} は, 一つの座標変数を固定することで得られる等値面上の無限小面積 dA_{α} = J_{α}dβdr, dA_{β} = J_{β}dαdr, dA_{r} = J_{r}dαdβ を決定する. ただし,

J_{α} = \frac{r(1 + Y^{2})\sqrt{1 + X^{2}}}{δ^{2}} ,  (A.6)
J_{β} = \frac{r(1 + X^{2})\sqrt{1 + Y^{2}}}{δ^{2}} ,  (A.7)
J_{α} = \frac{r^{2}(1 + X^{2})(1 + X^{2})}{δ^{3}}.  (A.8)

また, 他のベクトルの演算は計量とその副産物を介して定義される. 2 つのベクトルの外積は,

(u × v)^{i} = JG^{ij}ϵ_{jkl}u^{k}v^{l}  (A.9)

と定義される. ここで, ϵ_{ijk} はレヴィ=チヴィタ記号である.

A.2 球面立方座標系における微分演算子

立方球面座標系の自然基底 {e_{α}, e_{β}, e_{r}} を使って, 位置ベクトル r を,

r = x^{α}e_{α} + x^{β}e_{β} + x^{r}e_{r}  (A.10)

と定義する. このとき, 勾配演算子は,

∇^{i}ϕ = G^{ij} \frac{∂ϕ}{∂x^{j}} ,  (A.11)

また発散演算子は,

∇ · F = \frac{1}{J} \frac{∂}{∂x^{k}} (JF^{k})  (A.12)

の形式をとる. ただし, 計量 G^{ij}, 変換ヤコビアン J には, それぞれ (A.3), (A.5) を用いる.

また, 2 階テンソルの発散は,

∇_{j}T^{ij} = \frac{1}{J} \frac{∂}{∂x^{j}} (JT^{ij}) + Γ_{jk}^{i}T^{jk}  (A.13)

と書かれる. ここで, Γ_{jk}^{i} は第二種クリストッフェル記号であり, 一般に計量を使って,

Γ_{jk}^{i} = \frac{1}{2}G^{im} (\frac{∂G_{jm}}{∂x^{k}} + \frac{∂G_{km}}{∂x^{j}} − \frac{∂G_{jk}}{∂x^{m}} )  (A.14)

と書かれる. 今, (A.3) を用いることにより, 等角立方球面座標系に対する Γ_{jk}^{i} は次のように得られる.

Γ^{α} = (\begin{matrix}\frac{2XY^{2}}{δ^{2}} & \frac{−Y (1+Y^{2})}{δ^{2}} & \frac{1}{r} \\ \frac{−Y (1+Y^{2})}{δ^{2}} & 0 & 0 \\ \frac{1}{r} & 0 & 0 \end{matrix}) ,  (A.15)

Γ^{β} = (\begin{matrix}0 & \frac{−X(1+X^{2})}{δ^{2}} & 0 \\ \frac{−X(1+X^{2})}{δ^{2}} & \frac{2X^{2}Y}{δ^{2}} & \frac{1}{r} \\ 0 & \frac{1}{r} & 0 \end{matrix}) ,  (A.16)

Γ^{r} = \frac{r(1 + X^{2})(1 + Y^{2})}{δ^{4}} (\begin{matrix}−(1 + X^{2}) & XY & 0 \\ XY & −(1 + Y^{2}) & 0 \\ 0 & 0 & 0 \end{matrix}) .  (A.17)

A.3 地理座標系と球面立方座標系間の座標変換

地理座標系と球面立方座標系間における任意のベクトル v の座標変換を考える. 地理座標系および球面立方座標系の自然基底をそれぞれ, {e_{λ}, e_{ϕ}}, {e_{α}, e_{β}} とするとき, ベクトル v はこの 2 種類の基底を使って,

v = v^{λ}e_{λ} + v^{ϕ}e_{ϕ} = v^{α}e_{α} + v^{β}e_{β}  (A.18)

と書かれる. ２つの座標系間でのベクトル成分の変換は, 変数行列 A を使ってなされる. すなわち,

(\begin{matrix}v^{α} \\ v^{β} \end{matrix}) = A (\begin{matrix}v^{λ} \\ v^{ϕ} \end{matrix}) , (\begin{matrix}v^{λ} \\ v^{ϕ} \end{matrix}) = A^{−1} (\begin{matrix}v^{α} \\ v^{β} \end{matrix}) .  (A.19)

ここで, λ は経度 (λ = 0 がパッチ 1 の α = 0 と一致するするように選ばれる), ϕ は緯度である. 今, 地理座標系の基底 {e_{λ}, e_{ϕ}} は正規直交基底でなく, 自然基底であることに注意されたい.

A.3.1 赤道パッチにおける座標変換

赤道パッチは, インデックス n_{p} ∈ 1, 2, 3, 4 によって表される. 今, n_{p} = 1 は λ = 0 の緯線を含むパッチに対応するとする. 各赤道パッチにおける立方球面座標系と地理座標系間の位置座標の変換は, 次の関係によって定義される.

α = λ − \frac{π}{2} (n_{p} − 1), β = arctan (\frac{tan ϕ}{cos λ})  (A.20)
λ = α + \frac{π}{2} (n_{p} − 1), ϕ = arctan (tan β cos α)  (A.21)

赤道パッチにおいて, 変換行列 A は心射座標系 (X, Y) を使って,

 A = ( 1  0_{\frac{XY}{1+Y^{2}}  \frac{δ^{2}}{(1+Y^{2})\sqrt{1+X^{2}}}} )  (A.22)
 A^{−1} = ( 1  0_{\frac{−XY \sqrt{1+X^{2}}}{δ^{2}}  \frac{(1+Y^{2})\sqrt{1+X^{2}}}{δ^{2}}} )  (A.23)

と定義される.

A.3.2 極パッチにおける座標変換

極パッチは, インデックス n_{p} ∈ 5, 6 によって表される. 今, n_{p} = 5 は北極を含むパッチ, n_{p} = 6 は南極を含むパッチに対応するとする. ここで, パッチの指示子 s を導入する. n_{p} = 5 に対して s = 1, n_{p} = 6 に対して s = −1 とする.

極パッチにおける立方球面座標系と地理座標系間の位置座標の変換は, 次の関係によって定義される.

α = s arctan(cot ϕ sin λ), β = − arctan(cot ϕ cos λ)  (A.24)
λ = − arctan (\frac{tan α}{tan β}) , ϕ = s arccot(\sqrt{tan^{2} α + tan^{2} β})  (A.25)

赤道パッチにおいて, 変換行列 A は,

 A = (\begin{matrix}\frac{−sY}{1+X^{2}} & \frac{−sδ^{2}X}{(1+X^{2})\sqrt{X^{2}+Y^{2}}} \\ \frac{sX}{1+Y^{2}} & \frac{−sδ^{2}Y}{(1+Y^{2})\sqrt{X^{2}+Y^{2}}} \end{matrix} )  (A.26)
 A^{−1} = (\begin{matrix}−\frac{sY (1+X^{2})}{X^{2}+Y^{2}} & \frac{sX(1+Y^{2})}{X^{2}+Y^{2}} \\ −\frac{sX(1+X^{2})}{δ^{2}\sqrt{X^{2}+Y^{2}}} & −\frac{sY (1+Y^{2})}{δ^{2}\sqrt{X^{2}+Y^{2}}} \end{matrix} )  (A.27)

と定義される.

\overline{main.tex 2013/04/02(河合佑太)}

付録B {幾何モーメントの計算法に \atop ついて}

第 4 章で示した k-exact 法に基づく再構築において, 2 種類の幾何モーメントが現れる. (4.14) によって与えられる幾何モーメントは, 検査セルの中心が積分領域の一部であるため, 各検査セルに対して局所的である. 一方で, (4.18) の第一項目で定義される幾何モーメントは, その積分は一つの近傍検査セルにおいてなされ, 検査セル中心は異なるセルに含まれるために, 近傍セル間の関係を表す.

いずれの幾何モーメントも完全に検査セルの幾何のみに依存するため, あらかじめ計算しておき, 後に再利用するために保持しておくことは理論的には可能である. しかし, 幾何モーメントの総数は非常に多いため, 大規模計算において全ての幾何モーメントをメモリ上に保持できない場合も考えられる. たとえば, k = 2, N_{n} = 8 の場合, 幾何モーメントの総数は各検査セルに対し 54 個となる (より高次の再構築ではこの数はずっと大きくなる). ただし, 最終的に必要となるのは擬逆行列 𝔸^{−1} なので, メモリ容量が許容されるならば, 幾何モーメントの代わりに擬逆行列を保持しておけば, 非常に効率的に再構築を行うことができる. そのため, 本研究では前処理として擬逆行列を計算しておき, メモリに保持する方法を取っている.

しかし, メモリ容量の制限から擬逆行列 𝔸^{−1} を保持できない場合には, 各時間ステップの各検査セルに対する再構築において, 係数行列 𝔸 を求めるために全ての幾何モーメントを再計算しなければならず. この計算コストは非常に高い. また, 前処理の段階で擬逆行列を保持できる場合でも, 幾何モーメントの計算の負荷は大きいため, 幾何モーメントを効率的に計算する方法が必要である.

Ivan (2010) は, これらの問題を解決するために, 各検査セルに局所的な幾何モーメントのみメモリ上に保持し, 残りの幾何モーメントは幾何的な変換によって計算する方法を提案している. したがって, (4.14) で与えられる各検査セルに対する局所的な幾何モーメントの計算は前処理で行われ, その結果はメモリ上に保持される. 一方, (4.18) の第一項目で定義される幾何モーメントは, 局所的な幾何モーメ

\underline{56}

ントを用いて, 次のように計算される.

\begin{aligned}({ \atop α^{p_{1}}β^{p_{2}}})_{IJ} & = (\frac{1}{|Z_{J}|} ∫_{Z_{J}} (α − α_{I})^{p_{1}}(β − β_{I})^{p_{2}} dV ) − \overline{(α^{p_{1}}β^{p_{2}})}_{I} \\ & = ∑_{l=0}^{p_{2}} ∑_{q=0}^{p_{1}} [_{p_{1}}C_{q} _{p_{2}}C_{l} (Δα_{IJ})^{q}(Δβ_{IJ})^{l} \overline{(α^{p_{1}−q}β^{p_{2}−l})}_{J}] − \overline{(α^{p_{1}}β^{p_{2}})}_{I}. \end{aligned} ; (B.1)

ここで, _{p_{1}}C_{q}, _{p_{2}}C_{l} は 2 項係数であり, 次の漸化式によって効率的に求められる.

nC_{r} = \frac{n − r + 1}{r} _{n}C_{r−1}; _{n}C_{0} = 1.  (B.2)

また, 2 つの近傍セル間の各方向の距離を Δα_{IJ} = α_{J} − α_{I}, Δβ_{IJ} = β_{J} − β_{I} と表した.

参考文献

Adcroft, A., C. Hill, and J. Marshall, 2004: Implementation of an atmosphereocean general circulation model on expanded spherical cube. Mon. Weather Rev., 132, 2845–2863.

Davies, T., M. J. P. Cullen, A. J. Malcolm, M. H. Mawson, A. Staniforth, A. A. White, and N. Wood, 2005: A new dynamical core for the Met Office’s global and regional modelling of the atmosphere. Q. J. R. Meteorol. Soc., 131, 1759– 1782.

Giraldo, F. and M. Restelli, 2008: A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases. J. Comput. Phys., 227, 3849 – 3877.

Ivan, L., 2010: Development of High-Order CENO Finite-Volume Schemes with Block-Based Adaptive Mesh Refinement. Ph.D. thesis, University of Toronto.

Ivan, L., A. Susanto, H. D. Sterck, and C. Groth, 2012: High-Order Central ENO Finite-Volume Scheme for MHD on Three-Dimensional Cubed-Sphere Grids.

Jablonowski, C., P.H.Lauritzen, M.A.Taylor, and R.D.Nair, 2008: Idealized test cases for the dynamical cores of atmospheric general circulation models: a proposal for the NCAR ASP 2008 summer colloquium.

Jablonowski, C. and D. L. Williamson, 2006: A baroclinic instability test case for atmospheric model dynamical cores. Q. J. R. Meteorol. Soc., 132, 2943–2975.

Lauritzen, P. H., C. Jablonowski, M. A. Taylor, and R. D. Nair, 2010: Rotated versions of the jablonowski steady-state and baroclinic wave test cases: A dynamical core intercomparison. Journal of Advances in Modeling Earth Systems, 2, n/a–n/a.

Liou, M.-S., 1996: A Sequel to Ausm: AUSM+. J. Comput. Phys., 129, 364–382.

— 2006: A sequel to AUSM, Part II: AUSM^{+}-up for all speeds. J. Comput. Phys., 214, 137 – 170.

Liu, X.-D., S. Osher, and T. Chan, 1994: Weighted essentially non-oscillatory schemes. J. Comput. Phys., 115, 200–212.

Putman, W. M. and S.-J. Lin, 2007: Finite-volume transport on various cubedsphere grids. J. Comput. Phys., 227, 55–78.

— 2008: A Finite-Volume Dynamical Core on the Cubed-Sphere Grid. Numerical Modeling of Space Plasma Flows: ASTRONUM-2008 ASP Conference Series, Astronomical Society of the Pacific, volume 406, 268.

Roe, P., 1981: Approximate riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43, 357–372.

Ronchi, C., R. Iacono, and P. Paolucci, 1996: The “ cubed sphere ” : A new method for the solution of partial differential equations in spherical geometry. J. Comput. Phys., 124, 93–114.

Sadourny, R., 1972: Conservative finite-difference approximations of the primitive equations on quasi-uniform spherical grids. Mon. Weather Rev., 100, 136–144.

Satoh, M., 2002: Conservative Scheme for the Compressible Nonhydrostatic Models with the Horizontally Explicit and Vertically Implicit Time Integration Scheme. Mon. Weather Rev., 130, 1227–1245.

Skamarock, W. and J. Klemp, 1994: Efficiency and Accuracy of the Klemp-Wilhelmson Time-Splitting technique. Mon. Weather Rev., 122, 2623–2630.

Skamarock, W. C. and J. B. Klemp, 2008: A time-split nonhydrostatic atmospheric model for weather research and forecasting applications. J. Comput. Phys., 227, 3465–3485.

Taylor, M. A., J. Edwards, and A. S. Cyr, 2008: Petascale atmospheric models for the community climate system model: new developments and evaluation of scalable dynamical cores. Journal of Physics: Conference Series, 125, 012023.

Tomita, H. and M. Satoh, 2004: A new dynamical framework of nonhydrostatic global model using the icosahedral grid. Fluid Dynamics Research, 34, 357 – 400.

Ullrich, P. A. and C. Jablonowski, 2012a: Operator-split Runge-Kutta-Rosenbrock (RKR) methods for non-hydrostatic atmospheric models. Mon. Weather Rev., 140, 1257 – 1284.

— 2012b: Mcore: A non-hydrostatic atmospheric dynamical core utilizing highorder finite-volume methods. J. Comput. Phys., 231, 5078 – 5108.

Ullrich, P. A., C. Jablonowski, and B. van Leer, 2010: High-order finite-volume methods for the shallow-water equations on the sphere. J. Comput. Phys., 229, 6104 – 6134.

Wedi, N. P., K. Yessad, and A. Untch, 2009: The Nonhydrostatic Global IFS/ARPEGE: Model Formulation and Testing, volume 594 of ECMWF technical memorandum. European Centre for Medium-Range Weather Forecasts.

Weiss, J. M. and W. A. Smith., 1995: Preconditioning applied to variable and constant density flow. AIAA Journal, 33, 2050 – 2057.

Williamson, D. L., J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, 1992: A standard test set for numerical approximations to the shallow water equations in spherical geometry. J. Comput. Phys., 102, 211–224.