地球環境研究，Vol.17（2015）

二次元非発散渦度方程式の差分形について
吉﨑

正

憲＊

キーワード：過度方程式、差分形、保存則、有限要素法、Galerkin 法

１ ．はじめにアプローチとして、有限要素法（finite element method;FEM）を説明する。FEM では一連の基底関数を定義す

時間発展する非線型な微分方程式系を解くのは非常にるので形式上は差分形と思われるが、微分の特性を持つ難しい。それでも解こうとするには、数値的に解くのがのが特徴である。5節では、FEM のこの特性を用いる普通である。その場合、モデルの数式をグリッド化することにより導出される差分式が Arawaka スキームそのことにより差分形で表現する。しばしばその微分方程式ものになることを示す。6節はまとめである。

系はいくつかの保存式を含むことがある。できれば差分形でもそうした特性を持つのが望ましい。ところが、複数の保存式を同時に満たす差分形を導くのはなかなか難しい。ここでは、これをどうクリアするかという気象の問題からの解説である。

2 ．二次元非発散渦度方程式の特性
二次元平面を考えて、ｔを時間、x, y を空間座標、ψを流線関数とすると、速度 u, v はある量を保存する微分方程式系で、これらを保存しな（1）

いような差分形、例えば、二次元非発散渦度方程式（後述）の場合には、時間に関して長時間積分を試みると、と表される。ここで xy 平面に直交する渦度成分ζは数値不安定を起こしてしまう。ここで扱う二次元非発散（2)渦度方程式は大規模な気象学では基本となるもので、気候研究のように長期の時間積分を必要とする場合には、と定義される。その運動が移流だけで決まるとすると、時間と共に計算が発散するのは深刻な問題である。

支配する非線型方程式はある式を導出するために、論理的に積み上げる試みを通してより良いものにたどり着く方法と、もっと原理的なところから出発して式を導くという方法がある。同じ（3）式ならばできれば原理的なことを知ってからの導出がスマートで効率的であるけれども、学問の進展は必ずしもと書くことができる。この式を変形すると、そういう風になっていない。そうした例が、先に挙げた（3）´

二次元非発散渦度方程式の差分形の導出である。本稿は、前者のやり方として Arakawa （1960）（差分によるアプ（3）̋

ローチ）、後者のやり方として Jespersen（1974）（Galerkin 法によるアプローチ）の二つの論文の紹介という形をとりながら、二つのやり方を眺めてみる。
本稿の構成は以下のとおりである。2節では、二次元非発散渦度方程式の特性を述べる。3節では、論理的に差分形を積み上げてたどり着いた Arakawa スキームを紹介する。４節では、保存量を満たすような方向からの＊
立正大学 ・ 地球環境科学部 ・ 環境システム学科とも書き換えられる。
ここで、境界条件を周期条件あるいは0とおき領域積分することにより、以下に示すような保存が成り立つ。
全渦度の保存：（4）101

二次元非発散渦度方程式の差分形について（吉﨑）全エンストロフィの保存：

3 ．二次元非発散渦度方程式の差分形の導出（5）ψとςは同じグリッド点にあるとして、これらを変数A と書くと、差分点（i，j）の A は A i , j と表される（第全運動エネルギーの保存：（6）２図）。これを使って、例えば、x 方向の差分と平均を、これらは二次元非発散渦度方程式が持つ特性である。
二次元非発散渦度方程式で支配される系を非線型で時と定義すると、（3）、（3）́、（3）˝ の差分形（それぞれＦ、G、H と書間発展させる場合、全運動エネルギーと全エンストロフィが保存するためにどのような制約があるだろうか。簡単となる。

く）はのために、ある時間では一つのエネルギー E 0 をもつ塊を（7）考え、次の時間にこれが二つに分裂するとしよう（それ（8）ぞれのエネルギーを E 1 と E 2 とする）
。運動エネルギーの保存から E 0 =E １ ＋E ２ は満たす。また波数空間ではエンス（9）

トロフィはエネルギーに二乗の波数 k 2 がかかる量で表され、k i 2 E i と定義できるので、全エンストロフィも保存をする場合には、分裂後のエンストロフィはそれぞれ
ｋ １ 2 E １ と k ２ 2 E ２ となり、k 0 2 E 0 =k １ 2 E １ ＋k 2 2 E ２ が成り立となる。例えば、（0,0）では（7）´́つ。これら二つを満たすアナロジーとして、重力が働く中で質量のない棒を支点 O で、それにかかる物体を支えている状況を考える（第1図）。ここで E を質量、波数k 2 は O からの距離とすると、第1a 図では k 0 2 の距離にある E 0 の質量を F 0 という力で支えていることになる。次に、質量を二つに分けて距離を変えて再配分することを（8）´考えるが、そのとき上向きの力は変わらないとする。その場合、第1b 図では、質量とトルクの保存から、二つの質量のうち大きい（小さい）質量は距離が小さい（大きい）位置に来ざるを得ない。このように二つの保存量を持つ場合には勝手に距離を選べないというきびしい拘束条件がつく。

（9）´と書くことができる。このように、移流項の差分式は３通り書くことが可能となる。

第 1 図：質量とトルクの保存を表す力学モデル
点 O は支点であり質量のない棒を F ０ という力で支える。
（a）一つの質量の場合。
（b）二つの質量に分けた場
第 1 図 質量とトルクの保存を表す力学モデル．点
O は支点であり質量のない棒を F という力で支える．(a)一つの質量の場合．(b)二つの質量に分けた場合で，同じ F という力を保ったままで二つに合で、同じF ０ という力を保ったままで二つに分かれた物分かれた物体を支えることを考える．
体を支えることを考える。
００

102

9

第 2 図：（０,０）を中心とする格子点
第 2 図 (０,０)を中心とする格子点．(i,j)は省略してある．
（i,j）は省略してある。

地球環境研究，Vol.17（2015）

要な基礎知識となるので重要なステップとなる。FEM が一方、2変数関数ξ（x,y）の重積分の差分形を（10）微分の特性を保持するということを理解しているのであれば、この節はスキップして5節に進んでもかまわない。
まず領域と定義すると、
（7）
́、（8）
́、（9）
́ を用いた式（4）、（5）、（6）に関する保存性は第１表のようになる。これから、の1次元の空間変数 x によって支配される物理量 u（x）を考える。この系は差分方程式は必ずしもすべて保存するわけではないこと（12）がわかる。また、（7）、（8）、（9）の単独あるいは二つをとかける。ここでΛは x に関する線型微分演算子、f は強組み合わせるようなスキームでは、三つの保存式を満足制項とすると、Λと f（x）は既知として u（x）を求めるすることはできないことも分かる。ところが、二つの差のが問題となる。領域内で N 個のお互い直交する関数分形を組み合わせると×と×がキャンセルしあって○とφ（x）を定義して、iなる可能性がある。実際、三つの差分形を1/3として合（13）算することにより一つの差分形を作ると、三つの保存式を満足するものができる。これが Arakawa（1960）スと近似できるとする。ここでφ（x）は基底関数と呼ばれiキームである。ここで、改めて（0,0）における Arawakaる。これによる誤差をスキームの移流部分を J 0 , 0 とすると、（14）と定義すると、誤差はそれぞれの基底関数に直交すると仮定する。つまり、（15）（11）となる。

非発散渦度方程式の差分形として、Arakawa（1960）が成り立つ。そうすると、この系はスキームが提案されるまではさまざまな差分形が提案さ（16）れた。例えば、Miyakoda（1962）は、（7）のような差分形を用いて時間に関して中央差分で近似して時間積分をと、N 個の ui に関する代数式に帰着することができる。

行った。しかし、格子間隔スケールの大きさの渦が伸びこうして、微分方程式の階数を減ずることができる。

る“noodling”と呼ばれる現象により計算不安定を起こ実際に簡単な例として、した。当時としては、この簡単な微分方程式を差分形で（17）如何に長時間積分を行うかシビアな問題であった。こうした話題は新田ほか（1972）に解説があるので、関心のある方はそちらを参照してほしい。

で u（0）= u（π）= 0を満たす場合を考える。（12）のΛはにあたる。

二通りの基底関数を使ってみる。一つ目の解法として、

4 ．有限要素法（FEM）の概要（18）

二次元非発散渦度方程式からしばし離れて、Haltinerand Williams（1980）にしたがって、Galerkin 法の概要なる基底関数を考える。この場合は境界条件を満たしている。これから、を紹介する。本節は Jespersen（1974）を眺めるために必第 1 表： 差分形に関する保存則○は保存、×は非保存を表す。
渦度（19）と書くことができる。この uk を用いて、（13）のような解が得られる。

エンストロフィー 運動エネルギー

F○××G○○×H○×○

二つ目の解法として、図2a に示すような基底関数
φ（x）を用いる。つまり、
i（20）

103

二次元非発散渦度方程式の差分形について（吉﨑）を用いると、変形により、（16）から、（21）（27）となる。
（3）の右辺の計算は煩雑であるので、（0,0）にが得られる。右辺の一階微分は、関して、ς 1 , 1 がある項だけを示す。この場合、計算を必（22）

要とする領域は格子点（1,1）
、（1,0）、（0,0）、（0,1）で囲まれる部分であり、ψ 1 , 0 、ψ 1 , 1 、ψ 0 , 1 、ψ 0 , 0 だけの点を考慮する必要がある。

であるから、（23）が得られる。同様に、f（x）も基底関数で
（24）

（28）

から、第1項は、と展開して、整理すると、（25）（29）となり、第２項はが得られる。
（19）と（25）は一見異なって見える。これは、異なる（30）基底関数を使ったことによる解の表現の違いであって、（13）としてはどちらも同じとなる。このように、いろいとなる。これから、ろな基底関数を用いることにより、誤差が最小になるよ（31）

うに微分の階数を減らし、代数式として解くやり方がGalerkin 法である。Galerkin 法が優れている点は、基本的に微分方程式をもとにしているために、微分方程式で見られる保存式がそのまま成り立つことである。大循環数値モデルにおいて、経度方向にはフーリエ展開、緯度が求まる。他の項についても同様であり、複雑な計算の結果、（11）が得られる。時間項を単にだけとすると、この式は Arakawa スキームとなる。

方向にはルジャンドル展開という方法が良く使われるが、ちなみに、Jespersen（1974）の論文のタイトルを和訳

このスペルトル法も Galerkin 法の一つの形態であり、すると「荒川スキームは有限要素法である」とある。

Galerkin 法が持つ保存性を良く保持している。

Arakawa（1960）スキームは基底関数として矩形を用いたことになるが、Jespersen（1974）はさらに三角形 a’+５ ． Galerkin 法の 2 次元非発散渦度方程式への応用（3）の問題に Galarkin 法を適用してみる。Jespersenb’x+c’y（a’, b’, c’ は定数）など様々な基底関数を使うことにより、全渦度、全エンストロフィ、全運動エネルギーの保存を満たす差分スキームを一般化できることを示した。

（1974）にしたがって、（26）

６ ．まとめ
全渦度、全エンストロフィ、全運動エネルギーが保存という矩形の基底関数を仮定する。この場合は a＋bx＋するように、二次元非発散渦度方程式の差分形を二通りcy＋dxy（a, b, c, d は定数）の基底関数を使うことになのやり方で示した。一つは、論理的な思考から導く方法

る。第1図の（0,0）で（3）の時間項は、（Arakawa, 1960）であり、二つは Galerkin 法により一つの基底関数を用いて導出する方法（Jespersen, 1974）である。Arakawa スキームは数値モデルの開発における一つのブレークスルーをもたらしたもので、気象学の発展

104地球環境研究，Vol.17（2015）

参考文献

Φi (x)

Arakawa, A., 1960: Computational design for long-term
numerical integration of the equation of fluid motion: twodimensional incompressible flow. Part 1, J. Comp. Phys., 1,1

119-143.
Haltiner, G. J. and R. T. Williams, 1980: Numerical prediction
and dynamic meteorology. Second Edition. John Wiley &0X(i -1)ΔxiΔx(i+1)ΔxSons, 477pp.
Jespersen, D. C., 1974: Arakawa’s method is a finite-element
method. J. Comp. Phys., 16, 383-390.

第 3 図： １ 次元の不連続な基底関数
第 3 図： １次元の不連続な基底関数．

においてエポックメーキング的な出来事であった。これまで気象力学の発展に際してさまざまな危機に遭遇したが、こうした努力の積み重ねによって克服してきたので

Miyakoda, K., 1962： A trial of 500 hour barotropic forecast.
Proc. Int. Symp. Numerical Weather Prediction, Tokyo,Meteor. Soc. Japan, Tokyo, pp221-240.
新田尚 ・ 遠藤昌弘 ・ 大林智徳 ・ 菊池幸雄 ・ 近藤弘輝 ・ 岩嶋樹也，1972：気象力学に用いられる数値計算法 . 気象研究
ノート，110号 .

ある。本稿により気象学の発展を振り返るのも良いと考える。

Difference forms of 2-dimensional non-divergent vorticity equationYOSHIZAKI Masanori ＊＊Faculty of Geo-environmental Science, Rissho University
11105