2 次元非静力学モデルの離散化

山下達也, 杉山耕一朗, 石渡正樹

2010 年 3 月 17 日

要旨

理想気体を仮定した乾燥大気に対する 2 次元準圧縮性方程式系の離散化を行う. 変数の格子配置は Lorenz グリッドにしたがう. 空間方向の離散化は 2 次精度中心差分を用いて行う. 時間方向の離散化は時間分割法を用いて行う. 運動方程式と圧力方程式は短い時間刻み Δτ で時間積分を行う. 音波に関する項の離散化には HE-VI 法を採用し, u の式は前進差分, w, π の式は後退差分で評価する. 音波にかかわらない項についてはリープフロッグ法を用いて積分する. 熱力学の式とその他のトレーサの式は, リープフロッグ法を用いて長い時間刻み Δt で時間積分を行う.

\underline{2 次元非静力学モデルの離散化 1}

目次

第 1 章数値計算の概要 2 1.1 格子と変数の配置 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 空間, 時間方向の離散化の方法 . . . . . . . . . . . . . . . . . . . . . 2

第 2 章空間方向の離散化 4 2.1 平均操作 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 空間微分の離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2.1 2 次精度中心差分 . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.2 4 次精度中心差分 . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 空間離散化した基礎方程式 . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.1 静水圧の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.2 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.3 圧力方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.4 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.5 比湿の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.4 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.1 周期境界条件の与え方 . . . . . . . . . . . . . . . . . . . . . 10 2.4.2 すべりなし条件の与え方 . . . . . . . . . . . . . . . . . . . . 10 2.4.3 応力なし条件の与え方 . . . . . . . . . . . . . . . . . . . . . 11

第 3 章時間方向の離散化 12

3.1 運動方程式と圧力方程式 . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.1 音波に関連する項の時間方向の離散化 . . . . . . . . . . . . 13 3.1.2 音波に関連しない項の時間方向の離散化 . . . . . . . . . . . 18

3.2 熱力学の式と混合比の保存式の離散化 . . . . . . . . . . . . . . . . 19

3.2.1 湿潤飽和調節法 . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3 乱流運動エネルギーの式 . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4 時間フィルター . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.5 スポンジ層 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

第 4 章参考文献 27 付録 A 圧力方程式 (??) の左辺の空間微分の書き下し 28

A.1 下部境界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.2 上部境界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

付録 B 音波減衰項について 31

付録 C 差分式の導出と誤差 33 C.1 2 次精度中心差分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 C.2 4 次精度中心差分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

付録 D 化学物性値の計算法 36

D.1 相平衡条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 D.1.1 飽和蒸気圧 . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 D.1.2 圧平衡定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

D.2 生成のエンタルピー変化 . . . . . . . . . . . . . . . . . . . . . . . . 37 D.2.1 潜熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 D.2.2 反応熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

付録 E 単位の換算等の計算 39 E.1 混合比とモル比 (分圧) との換算式 . . . . . . . . . . . . . . . . . . . 39

謝辞 \underline{2 次元非静力学モデルの離散化第 1 章数値計算の概要 3}

第1章数値計算の概要

1.1 格子と変数の配置

本モデルでは図 1.1 に示すように, 水平・鉛直方向の変数をずらして互い違いに配置する. 配置方法は, 水平方向に Arakawa C グリッド, 鉛直方向には Lorenz グリッドを用いている.

すべてのスカラー量 (ϕ : 気圧偏差, 温位偏差, 水蒸気混合比, 雲水混合比, 雨水混合比) を格子の中心に配置する. ベクトル量 (速度 u, w など) は中心から次のように半格子ずらして配置する. x 方向ベクトル量を評価する点は, 鉛直方向にはスカラー量と同じとし水平方向に半格子ずらす. x 方向ベクトル量を評価する点は, 水平方向にはスカラー量と同じとし鉛直方向に半格子ずらす. 格子点のラベルづけおよび境界は図 1.2 のように設定する.

1.2 空間, 時間方向の離散化の方法

空間方向の離散化は 2 次精度中心差分を用いて行う. 時間方向の離散化は時間分割法を用いて行う. 運動方程式と圧力方程式は短い時間刻み Δτ で時間積分を行う. 音波に関する項の離散化には HE-VI 法を採用し, u の式は前進差分, w, π の式は後退差分で評価する. 音波にかかわらない項についてはリープフロッグ法を用いて積分する. 熱力学の式とその他のトレーサの式は, リープフロッグ法を用いて長い時間刻み Δt で時間積分を行う.

\underline{2 次元非静力学モデルの離散化第 1 章数値計算の概要 4}

図 1.1: 格子点の配置. 図 1.2: 添字と境界の設定. 実際に計算する添字の範囲は 1 から im とする.

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化第 2 章空間方向の離散化 5}

第2章空間方向の離散化

この節では空間微分の離散化の方法とそのために必要となる平均操作, 境界条件の与え方について説明する. 離散化は 2 次精度差分または 4 次精度差分を用いて行う.

図 1.1 の空間の位置を表す添字として, x 方向フラックスの格子点を (i(u), k), z 方向フラックスの格子点を (i, k(w)), スカラー量の格子点を (i, k), 格子の角に当たる点を (i(u), k(w)) とする (図 1.1 参照). 但し 1 ≤ i(u) ≤ im, 1 ≤ i ≤ im, 1 ≤ k(w) ≤ km, 1 ≤ k ≤ km である.

2.1 平均操作

空間微分の離散化を行う前に, そのために必要となる平均操作を定義しておく. 例えば x 方向フラックス格子点で評価される変数をスカラー量の格子点で評価する場合は, フラックス格子点の値を平均してスカラー格子点での値とみなす.

必要となる平均操作を以下に示す. ここでは x 方向のフラックス格子点の変数を u_{i(u),k}, z 方向のフラックス格子点の変数を w_{i,k(w)}, スカラー格子点の変数を Π_{i,k} としている.

Π_{i(u),k} ≡ \frac{Π_{i+1,k} + Π_{i,k}}{2}  (2.1)
Π_{i,k(w)} ≡ \frac{Π_{i,k+1} + Π_{i,k}}{2}  (2.2)
Π_{i(u),k(w)} ≡ \frac{Π_{i,k} + Π_{i+1,k} + Π_{i,k+1} + Π_{i+1,k+1}}{4}  (2.3)
u_{i,k} ≡ \frac{u_{i(u),k} + u_{i−1(u),k}}{2}  (2.4)
u_{i,k(w)} ≡ \frac{u_{i(u),k+1} + u_{i−1(u),k+1} + u_{i(u),k} + u_{i−1(u),k}}{4}  (2.5)
u_{i(u),k(w)} ≡ \frac{u_{i(u),k+1} + u_{i(u),k}}{2}  (2.6)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

w_{i,k} ≡ \frac{w_{i,k(w)} + w_{i,k−1(w)}}{2}  (2.7)
w_{i(u),k} ≡ \frac{w_{i+1,k(w)} + w_{i,k(w)} + w_{i+1,k−1(w)} + w_{i,k−1(w)}}{4}  (2.8)
w_{i(u),k(w)} ≡ \frac{w_{i+1,k(w)} + w_{i,k(w)}}{2}  (2.9)

2.2 空間微分の離散化

2.2.1 2 次精度中心差分

空間微分を 2 次精度差分で離散化する際に必要となる微分操作を以下に示す. ここでは x 方向のフラックス格子点の変数を u_{i(u),k}, z 方向のフラックス格子点の変数を w_{i,k(w)}, スカラー格子点の変数を Π_{i,k} としている. x, z 方向ともにフラックス格子点の変数を ϕ_{i(u),k(w)} としている.

それぞれの変数に対して微分を評価する格子点は一意に決まる. そのため, 他の格子点において微分を評価する場合には平均操作を用いる.

[\frac{∂Π}{∂x}]_{i(u),k} ≡ \frac{Π_{i+1,k} − Π_{i,k}}{Δx}  (2.10)
[\frac{∂Π}{∂z}]_{i,k(w)} ≡ \frac{Π_{i,k+1} − Π_{i,k}}{Δz}  (2.11)
[\frac{∂u}{∂x}]_{i,k} ≡ \frac{u_{i(u),k} − u_{i−1(u),k}}{Δx}  (2.12)
[\frac{∂u}{∂z}]_{i(u),k(w)} ≡ \frac{u_{i(u),k+1} − u_{i(u),k}}{Δz}  (2.13)
[\frac{∂w}{∂x} ]_{i(u),k(w)} ≡ \frac{w_{i+1,k(w)} − w_{i,k(w)}}{Δx}  (2.14)
[\frac{∂w}{∂z} ]_{i,k} ≡ \frac{w_{i,k(w)} − w_{i,k−1(w)}}{Δz}  (2.15)
[\frac{∂ϕ}{∂x}]_{i,k(w)} ≡ \frac{ϕ_{i(u),k(w)} − ϕ_{i−1(u),k(w)}}{Δx}  (2.16)
[\frac{∂ϕ}{∂z}]_{i(u),k} ≡ \frac{ϕ_{i(u),k(w)} − ϕ_{i(u),k−1(w)}}{Δz}  (2.17)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化第 2 章空間方向の離散化 7}

2.2.2 4 次精度中心差分

2 次精度中心差分の場合と同様に, 空間微分を 4 次精度差分で離散化する際に必要となる微分操作を以下に示す.

[\frac{∂Π}{∂x}]_{i(u),k} ≡ \frac{9}{8} (\frac{Π_{i+1,k} − Π_{i,k}}{Δx} ) − \frac{1}{24} (\frac{Π_{i+2,k} − Π_{i−1,k}}{Δx} )  (2.18)
[\frac{∂Π}{∂z}]_{i,k(w)} ≡ \frac{9}{8} (\frac{Π_{i,k+1} − Π_{i,k}}{Δx} ) − \frac{1}{24} (\frac{Π_{i,k+2} − Π_{i,k−1}}{Δx} )  (2.19)
[\frac{∂u}{∂x}]_{i,k} ≡ \frac{9}{8} (\frac{u_{i(u),k} − u_{i−1(u),k}}{Δx} ) − \frac{1}{24} (\frac{u_{i(u)+1,k} − u_{i−2(u),k}}{Δx} )  (2.20)
[\frac{∂u}{∂z}]_{i(u),k(w)} ≡ \frac{9}{8} (\frac{u_{i(u),k+1} − u_{i(u),k}}{Δx} ) − \frac{1}{24} (\frac{u_{i(u),k+2} − u_{i(u),k−1}}{Δx} )  (2.21)
[\frac{∂w}{∂x} ]_{i(u),k(w)} ≡ \frac{9}{8} (\frac{w_{i+1,k(w)} − w_{i,k(w)}}{Δx} ) − \frac{1}{24} (\frac{w_{i+2,k(w)} − w_{i−1,k(w)}}{Δx} ) (2.22)
[\frac{∂w}{∂z} ]_{i,k} ≡ \frac{9}{8} (\frac{w_{i,k(w)} − w_{i,k−1(w)}}{Δz} ) − \frac{1}{24} (\frac{w_{i,k+1(w)} − w_{i,k−2(w)}}{Δz} )  (2.23)

[\frac{∂ϕ}{∂x}]_{i,k(w)} ≡ \frac{9}{8} (\frac{ϕ_{i(u),k(w)} − ϕ_{i−1(u),k(w)}}{Δx} ) − \frac{1}{24} (\frac{ϕ_{i+1(u),k(w)} − ϕ_{i−2(u),k(w)}}{Δx} )

(2.24)

[\frac{∂ϕ}{∂z}]_{i(u),k} ≡ \frac{9}{8} (\frac{ϕ_{i(u),k(w)} − ϕ_{i(u),k−1(w)}}{Δz} ) − \frac{1}{24} (\frac{ϕ_{i(u),k+1(w)} − ϕ_{i(u),k−2(w)}}{Δz} )

(2.25)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化第 2 章空間方向の離散化 8}

2.3 空間離散化した基礎方程式

2.3.1 静水圧の式

[\frac{∂\overline{Π}}{∂z}]_{i,k} = −\frac{g}{⟨c_{p}⟩[\overline{θ_{v}}]_{i,k}}  (2.26)

基本場の密度\overline{ρ}_{i,k} は以下のように計算する.

\overline{ρ}_{i,k} = \frac{p_{0}}{⟨R⟩} \frac{[\overline{Π}^{⟨c_{v}⟩/⟨R⟩}]_{i,k}}{[\overline{θ_{v}}]_{i,k}}  (2.27)

2.3.2 運動方程式

\frac{∂u_{i(u),k}}{∂t} = −u_{i(u),k} [\frac{∂u}{∂x}]_{i(u),k} − w_{i(u),k} [\frac{∂u}{∂z}]_{i(u),k}
−⟨c_{p}⟩[\overline{θ_{v}}]_{i(u),k} [\frac{∂Π}{∂x}]_{i(u),k} + [D\overline{u}]_{i(u),k} + [D_{u}]_{i(u),k}  (2.28)
\frac{∂w_{i,k(w)}}{∂t} = −u_{i,k(w)} [\frac{∂w}{∂x} ]_{i,k(w)} − w_{i,k(w)} [\frac{∂w}{∂z} ]_{i,k(w)}

−⟨c_{p}⟩[\overline{θ_{v}}]_{i,k(w)} [\frac{∂Π}{∂z}]_{i,k(w)} + [D_{w}]_{i,k(w)}

+g \frac{θ_{i,k(w)}}{\overline{θ}_{i,k(w)}} g

+\frac{[q_{d}]_{i,k(w)}⟨M⟩/M_{d} + ∑[q_{v}]_{i,k(w)}⟨M⟩/M_{v}}{[\overline{q_{d}}]_{i,k(w)}⟨M⟩/M_{d} + ∑[\overline{q_{v}}]_{i,k(w)}⟨M⟩/M_{v}} g  (2.29)

2.3.3 圧力方程式

\frac{∂Π_{i,k}}{∂t} = −\frac{\overline{⟨c_{s}⟩^{2}}_{i,k}}{⟨c_{p}⟩ρ̄_{i,k}[θ̄_{v}^{2}]_{i,k}} [∇· (ρ̄θ̄_{v}u)]_{i,k}

+\frac{\overline{⟨c_{s}⟩^{2}}_{i,k}}{⟨c_{p}⟩[θ̄_{v}]_{I,k}} [\frac{1}{ρ̄_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

+\frac{1}{θ̄_{i,k}} \{\frac{1}{Π̄_{i,k}} ([Q_{cond}]_{i,k} + [Q_{rad}]_{i,k} + [Q_{dis}]_{i,k}) + [D_{θ}]_{i,k}\}

+\frac{1}{[q̄_{d}]_{i,k}⟨M⟩/M_{d} + ∑[q̄_{v}]_{i,k}⟨M⟩/M_{v}}

\{\frac{⟨M⟩}{M_{d}} (−\frac{[q̄_{d}]_{i,k}}{ρ̄_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k} + [D_{q_{d}}]_{i,k})

+ ∑ \frac{⟨M⟩}{M_{v}} (−\frac{[q̄_{v}]_{i,k}}{ρ̄_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k} + [M_{src}(q_{v})]_{i,k} + [D_{q_{v}}]_{i,k}) \}].

(2.30)

基本場の音速\overline{⟨c_{s}⟩} は以下のように計算する.

\overline{⟨c_{s}⟩}_{i,k}^{2} = \frac{⟨c_{p}⟩⟨R⟩}{⟨c_{v}⟩}\overline{Π}_{i,k}[\overline{θ_{v}}]_{i,k}.  (2.31)

2.3.4 熱力学の式

\frac{∂θ_{i,k}}{∂t} = −u_{i,k} [\frac{∂θ}{∂x} ]_{i,k} − w_{i,k} [\frac{∂θ}{∂z} ]_{i,k} − w_{i,k} [\frac{∂\overline{θ}}{∂z} ]_{i,k}
+\frac{1}{\overline{Π}_{i,k}} ([Q_{cond}]_{i,k} + [Q_{rad}]_{i,k} + [Q_{dis}]_{i,k}) + [D_{θ}]_{i,k} + [D\overline{θ}]_{i,k} (2.32)

2.3.5 比湿の式

\frac{∂[q_{d}]_{i,k}}{∂t} = −u_{i,k} [\frac{∂q_{d}}{∂x} ]_{i,k} − w_{i,k} [\frac{∂q_{d}}{∂z} ]_{i,k} − w_{i,k} [\frac{∂\overline{q_{d}}}{∂z} ]_{i,k}
−\frac{[\overline{q_{d}}]_{i,k}}{\overline{ρ}_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k} + [D_{q_{d}}]_{i,k},  (2.33)
\frac{∂[q_{v}]_{i,k}}{∂t} = −u_{i,k} [\frac{∂q_{v}}{∂x} ]_{i,k} − w_{i,k} [\frac{∂q_{v}}{∂z} ]_{i,k} − \overline{u}_{i,k} [\frac{∂q_{v}}{∂x} ]_{i,k} − w_{i,k} [\frac{∂\overline{q_{v}}}{∂z} ]_{i,k}

−\frac{[\overline{q_{v}}]_{i,k}}{\overline{ρ}_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k} + \frac{1}{\overline{ρ}_{i,k}} [M_{src}(ρ_{v})]_{i,k} + [D\overline{q_{v}}]_{i,k} + [D_{q_{v}}]_{i,k},

(2.34)

\frac{∂[q_{c}]_{i,k}}{∂t} = −u_{i,k} [\frac{∂q_{c}}{∂x} ]_{i,k} − w_{i,k} [\frac{∂q_{c}}{∂z} ]_{i,k} − \frac{[\overline{q_{c}}]_{i,k}}{\overline{ρ}_{i,k}} ∑[M_{fall}(ρ_{r})]_{i,k}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

+\frac{1}{\overline{ρ}_{i,k}} [M_{src}(ρ_{c})]_{i,k} + [D_{q_{c}}]_{i,k},  (2.35)
\frac{∂[q_{r}]_{i,k}}{∂t} = −u_{i,k} [\frac{∂q_{r}}{∂x} ]_{i,k} − w_{i,k} [\frac{∂q_{r}}{∂z} ]_{i,k}
+\frac{1}{\overline{ρ}_{i,k}} [M_{fall}(ρ_{r})]_{i,k} + \frac{1}{\overline{ρ}_{i,k}} [M_{src}(ρ_{r})]_{i,k} + [D_{q_{r}}]_{i,k}.  (2.36)

2.4 境界条件

ここでは離散化した変数に対する境界条件の与え方をまとめる. 考慮する境界条件は, 周期境界条件, 境界ですべりなし条件と応力なし条件である.

2.4.1 周期境界条件の与え方

例として, x 方向フラックス格子点に配置された変数 u_{i(u),k} を考える. 計算領域内の x 方向の添字を 1(u) ∼ im(u) とし, 糊代部分の格子点数を 2 とする (図 1.2 参照). このとき周期境界条件は以下のように与えられる.

u_{0(u),k} = u_{im(u),k}  (2.37)
u_{−1(u),k} = u_{im−1(u),k}  (2.38)
u_{im+1(u),k} = u_{1(u),k}  (2.39)
u_{im+2(u),k} = u_{2(u),k}  (2.40)

ただし k は任意の整数であり, その範囲は −1 ≤ k ≤ km + 2 である.

z 方向フラックス格子点に配置された変数, スカラー格子点に配置された変数に対しても同様に与えることができる.

2.4.2 すべりなし条件の与え方

境界で速度を 0 とする. この場合, 境界をはさんで変数の値が反対称になるように与える.

例として x 方向に境界を与えた場合を考える. x 方向フラックス格子点に配置された変数に対しては.

u_{0(u),k} = u_{im(u),k} = 0  (2.41)
u_{−1(u),k} = −u_{1(u),k}  (2.42)
u_{im+1(u),k} = −u_{im−1(u),k}  (2.43)
u_{im+2(u),k} = −u_{im−2(u),k}  (2.44)

とする. 境界上に配置されていない変数に対しては,

Π_{0,k} = −Π_{1,k}  (2.45)

Π_{−1,k} = −Π_{1,k}  (2.46)
Π_{im+1,k} = −Π_{im,k}  (2.47)
Π_{im+2,k} = −Π_{im−1,k}  (2.48)

とする.

2.4.3 応力なし条件の与え方

境界上で法線方向速度を 0, 接線方向速度の法線方向微分を 0 とする. この場合, 境界上で配置された速度成分は境界をはさんで変数の値が反対称になるように与え, 境界上に配置されていない変数に対しては壁をはさんで変数の値が対称になるように与える.

例として x 方向に境界を与えた場合を考える. x 方向フラックス格子点に配置された変数に対しては.

u_{0(u),k} = u_{im(u),k} = 0  (2.49)
u_{−1(u),k} = −u_{1(u),k}  (2.50)
u_{im+1(u),k} = −u_{im−1(u),k}  (2.51)
u_{im+2(u),k} = −u_{im−2(u),k}  (2.52)

とする. 境界上に配置されていない変数に対しては,

Π_{0,k} = Π_{1,k}  (2.53)
Π_{−1,k} = Π_{1,k}  (2.54)
Π_{im+1,k} = Π_{im,k}  (2.55)
Π_{im+2,k} = Π_{im−1,k}  (2.56)

とする.

\underline{2 次元非静力学モデルの離散化第 3 章時間方向の離散化 13}

第3章時間方向の離散化

3.1 運動方程式と圧力方程式

空間離散化された運動方程式 (2.28), (2.29) と圧力方程式 (2.30) を時間方向に離散化する. 音波に関連する項は短いタイムステップ Δτ で離散化し, その他の項は長いタイムステップ Δt で離散化する. 音波に関連する項の離散化には HE-VI 法を採用し, u の式は前進差分, w, Π の式は後退差分 (クランク・ニコルソン法) で離散化する. その他の項の離散化にはリープフロッグ法を用いる. 離散化した式の計算はまず u の式から行う. 得られた τ + Δτ の u を用いて Π を計算し, u, Π を用いて w を計算する.

運動方程式の各項のうち, 音波に関係しない項を F_{u}, F_{w} としてまとめると, 運動方程式と圧力方程式は以下のように書ける.

\frac{∂u_{i(u),k}}{∂t} = − [⟨c_{p}⟩θ̄_{v} \frac{∂(Π − αDiv)}{∂x}]_{i(u),k} + [F_{u}]_{i(u),k}^{t},  (3.1)
\frac{∂w_{i,k(w)}}{∂t} = − [⟨c_{p}⟩θ̄_{v} \frac{∂(Π − αDiv)}{∂z}]_{i,k(w)} + [F_{w}]_{i,k(w)}^{t},  (3.2)
\frac{∂Π_{i,k}}{∂t} + [\frac{\overline{⟨c_{s}⟩}^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} \frac{∂(ρ̄θ̄_{v}w)}{∂z}]_{i,k} = − [\frac{\overline{⟨c_{s}⟩}^{2}}{⟨c_{p}⟩θ̄_{v}} \frac{∂u}{∂x}]_{i,k} .  (3.3)

ただし u, w の式には音波減衰項 αDiv を加えてある(Skamarock and Klemp, 1992). 音波に関連しない項 F_{u}, F_{w} は,

[F_{u}]_{i(u),k}^{t} = − [Adv.u]_{i(u),k}^{t} + [Turb.u]_{i(u),k}^{t−Δt} + [Diff.u]_{i(u),k}^{t−Δt}  (3.4)

[F_{w}]_{i,k(w)}^{t} = − [Adv.w]_{i,k(w)}^{t} + [Buoy]_{i,k(w)}^{t} + [Turb.w]_{i,k(w)}^{t−Δt} . + [Diff.w]_{i,k(w)}^{t−Δt} . (3.5)

であり, 時刻 t で評価することにする. 但し, 中心差分でリープフロッグ法を用いるため, 数値粘性項 Diff を追加してある.

\underline{2 次元非静力学モデルの離散化第 3 章時間方向の離散化 14}

3.1.1 音波に関連する項の時間方向の離散化

水平方向の運動方程式の離散化 (3.1) を時間方向に離散化すると以下のようになる.

u_{i(u),k}^{τ+Δτ} = u_{i(u),k}^{τ} − [⟨c_{p}⟩θ̄_{v}Δτ \{\frac{∂Π^{τ}}{∂x} − \frac{∂(αDiv)^{τ}}{∂x} \}]_{i(u),k} + F_{u,i(u),k}^{t}Δτ (3.6)

鉛直方向の運動方程式と圧力方程式の離散化 HE-VI 法を用いるので, w と Π の式を連立して解く. w の式において音波減衰項は前進差分, 圧力項は後退差分で離散化する. Π の式において水平微分項は (3.6) で求めた u^{τ+Δτ} を用いて離散化し, 鉛直微分項は後退差分で離散化する.

w_{i,k(w)}^{τ+Δτ} = w_{i,k(w)}^{τ} − ⟨c_{p}⟩θ̄_{v}Δτ \{β \frac{∂Π^{τ+Δτ}}{∂z} + (1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \}
+F_{w,i,k(w)}^{t}Δτ.  (3.7)
Π^{τ+Δτ} + β \frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} \frac{∂(ρ̄θ̄_{v}w^{τ+Δτ})}{∂z} = Π^{τ} − (1 − β)\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} \frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}

−\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} \frac{∂u^{τ+Δτ}}{∂x} − \overline{u}Δτ \frac{∂Π^{τ}}{∂x} + H_{i,k}^{τ}Δτ.

(3.8)

但し

H_{i,k}^{τ} = \frac{c̄_{s}^{2}}{⟨c_{p}⟩θ̄_{v}} [\frac{1}{ρ̄} ∑ M_{fall}^{τ} + \frac{1}{θ̄}\{\frac{1}{Π̄}(Q_{cond}^{τ} + Q_{rad}^{τ} + Q_{dis}^{τ}) + D_{θ}^{t−Δt}\}

\frac{1}{q̄_{d}⟨M⟩/M_{d} + ∑ q̄_{v}⟨M⟩/M_{v}} \{\frac{⟨M⟩}{M_{d}}( − \frac{q̄_{d}}{ρ̄} ∑ M_{fall}^{τ} + D_{q_{d}}^{t−Δt})

+ ∑ \frac{⟨M⟩}{M_{v}}( − \frac{q̄_{v}}{ρ̄} ∑ M_{fall}^{τ} + M_{src}^{τ}(q_{v})D_{q_{v}}^{t−Δt})\}]_{i,k}  (3.9)

である. (3.8) 式に (3.7) を代入して w^{τ+Δτ} を消去する.

Π^{τ+Δτ} − β^{2} \frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} \frac{∂}{∂z} \{(⟨c_{p}⟩ρ̄θ̄_{v}^{2}) (\frac{∂Π^{τ+Δτ}}{∂z} )\}

= Π^{τ} − (1 − β)\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} \frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z} − \frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} \frac{∂u^{τ+Δτ}}{∂x}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

−β \frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} [\frac{∂}{∂z} ρ̄θ̄_{v} \{w^{τ} − ⟨c_{p}⟩θ̄_{v}Δτ \{(1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \} + F_{w}^{t}Δτ \}]

−\overline{u}Δτ \frac{∂Π^{τ}}{∂x} + H_{i,k}^{τ}Δτ.  (3.10)

(3.10) 式右辺を空間方向に離散化し, 格子点位置を表す添字を付けて表すと以下のようになる (計算の詳細は第 A 章参照).

\{−β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k(w)}\} Π_{i,k+1}^{τ+Δτ}

+ [1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{2}} \{(⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k(w)} + (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k−1(w)}\}] Π_{i,k}^{τ+Δτ}

+ \{−β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k−1(w)}\} Π_{i,k−1}^{τ+Δτ}

= Π_{i,k}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,k} − (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{k} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,k}

−β (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} [\frac{∂}{∂z} (ρ̄θ̄_{v})_{i,k(w)} \{w_{i,k(w)}^{τ}

− (⟨c_{p}⟩θ̄_{v})_{i,k(w)} Δτ \{(1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \}_{i,k(w)} + (F_{w}^{t})_{i,k(w)} Δτ\}]_{i,k}

− (\overline{u}Δτ \frac{∂Π^{τ}}{∂x} )_{i,k} + H_{i,k}^{τ}Δτ.

(3.11)

但し平均場の量は鉛直方向にしか依存しないので z 方向の添字のみ付けてある.

境界条件

上下境界を固定壁とする場合, 境界条件は上部下部境界で,

w(i, 0(w)) = 0,  (3.12)
w(i, km(w)) = 0  (3.13)

である.

下部境界:

下部境界 (k(w) = 0(w)) について考える. この時 (3.7) 式に添字を付けて書き下 \underline{2 次元非静力学モデルの離散化第 3 章時間方向の離散化 16} すと,

β (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,0(w)} = (\frac{∂(αDiv)^{τ}}{∂z} )_{i,0(w)} − (1 − β) (\frac{∂Π^{τ}}{∂z} )_{i,0(w)} + (\frac{1}{⟨c_{p}⟩θ̄_{v}} F_{w}^{t})_{i,0(w)}
≡ E_{i,0(w)}  (3.14)

となる. したがって (3.11) 式は以下のようになる.

\{−β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{1(w)}\} Π_{i,2}^{τ+Δτ}

+ \{1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{1(w)}\} Π_{i,1}^{τ+Δτ}

= Π_{i,1}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,1} − (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{1} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,1}

−β (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1}

[\frac{∂}{∂z} ρ̄θ̄_{v} \{w^{τ} − ⟨c_{p}⟩θ̄_{v}Δτ \{(1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \} + F_{w}^{t}Δτ \}]_{i,1}

−β (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \frac{1}{Δz} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{i,0(w)} E_{i,0(w)} − (\overline{u}\frac{∂Π^{τ}}{∂x} )_{i,1} Δτ + H_{i,1}^{τ}Δτ.

(3.15)

上部境界:

上部境界 (k(w) = km(w)) について考える. この時 (3.7) 式を添字を付けて書き下すと,

β (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,km(w)} = (\frac{∂(αDiv)^{τ}}{∂z} )_{i,km(w)} − (1 − β) (\frac{∂Π^{τ}}{∂z} )_{i,km(w)}

+ (\frac{1}{⟨c_{p}⟩θ̄_{v}} F_{w}^{t})_{i,km(w)}

≡ E_{i,km(w)}  (3.16)

となる. したがって (3.11) 式は以下のようになる.

\{1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{km−1(w)}\} Π_{i,km}^{τ+Δτ}

+ \{−β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{km−1(w)}\} Π_{i,km−1}^{τ+Δτ}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

= Π_{i,km}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,km} − (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{km} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,km}

−β (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km}

[\frac{∂}{∂z} ρ̄θ̄_{v} \{w^{τ} − ⟨c_{p}⟩θ̄_{v}Δτ \{(1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \} + F_{w}^{t}Δτ \}]_{i,km}

+\frac{β}{Δz} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{km(w)} E_{i,km(w)} − (\overline{u}\frac{∂Π^{τ}}{∂x} )_{i,km} Δτ + H_{i,km}^{τ}Δτ.

(3.17)

圧力方程式の時間積分方法

(3.11), (3.15), (3.17) 式を連立すると, 以下のような行列式の形式で書くことができる.

( \begin{matrix}A_{1} & B_{2} & & 0 \\ C_{1} & ... & ... \\ & ... & ... & B_{km} \\ 0 & & C_{km−1} & A_{km} \end{matrix} ) ( \begin{matrix}Π_{1,1} & Π_{2,1} & ⋯ & Π_{im,1} \\ Π_{1,2} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ Π_{1,km} & ⋯ & ⋯ & Π_{im,km} \end{matrix} )^{τ+Δτ}

= ( \begin{matrix}D_{1,1} & D_{2,1} & ⋯ & D_{im,1} \\ D_{1,2} & ... & ... & ⋮ \\ ⋮ & & ... & ⋮ \\ D_{1,km} & ⋯ & ⋯ & D_{im,km} \end{matrix} )^{τ} .  (3.18)

この連立方程式を解くことで Π_{i,k} を求める. この連立方程式の係数は以下のように書ける.

A_{k} = 1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \frac{1}{Δz^{2}} \{(⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k(w)} + (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k−1(w)}\}
(k = 2, 3, ⋯ km − 1),
A_{1} = 1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{1(w)} ,
A_{km} = 1 + β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{km−1(w)} ,
B_{k} = −β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k−1} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k−1(w)} ,
(k = 2, 3, ⋯ km),

C_{k} = −β^{2} (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k+1} \frac{1}{Δz^{2}} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{k(w)} ,

(k = 1, 2, ⋯ km − 1),

D_{i,k} = Π_{i,k}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,k}

− (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{k} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,k} + F_{i,k} (k = 2, 3, ⋯ km − 1),

D_{i,1} = Π_{i,1}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,1}

− (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{1} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,1} + F_{i,1}

−β (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{1} \frac{1}{Δz} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{i,0(w)} E_{i,0(w)},

D_{i,km} = Π_{i,km}^{τ} − (1 − β) (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \{\frac{∂(ρ̄θ̄_{v}w^{τ})}{∂z}\}_{i,km}

− (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩θ̄_{v}} )_{km} (\frac{∂u^{τ+Δτ}}{∂x} )_{i,km} + F_{i,km}

+β (\frac{c̄_{s}^{2}Δτ^{2}}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{km} \frac{1}{Δz} (⟨c_{p}⟩ρ̄θ̄_{v}^{2})_{km(w)} E_{i,km(w)}.

ただし,

E_{i,k(w)} ≡ (\frac{∂(αDiv)^{τ}}{∂z} )_{i,k(w)} − (1 − β) (\frac{∂Π^{τ}}{∂z} )_{i,k(w)} + (\frac{1}{⟨c_{p}⟩θ̄_{v}} F_{w}^{t})_{i,k(w)}

F_{i,k} ≡ −β (\frac{c̄_{s}^{2}Δτ}{⟨c_{p}⟩ρ̄θ̄_{v}^{2}} )_{k}

[\frac{∂}{∂z} (ρ̄θ̄_{v})_{i,k(w)} \{w_{i,k(w)}^{τ}

− (⟨c_{p}⟩θ̄_{v})_{i,k(w)} Δτ \{(1 − β)\frac{∂Π^{τ}}{∂z} − \frac{∂(αDiv)^{τ}}{∂z} \}_{i,k(w)} + (F_{w}^{t})_{i,k(w)} Δτ \}]_{i,k}.

− (\overline{u}\frac{∂Π^{τ}}{∂x} )_{i,k} Δτ + H_{i,k}^{τ}

H_{i,k}^{τ} ≡ \frac{c̄_{s}^{2}}{⟨c_{p}⟩θ̄_{v}} [\frac{1}{ρ̄} ∑ M_{fall}^{τ} + \frac{1}{θ̄}\{\frac{1}{Π̄}(Q_{cond}^{τ} + Q_{rad}^{τ} + Q_{dis}^{τ}) + D_{θ}^{t−Δt}\}

\frac{1}{q̄_{d}⟨M⟩/M_{d} + ∑ q̄_{v}⟨M⟩/M_{v}} \{\frac{⟨M⟩}{M_{d}}( − \frac{q̄_{d}}{ρ̄} ∑ M_{fall}^{τ} + D_{q_{d}}^{t−Δt})

+ ∑ \frac{⟨M⟩}{M_{v}}( − \frac{q̄_{v}}{ρ̄} ∑ M_{fall}^{τ} + M_{src}^{τ}(q_{v})D_{q_{v}}^{t−Δt})\}]_{i,k}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化第 3 章時間方向の離散化 19} である.

3.1.2 音波に関連しない項の時間方向の離散化

運動方程式の音波に関連しない項 (3.1), (3.2) 式を離散化する.

F_{u,i(u),k}^{t} = − [Adv.u]_{i(u),k}^{t} + [Turb.u]_{i(u),k}^{t−Δt} + [Diff.u]_{i(u),k}^{t−Δt} ,  (3.19)
F_{w,i,k(w)}^{t} = − [Adv.w]_{i(u),k}^{t} + [Buoy]_{i,k(w)}^{t} + [Turb.w]_{i,k(w)}^{t−Δt} + [Diff.w]_{i,k(w)}^{t−Δt} .

(3.20)

ここで, Adv は移流項, Turb は粘性拡散項, Buoy は浮力項, Diff は数値粘性項である. それぞれの項を書き下すと,

[Adv.u]_{i(u),k}^{t} = u_{i(u),k}^{t} [\frac{∂u}{∂x}]_{i(u),k}^{t} + w_{i(u),k}^{t} [\frac{∂u}{∂z}]_{i(u),k}^{t}
+\overline{u}_{i(u),k} [\frac{∂u}{∂x}]_{i(u),k}^{t} + w_{i(u),k}^{t} [\frac{∂\overline{u}}{∂z}]_{i(u),k} ,  (3.21)
[Adv.w]_{i,k(w)}^{t} = u_{i,k(w)}^{t} [\frac{∂w}{∂x} ]_{i,k(w)}^{t} + w_{i,k(w)}^{t} [\frac{∂w}{∂z} ]_{i,k(w)}^{t}
+\overline{u}_{i,k(w)} [\frac{∂w}{∂x} ]_{i,k(w)}^{t}  (3.22)

であり, 浮力項は,

[Buoy]_{i,k(w)}^{t} = g \frac{θ_{i,k(w)}^{t}}{\overline{θ}_{i,k(w)}}
+g \frac{[q_{d}]_{i,k(w)}^{t}⟨M⟩/M_{d} + ∑[q_{v}]_{i,k(w)}^{t}⟨M⟩/M_{v}}{[\overline{q_{d}}]_{i,k(w)}⟨M⟩/M_{d} + ∑[\overline{q_{v}}]_{i,k(w)}^{t}⟨M⟩/M_{v}}  (3.23)

であり, 粘性拡散項は,

[Turb.u]_{i(u),k}^{t−Δt}
= 2 [\frac{∂}{∂x} \{(K_{m})_{i,k} (\frac{∂u}{∂x})_{i,k}\}]_{i(u),k}^{t−Δt}

+ [\frac{∂}{∂z} \{(K_{m})_{i(u),k(w)} (\frac{∂w}{∂x} )_{i(u),k(w)} + (K_{m})_{i(u),k(w)} (\frac{∂u}{∂z})_{i(u),k(w)}\}]_{i(u),k}^{t−Δt}

−\frac{2}{3C_{m}^{2}l^{2}} (\frac{∂K_{m}^{2}}{∂x} )_{i(u),k}^{t−Δt}  (3.24)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

[Turb.w]_{i,k(w)}^{t−Δt}
= 2 [\frac{∂}{∂z} \{(K_{m})_{i,k} (\frac{∂w}{∂z} )_{i,k}\}]_{i,k(w)}^{t−Δt}

+ [\frac{∂}{∂x} \{(K_{m})_{i(u),k(w)} (\frac{∂w}{∂x} )_{i(u),k(w)} + (K_{m})_{i(u),k(w)} (\frac{∂u}{∂z})_{i(u),k(w)}\}]_{i,k(w)}^{t−Δt}

−\frac{2}{3C_{m}^{2}l^{2}} (\frac{∂K_{m}^{2}}{∂z} )_{i,k(w)}^{t−Δt}  (3.25)

である. 数値粘性項は,

[Diff.u]_{i(u),k}^{t−Δt} = ν_{h} \{\frac{∂}{∂x} (\frac{∂u}{∂x})_{i,k}\}_{i(u),k}^{t−Δt} + ν_{v} \{\frac{∂}{∂z} (\frac{∂u}{∂z})_{i(u),k(w)}\}_{i(u),k}^{t−Δt} (3.26)

[Diff.w]_{i,k(w)}^{t−Δt} = ν_{h} \{\frac{∂}{∂x} (\frac{∂w}{∂x} )_{i(u),k(w)}\}_{i,k(w)}^{t−Δt} + ν_{v} \{\frac{∂}{∂z} (\frac{∂w}{∂z} )_{i,k}\}_{i,k(w)}^{t−Δt}(3.27)

である. K_{m} は乱流エネルギーの時間発展方程式から計算し (詳細は後述), ν_{h}, ν_{v} は以下のように定める.

ν_{h} = \frac{α_{h}Δx^{2}}{Δt}  (3.28)
ν_{v} = \frac{α_{v}Δz^{2}}{Δt}  (3.29)

ここで Δx, Δz は水平・鉛直方向の格子間隔を意味し, α_{h}, α_{v} はそれぞれ,

α_{h} ≤ \frac{1}{8},  α_{v} ≤ \frac{1}{8}  (3.30)

とする^{1}.

3.2 熱力学の式と混合比の保存式の離散化

熱の式と混合比の保存式の右辺をまとめて F で表し, 時間方向にリープフロッグ法を用いて離散化する.

θ_{i,k}^{t+Δt} = θ_{i,k}^{t−Δt} + 2Δt[F_{θ}]_{i,k}^{t}  (3.31)
[q_{v}]_{i,k}^{t+Δt} = [q_{v}]_{i,k}^{t−Δt} + 2Δt[F_{q_{v}}]_{i,k}^{t}  (3.32)
[q_{c}]_{i,k}^{t+Δt} = [q_{c}]_{i,k}^{t−Δt} + 2Δt[F_{q_{c}}]_{i,k}^{t}  (3.33)
[q_{r}]_{i,k}^{t+Δt} = [q_{r}]_{i,k}^{t−Δt} + 2Δt[F_{q_{r}}]_{i,k}^{t}  (3.34)

1(3.30) は CreSS のドキュメントにも現れるが, その導出方法は良く分からない.

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化第 3 章時間方向の離散化 21} ここで,

[F_{θ}]_{i,k} = − [Adv.θ]_{i,k}^{t} − [Adv.θ̄]_{i,k}^{t} + [Turb.θ]_{i,k}^{t−Δt} + [Turb.θ̄]_{i,k}^{t−Δt} + [Diff.θ]_{i,k}^{t−Δt}
+[Q_{cnd}]_{i,k}^{t} + [Q_{rad}]_{i,k}^{t−Δt} + [Q_{dis}]_{i,k}^{t−Δt}  (3.35)
[F_{q_{v}}]_{i,k}^{t} = − [Adv.q_{v}]_{i,k}^{t} − [Adv.q̄_{v}]_{i,k}^{t} + [Turb.q_{v}]_{i,k}^{t−Δt} + [Turb.q̄_{v}]_{i,k}^{t−Δt}
+ [Diff.q_{v}]_{i,k}^{t−Δt} + [EV_{rv}]_{i,k}^{t}  (3.36)
[F_{q_{c}}]_{i,k}^{t} = − [Adv.q_{c}]_{i,k}^{t} + [Turb.q_{c}]^{t−Δt} + [Diff.q_{c}]^{t−Δt}
− [CN_{cr} + CL_{cr}]_{i,k}^{t}  (3.37)
[F_{q_{r}}]_{i,k}^{t} = − [Adv.q_{r}]_{i,k}^{t} + [Turb.q_{r}]_{i,k}^{t−Δt} + [Diff.q_{c}]^{t−Δt}
+ [CN_{cr} + CL_{cr} − EV_{rv}]_{i,k}^{t} + [P R_{r}]_{i,k}^{t}  (3.38)

である. 移流を中心差分で安定して解くために, 数値粘性項 Diff を追加してある. また, CN_{vc}, EV_{cv} 項は湿潤飽和調節法より決めるため, それらの項を含めない. θ, q_{v}, q_{c}, q_{r} をまとめて ϕ で表し, それぞれの項を書き下す. 移流項は,

[Adv.ϕ]_{i,k}^{t} = [u_{i(u),k} [\frac{∂ϕ}{∂x}]_{i(u),k}]_{i,k}^{t} + [w_{i,k(w)} [\frac{∂ϕ}{∂z}]_{i,k(w)}]_{i,k}^{t}  (3.39)

であり, 基本場の移流項は,

[Adv.ϕ̄]_{i,k}^{t} = [w_{i,k(w)} [\frac{∂\overline{ϕ}}{∂z}]_{i,k(w)}]_{i,k}^{t}  (3.40)

である. 粘性拡散項は CReSS と同様に 1.5 次のクロージャーを用いることで,

[Turb.ϕ]_{i,k}^{t−Δt} = [\frac{∂}{∂x} \{(K_{h})_{i(u),k} (\frac{∂ϕ}{∂x})_{i(u),k}\}]_{i,k}^{t−Δt}
+ [\frac{∂}{∂z} \{(K_{h})_{i,k(w)} (\frac{∂ϕ}{∂z})_{i,k(w)}\}]_{i,k}^{t−Δt}  (3.41)

となり, 基本場の粘性拡散項は,

[Turb.ϕ̄]_{i,k}^{t−Δt} = [\frac{∂}{∂z} \{(K_{h})_{i,k(w)} (\frac{∂\overline{ϕ}}{∂z})_{i,k(w)}\}]_{i,k}^{t−Δt}  (3.42)

となる. 数値粘性項は,

[Diff_{ϕ}]_{i,k}^{t−Δt} = ν_{h} \{\frac{∂}{∂x} (\frac{∂ϕ}{∂x})_{i(u),k}\}_{i,k}^{t−Δt} + ν_{v} \{\frac{∂}{∂z} (\frac{∂ϕ}{∂z})_{i,k(w)}\}_{i,k}^{t−Δt} (3.43)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

である. K_{h} は乱流エネルギーの時間発展方程式から計算する (詳細は後述). ν_{h}, ν_{v} は (3.30) 式を利用する.

凝縮加熱項 Q_{cnd} は

[Q_{cnd}]_{i,k}^{t} = − [\frac{L}{c_{p}Π̄_{d}Π̄}EV_{rv}]_{i,k}^{t} = −\frac{L}{c_{p}_{d}Π̄} \{4.85 × 10^{−2}([q_{vsw}]_{i,k}^{t} − [q_{v}]_{i,k}^{t})(ρ̄_{i,k}[q_{r}]_{i,k})^{0.65}\}

(3.44)

である.

散逸加熱項 Q_{dis} は

[Q_{dis}]_{i,k}^{t−Δt} = \frac{1}{c_{p}_{d}Π̄} \frac{C_{ε}}{l} \frac{(K_{m,i,k}^{t−Δt})^{3}}{(C_{m}l)^{3}} = \frac{1}{c_{p}_{d}Π̄} \frac{(K_{m,i,k}^{t−Δt})^{3}}{C_{m}^{2}l^{4}}  (3.45)

と与える. ここで l = (ΔxΔz)^{1/2} である.

放射強制 [Q_{rad}]_{i,k} は計算設定ごとに与える.

雲水から雨水への変換を表す CN_{cr}, CL_{cr} は以下のようになる.

[CN_{cr}]_{i,k}^{t} = (q_{c,i,k}^{t} − q_{c0})/τ_{ac}  (3.46)
[CL_{cr}]_{i,k}^{t} = 2.2[q_{c}]_{i,k}^{t} ([ρ̄]_{i,k} [q_{r}]_{i,k}^{t})^{0.875}  (3.47)

雨水の蒸発を表す EV_{rv} は以下のようになる.

[EV_{rv}]_{i,k}^{t} = 4.85 × 10^{−2}([q_{vsw}]_{i,k}^{t} − [q_{v}]_{i,k}^{t})([ρ̄]_{i,k}[q_{r}]_{i,k}^{t})^{0.65}  (3.48)

降水による雨水フラックスを表す P R_{r} は以下のように書ける.

[P R_{r}]_{i,k}^{t} = \frac{1}{[ρ̄]_{i,k}} \frac{∂}{∂z} ([ρ̄]_{i,k}[U_{r}]_{i,k}^{t}[q_{r}]_{i,k}^{t}).  (3.49)
[U_{r}]_{i,k}^{t} = 12.2([q_{r}]_{i,k}^{t})^{0.125}  (3.50)

3.2.1 湿潤飽和調節法

Klemp and Wilhelmson (1978), CReSS ユーザーマニュアル (坪木と榊原, 2001) では, 水蒸気と雲水の間の変換を表す −CN_{vc} + EV_{cv} は, Soong and Ogura (1973) において開発された湿潤飽和調節法を用いている. この方法は dS = 0 の断熱線と, μ_{vapar} = μ_{condensedphase} の平衡条件 (μ は化学ポテンシャル) の交わる温度・圧力・組成を反復的に求める数値解法である. 以下ではそのやり方を解説する.

飽和蒸気圧を用いる場合

湿潤飽和調節法を用いる場合, まず始めに (3.31) – (3.34) 式から求まる量に ∗ を添付し, [θ]^{∗}, [q_{v}]^{∗}, [q_{c}]^{∗}, [q_{r}]^{∗} とする. 水に対する過飽和混合比

Δq_{c} = MAX{0, [q_{v}]^{∗} − q_{vsw}([θ]^{∗})}  (3.51)

が Δq_{c} > 0, もしくは雲粒混合比が q_{c}^{∗} > 0 ならば, 次式を用いて暫定的に θ, q_{v}, q_{c} を求める.

[θ]^{t+Δt} = θ^{∗} + \frac{γ([q_{v}]^{∗} − q_{vsw}([θ]^{∗}))}{1 + γ \frac{∂q_{vsw}([θ]^{∗})}{∂θ}}  (3.52)
[q_{v}]^{t+Δt} = [q_{v}]^{∗} + \frac{[θ]^{∗} − [θ]^{t+Δt}}{γ} ,  (3.53)
[q_{c}]^{t+Δt} = [q_{v}]^{∗} + [q_{c}]^{∗} − [q_{v}]^{t+Δt}.  (3.54)

ただし, γ = L_{v}/(c_{p}Π) である. もしも [q_{c}]^{t+Δt} > 0 ならば, 暫定的に得られた値を ∗ 付きのものに置き換え, (3.52) – (3.54) 式の値が収束するまで繰り返し適用する. Soong and Ogura(1973) によると, 通常高々数回繰り返せば収束し, 調整後の値が得られるそうである.

もしも q_{c}^{t+Δt} < 0 の場合には,

[θ]^{t+Δt} = [θ]^{∗} − γ[q_{c}]^{∗},  (3.55)
[q_{v}]^{t+Δt} , = [q_{v}]^{∗} + [q_{c}]^{∗}  (3.56)
[q_{c}]^{t+Δt} = 0  (3.57)

とし, 繰り返しを中止する.

圧平衡定数を用いる場合

硫化アンモニウムの生成反応

NH_{3} + H_{2}S → NH_{4}SH  (3.58)

のような, 2 種類の気体 1 モルづつから凝縮物質 1 モルが生成されるような生成反応の場合の, 湿潤飽和調節法を考える.

硫化アンモニウムの生成反応の圧平衡定数は,

K_{p} ≡ ln(p_{NH_{3}} · p_{H_{2}S}) = 61.781 − \frac{10834}{T} − ln 10^{2}  (3.59)

である. 圧平衡定数を用いることで, 任意の温度に対するアンモニアと硫化水素のモル比の積を求めることができる.

任意の温度 T における NH_{4}SH の生成量を X とすると, 圧平衡定数の式は以下のように書ける.

(p_{NH_{3}} − X)(p_{H_{2}S} − X) = e^{k_{p}}
X^{2} − (p_{NH_{3}} + p_{H_{2}S})X + p_{NH_{3}} · p_{H_{2}S} − e^{k_{p}} = 0  (3.60)

解の公式を使うと, 生成量 X は以下となる.

X = \frac{1}{2} \{(p_{NH_{3}} + p_{H_{2}S}) ± \sqrt{(p_{NH_{3}} + p_{H_{2}S})^{2} − 4(p_{NH_{3}} · p_{H_{2}S} − e^{K_{p}})}\}

X = \frac{1}{2} \{(p_{NH_{3}} + p_{H_{2}S}) ± \sqrt{(p_{NH_{3}} − p_{H_{2}S})^{2} + 4e^{K_{p}}}\}  (3.61)

根号の符号は exp (K_{p}) ≈ 0 の場合にとりうる X の値を仮定することで決める. exp (K_{p}) ≈ 0 の場合, 明らかに

X = min(P_{NH_{3}}, P_{H_{2}S})  (3.62)

である. ここで木星大気を想定し, P_{NH_{3}} > P_{H_{2}S} であることを仮定すると X = P_{H_{2}S} である. そこで (3.61) の根号の符号は exp (K_{p}) ≈ 0 のとき X = P_{H_{2}S} となるよう, 負を選択する.

X = \frac{1}{2} \{(p_{NH_{3}} + p_{H_{2}S}) − \sqrt{(p_{NH_{3}} − p_{H_{2}S})^{2} + 4e^{K_{p}}}.\}  (3.63)

X の満たすべき条件は,

0 ≤ X ≤ min(P_{NH_{3}}, P_{H_{2}S})  (3.64)

である. 上記の条件を満たさない場合には X = 0 とする.

X が (3.64) 式の条件を満たすならば, 次式を用いて暫定的に θ, q_{v}, q_{c} を求める.

[q_{NH_{3}}]^{t+Δt} = [q_{NH_{3}}]^{∗} + Δq_{NH_{3}},  (3.65)
[q_{H_{2}S}]^{t+Δt} = [q_{H_{2}S}]^{∗} + Δq_{H_{2}S},  (3.66)
[q_{NH_{4}SH}]^{t+Δt} = [q_{NH_{3}}]^{∗} + [q_{H_{2}S}]^{∗} − [q_{NH_{3}}]^{t+Δt} − [q_{H_{2}S}]^{t+Δt},  (3.67)
[θ]^{t+Δt} = θ^{∗} + γ ([q_{NH_{3}}]^{∗} + [q_{H_{2}S}]^{∗} − [q_{NH_{3}}]^{t+Δt} − [q_{H_{2}S}]^{t+Δt}) . (3.68)

ただし, γ = L_{NH_{4}SH}/(c_{p}_{d}Π) であり, Δq_{NH_{3}} と Δq_{H_{2}S} はそれぞれ, 生成量 X に対応する NH_{3} と H_{2}S の混合比である. 温位が収束するまで反復改良を行う.

3.3 乱流運動エネルギーの式

Klemp and Wilhelmson (1978) および CReSS (坪木と榊原篤志, 2001) と同様に, 1.5 次のクロージャーを用いる. 乱流エネルギーの時間発展方程式をリープフロッグ法を用いて時間方向に離散化すると, 以下のようになる.

[E_{turb}]_{i,k}^{t+Δt} = [E_{turb}]_{i,k}^{t−Δt} + 2Δt[F_{E_{turb}}]_{i,k}^{t}  (3.69)

ここで,

[F_{E_{turb}}]_{i,k}^{t} = −[Adv.E_{turb}]_{i,k}^{t} + [Buoy.E_{turb}]_{i,k}^{t−Δt} + [Shear.E_{turb}]_{i,k}^{t−Δt}
+[Turb.E_{turb}]_{i,k}^{t−Δt} + [Disp.E_{turb}]_{i,k}^{t−Δt}  (3.70)

である. CReSS にならい, 移流項を t で, 移流項以外を t − Δt で評価した.

F_{E_{turb}} に含まれる各項は以下のように書き下すことができる.

[Adv.E_{turb}]_{i,k}^{t} = \{u_{i(u),k} (\frac{∂E_{turb}}{∂x} )_{i(u),k}\}_{i,k}^{t} + \{w_{i,k(w)} (\frac{∂E_{turb}}{∂z} )_{i,k(w)}\}_{i,k}^{t} (3.71)

[Buoy.E_{turb}]_{i,k}^{t−Δt} = − \{\frac{3gC_{m}l}{2\overline{θ_{v}}}E_{turb_{i,k}} (\frac{∂θ_{el}}{∂z} )_{i,k(w)}\}_{i,k}^{t−Δt}  (3.72)

[Shear.K_{m}]_{i,k}^{t−Δt} = 2C_{m}l[E_{turb}^{1/2}]_{i,k}^{t−Δt} \{[(\frac{∂u}{∂x})^{2}]_{i,k}^{t−Δt} + [(\frac{∂w}{∂z} )^{2}]_{i,k}^{t−Δt}\}
+C_{m}l[E_{turb}^{1/2}]_{i,k}^{t−Δt}
[\{(\frac{∂u}{∂z})_{i(u),k(w)}\}_{i,k}^{t−Δt} + \{(\frac{∂w}{∂x} )_{i(u),k(w)}\}_{i,k}^{t−Δt}]^{2}

−\frac{2}{3}[E_{turb}]_{i,k}^{t−Δt} \{(\frac{∂u}{∂x})_{i,k}^{t−Δt} + (\frac{∂w}{∂z} )_{i,k}^{t−Δt}\}  (3.73)

[Turb.E_{turb}]_{i,k}^{t−Δt} = C_{m}l \{\frac{∂}{∂x} [[E_{turb}^{1/2}]_{i,k} (\frac{∂E_{turb}}{∂x} )_{i(u),k}]\}_{i,k}^{t−Δt}

+C_{m}l \{\frac{∂}{∂z} [[E_{turb}^{1/2}]_{i,k} (\frac{∂E_{turb}}{∂z} )_{i,k(w)}]\}_{i,k}^{t−Δt}  (3.74)
[Disp.E_{turb}]_{i,k}^{t−Δt} = −\frac{C_{ε}}{l} (E_{turb}^{3/2})_{i,k}^{t−Δt}  (3.75)

ここで C_{ε} = C_{m} = 0.2, 混合距離 l = (ΔxΔz)^{1/2} とする. また θ_{el} は以下で与えられる.

θ_{el} = \overline{θ_{v}} + θv^{′} (for q_{c} = 0)  (3.76)

θ_{el} = \overline{θ_{v}} + θv^{′} + \frac{∑ Lq_{v}}{c_{p}_{d}π̄} (for q_{c} > 0)  (3.77)

ただし,

\overline{θ_{v}} + θv^{′} = (\frac{⟨M⟩}{M_{d}} \overline{q_{d}} + ∑ \frac{⟨M⟩}{M_{v}} \overline{q_{v}}) (\overline{θ} + θ) + (\frac{⟨M⟩}{M_{d}}q_{d} + ∑ \frac{⟨M⟩}{M_{v}}q_{v}) \overline{θ} (3.78)

である.

3.4 時間フィルター

リープフロッグ法を用いたことによって生じる計算モードの増幅を抑制するため, Asselin (1972) の時間フィルターを長い時間刻みで 1 ステップ計算する毎に (実際には短い時間刻みの計算を N_{τ} ≡ 2Δt/Δτ ステップ計算する毎に) 適用する.

たとえば (3.6) を用いて u_{i(u),k}^{t+Δt} を計算する場合, 以下のように時間フィルターを適用する.

u_{i(u),k}^{∗} = u_{i(u),k}^{τ+(N_{τ} −1)Δτ} − [c̄_{p}θ̄_{v}Δτ \{\frac{∂π^{τ+(N_{τ} −1)Δτ}}{∂x} − \frac{∂(αDiv)^{τ+(N_{τ} −1)Δτ}}{∂x} \}]_{i(u),k}
+F_{u,i(u),k}^{t}Δτ,
u_{i(u),k}^{t+Δt} = (1 − 2γ)u_{i(u),k}^{t} + γ(u_{i(u),k}^{∗} + u_{i(u),k}^{t−Δt})  (3.79)

ここで γ はフィルターの係数であり, その値は 0.05 を用いる. (3.7), (3.8) の計算に対しても同様に時間フィルターを適用する.

3.5 スポンジ層

境界面付近での波の反射を抑えるために, 基礎方程式の付加的な項を付け加える.

\frac{∂ϕ}{∂t} = −Adv.ϕ + ⋯ + γ_{h}(x)(ϕ − ϕ_{e}) + γ_{v}(z)(ϕ − ϕ_{e})  (3.80)

ただし, ϕ は任意の予報変数であり, ϕ_{e} は客観解析値等の既知の値である. この項は 1 つ前のタイムステップ t − Δt で計算され, 小さいタイムステップで扱われる予報変数に対しても, 移流項や数値粘性項と同様に 2Δt の大きなタイムステップ間の値として評価される。具体的には,

[π]^{t+Δt} = 2Δt {[Adv.π]^{t} + ⋯ + {γ_{h}(x) + γ_{v}(z)} (π − π̄)^{t−Δt}} (3.81)
[u]^{t+Δt} = 2Δt {[Adv.u]^{t} + ⋯ + {γ_{h}(x) + γ_{v}(z)} [u]^{t−Δt}}  (3.82)
[w]^{t+Δt} = 2Δt {[Adv.w]^{t} + ⋯ + {γ_{h}(x) + γ_{v}(z)} [w]^{t−Δt}}  (3.83)
[θ]^{t+Δt} = 2Δt {[Adv.θ]^{t} + ⋯ + {γ_{h}(x) + γ_{v}(z)} [θ]^{t−Δt}}  (3.84)

とする. 但し π̄ はエクスナー関数の基本場である.

γ_{h}, γ_{v} はそれぞれ水平方向には各境界面に向かって, 鉛直方向には上境界面に向かって小さくなる減衰係数である. これらの減衰係数は, 水平方向には吸収層の厚みを d_{h} とし, x の範囲を 0 ≤ x ≤ x_{max} とすれば,

γ_{h} = α_{h} (1 − \frac{x}{d_{h}} )^{3} (x < d_{h}),

γ_{h} = 0 (d_{h} ≤ x ≤ x_{max} − d_{h}),

γ_{h} = α_{h} (1 − \frac{(x_{max} − x)}{d_{h}})^{3} (x > x_{max} − d_{h}),  (3.85)

であり, 鉛直方向には吸収層の厚さを d_{v} とし, z の範囲を 0 ≤ z ≤ z_{max} とすれば,

γ_{v} = 0 (z ≤ z_{max} − d_{v}),
γ_{v} = α_{v} (1 − cos \frac{π(z − z_{max} − d_{v})}{d_{v}})^{3} (z > z_{max} − d_{v}),  (3.86)

である. ここで, α_{h}, α_{v} はそれぞれ水平・鉛直方向の減衰定数である. α_{h}, α_{v} は時間の逆数の次元を持ち, それらの逆数 1/α_{h}, 1/α_{v} は e-folding time と呼ばれる. e-folding time は通常 100 – 300 s に設定する. また吸収層の厚み d_{h}, d_{v} はそれぞれ, 水平方向には数格子分, 鉛直方向には上面から 1/3 程度設定すれば良い.

\underline{2 次元非静力学モデルの離散化第 4 章参考文献 28}

第4章参考文献

Asselin, R. A., 1972: Frequency filter for time integrations. Mon. Wea. Rev., 100, 487–490.

Klemp J. B. and R. B. Wilhelmson, 1978: The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070-1096.

Skamarock, W. C., and J. B. Klemp, 1992: The stability of time-split numerical methods for the hydrostatic and the nonhydrostatic elastic equations, Mon. Wea. Rev., 120, 2109–2127.

坪木和久, 榊原篤志, 2001: CReSS ユーザーズガイド第 2 版.

http://www.tokyo.rist.or.jp/CReSS Fujin/CReSS.top.html \underline{{2 次元非静力学モデルの離散化 \atop 付録 A 圧力方程式 (3.10) の左辺の空間微分の書き下し 29}}

付録A {圧力方程式 (3.10) の左辺 \atop の空間微分の書き下し}

(3.10) 左辺の変形を行う.

(3.10)left side = Π_{i,k}^{τ+Δτ}

−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{k(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,k(w)}\}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{k−1(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,k−1(w)}\}

= Π_{i,k}^{τ+Δτ}

−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{k(w)} (\frac{Π_{i,k+1}^{τ+Δτ} − Π_{i,k}^{τ+Δτ}}{Δz} )\}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{k−1(w)} (\frac{Π_{i,k}^{τ+Δτ} − Π_{i,k−1}^{τ+Δτ}}{Δz} )\}

= \{−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{k(w)}\} Π_{i,k+1}^{τ+Δτ}

+ [1 + β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz^{2}} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{k(w)} + (c̄_{p}ρ̄θ̄_{v}^{2})_{k−1(w)}\}] Π_{i,k}^{τ+Δτ}

+ \{−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{k} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{k−1(w)}\} Π_{i,k−1}^{τ+Δτ}.  (A.1)

A.1 下部境界

下部境界 (k(w) = 0(w)) について考える. この時 (3.7) 式は,

β (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,0(w)} = [(\frac{∂(αDiv)^{τ}}{∂z} ) − (1 − β) (\frac{∂Π^{τ}}{∂z} ) + (\frac{1}{c̄_{p}θ̄_{v}} F_{w}^{t})]_{i,0(w)}
≡ E_{i,0(w)}  (A.2)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{{2 次元非静力学モデルの離散化 \atop 付録 A 圧力方程式 (3.10) の左辺の空間微分の書き下し 30}} となるので, (3.10) 式の左辺は, k = 1 の場合には,

(3.10)left side = Π_{i,1}^{τ+Δτ}

−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{1(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,1(w)}\}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{0(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,0(w)}\}

= Π_{i,1}^{τ+Δτ}

−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{1(w)} (\frac{Π_{i,2}^{τ+Δτ} − Π_{i,1}^{τ+Δτ}}{Δz} )\}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{0(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,0(w)}\}

= \{−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{1(w)}\} Π_{i,2}^{τ+Δτ}

+ \{1 + β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{1(w)}\} Π_{i,1}^{τ+Δτ}

+β (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{1} \frac{1}{Δz} (c̄_{p}ρ̄θ̄_{v}^{2})_{0(w)} E_{i,0(w)}

A.2 上部境界

上部境界 (k(w) = km(w)) について考える. (3.10) 式の左辺は,

β (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,km(w)} = [(\frac{∂(αDiv)^{τ}}{∂z} ) − (1 − β) (\frac{∂Π^{τ}}{∂z} ) + (\frac{1}{c̄_{p}θ̄_{v}} F_{w}^{t})]_{i,km(w)}
≡ E_{i,km(w)}  (A.3)

となるので, (3.10) 式の左辺は, k(w) = km(w) の場合には,

(3.10)left side = Π_{i,km}^{τ+Δτ}

−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{km(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,km(w)}\}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{km−1(w)} (\frac{∂Π^{τ+Δτ}}{∂z} )_{i,km−1(w)}\}

= Π_{i,km}^{τ+Δτ}

−β (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz} (c̄_{p}ρ̄θ̄_{v}^{2})_{km(w)} E_{i,km(w)}

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{{2 次元非静力学モデルの離散化 \atop 付録 A 圧力方程式 (3.10) の左辺の空間微分の書き下し 31}}

+β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz} \{(c̄_{p}ρ̄θ̄_{v}^{2})_{km−1(w)} (\frac{Π_{i,km}^{τ+Δτ} − Π_{i,km−1}^{τ+Δτ}}{Δz} )\}

= \{1 + β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{km−1(w)}\} Π_{i,km}^{τ+Δτ}

+ \{−β^{2} (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz^{2}} (c̄_{p}ρ̄θ̄_{v}^{2})_{km−1(w)}\} Π_{i,km−1}^{τ+Δτ}

−β (\begin{matrix}c̄^{2} & Δτ & 2 \\ c̄_{p} & ρ̄ & θ̄_{v}^{2} \end{matrix} )_{km} \frac{1}{Δz} (c̄_{p}ρ̄θ̄_{v}^{2})_{km(w)} E_{i,km(w)}

\underline{2 次元非静力学モデルの離散化付録 B 音波減衰項について 32}

付録B 音波減衰項について

本モデルで用いている時間方向の離散化方法 (時間分割法) は, Klemp and Wilhelmson (1978) によって最初に提案された方法である. この方法を用いると, 音速および移流に対する CFL 条件をそれぞれ満たしている場合でも計算不安定を起こす場合がある. この影響は短い時間間隔で積分するステップ数 2Δt/Δτ を増加させていくにつれ大きくなる (Skamarock and Klemp, 1992, 図 1 を参照). Skamarock and Klemp (1992) はこの計算不安定を回避するため, 音波を選択的に減衰させるフィルターとして運動方程式に音波減衰項を導入することを提案した.

音波減衰項は速度場の発散に対する拡散として作用する. このことは音波減衰項を含む線形化された基礎方程式

\frac{∂u}{∂t} = −\overline{c}_{p}\overline{θ}_{v} \frac{∂Π}{∂x} + α\frac{∂D}{∂x},  (B.1)
\frac{∂w}{∂t} = −\overline{c}_{p}\overline{θ}_{v} \frac{∂Π}{∂z} + α\frac{∂D}{∂z} + g \frac{θ}{\overline{θ}_{v}} ,  (B.2)
\frac{∂Π}{∂t} = −\begin{matrix}\overline{c} & 2 \\ \overline{c} & p & \overline{ρ} & \overline{θ} & {2 \atop v} \end{matrix} [\frac{∂\overline{ρ}\overline{θ}_{v}u}{∂x} + \frac{∂\overline{ρ}\overline{θ}_{v}w}{∂z} ] ,  (B.3)
\frac{∂θ}{∂t} = −w\frac{∂\overline{θ}_{v}}{∂z} ,  (B.4)
D = \frac{∂u}{∂x} + \frac{∂w}{∂z}

から発散方程式

\frac{∂D}{∂t} = −\overline{c}_{p}\overline{θ}_{v}∇^{2}Π + g \frac{∂}{∂z} (\frac{θ}{\overline{θ}_{v}} ) + α∇^{2}D,  (B.5)
∇^{2} = \frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂z^{2}}

を求めることで理解することができる.

音波減衰項は速度場の発散に対する拡散なので, 重力波成分も拡散させる可能性が考えられる. しかし Skamarock and Klemp (1992) は線形化された基礎方程式の分散関係を用いて, 音波減衰項の重力波成分への影響は小さいことを示している.

以下では Skamarock and Klemp (1992) における議論の概要を示す. 簡単のため圧力方程式に現れる基本場の量は定数とし, 全ての変数を u = ûe^{i(kx+lz−ωt)} のような解を持つと仮定して分散関係を求めると,

ω^{4} + iα(k^{2} + l^{2})ω^{3} − [\overline{c}^{2}(k^{2} + l^{2}) + N^{2}]ω^{2} − iαk^{2}N^{2}ω + \overline{c}^{2}k^{2}N^{2} = 0, (B.6)

となる. ここで N^{2} = \frac{g}{\overline{θ}_{v}} \frac{∂\overline{θ}_{v}}{∂z} である. 純粋な内部重力波の分散関係式

ω^{2} = \frac{N^{2}k^{2}}{k^{2} + l^{2}} ,

を α を含む項に代入すると α を含む項は互いに打ち消しあうことから, 音波減衰項の重力波成分への影響は小さいであろうと予想される. 実際に, ϵ ≡ α\sqrt{k^{2} + l^{2}}/\overline{c} が小さいと仮定して ω を ϵ の巾で展開して音波減衰項の影響を評価すると, 音波減衰項は重力波成分を減衰させる方向にはたらき, その減衰率は非常に小さいことがわかる (Skamarock and Klemp, 1992, 図 5 参照).

音波減衰項の係数 α の値を決めるにあたり, 考慮しなければならないことが 2 つある. 1 つは音波減衰項そのものが計算不安定の原因とならないようすることであり, もう 1 つは音波減衰項の重力波成分への影響が大きくならないようすることである. 前者は (B.5) に示したように, 音波減衰項が発散方程式の拡散項となることから要請される. 音波減衰項の時間積分は前進差分を用いて行われるので ((3.6), (3.7) を参照), 計算不安定を起こさないためには拡散項を前進差分で時間積分する場合の安定性条件

\frac{αΔτ}{Min(Δx^{2}, Δz^{2})} ≤ \frac{1}{2}  (B.7)

を満たさなければならない. 後者は前段落で紹介した Skamarock and Klemp (1992) の議論で用いた ϵ が小さいという仮定から要請される^{1}. ϵ の最大値は

Max(ϵ) = \frac{2α}{Min(Δx, Δz)\overline{c}}
= \frac{2αΔτ}{Min(Δx^{2}, Δz^{2})} /\frac{\overline{c}Δτ}{Min(Δx, Δz)}

と与えられるので, Max(ϵ) ≤ 1 とするためには

\frac{αΔτ}{Min(Δx^{2}, Δz^{2})} ≤ \frac{1}{2} \frac{\overline{c}Δτ}{Min(Δx, Δz)}  (B.8)

でなければならない.(B.8) 右辺に現れる音速に対するクーラン数は 1 より小さい値とするので, (B.8) の条件を考慮すれば (B.7) は自動的に満たされることになる.

1この条件に関する議論は Skamarock and Klemp (1992) ではなされていない \overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化付録 C 差分式の導出と誤差 34}

付録C 差分式の導出と誤差

ここでは交互格子を用いた場合の空間微分の差分式の導出と, その誤差についてまとめる. 具体例としてフラックス格子点の変数 u の, (i, j) 格子点上における x 方向一階微分

[\frac{∂u}{∂x}]_{i,k}

の差分式と, その誤差を考える.

C.1 2 次精度中心差分

フラックス格子点 (i(u), j) 上の u (u_{i(u),j}) を x 方向に −Δx/2 だけずれたスカラー格子点 (i, j) 上の u (u_{i,j}) のテーラー展開として表すと, 以下のようになる.

u_{i(u),j} = u_{i,j} + [\frac{∂u}{∂x}]_{i,k} \frac{Δx}{2} + \frac{1}{2!} [\frac{∂^{2}u}{∂x^{2}} ]_{i,k} (\frac{Δx}{2})^{2} + \frac{1}{3!} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{Δx}{2})^{3}
+\frac{1}{4!} [\frac{∂^{4}u}{∂x^{4}} ]_{i,k} (\frac{Δx}{2})^{4} + \frac{1}{5!} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{Δx}{2})^{5} + · · · .  (C.1)

同様に, フラックス格子点 (i − 1(u), j) 上の u (u_{i−1(u),j}) を u_{i,j} のテーラー展開として表すと, 以下のようになる.

u_{i−1(u),j} = u_{i,j} − [\frac{∂u}{∂x}]_{i,k} \frac{Δx}{2} + \frac{1}{2!} [\frac{∂^{2}u}{∂x^{2}} ]_{i,k} (\frac{Δx}{2})^{2} − \frac{1}{3!} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{Δx}{2})^{3}
+\frac{1}{4!} [\frac{∂^{4}u}{∂x^{4}} ]_{i,k} (\frac{Δx}{2})^{4} − \frac{1}{5!} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{Δx}{2})^{5} + · · · .  (C.2)

(C.1) − (C.2) より,

u_{i(u),j} − u_{i−1(u),j} = [\frac{∂u}{∂x}]_{i,k} Δx + \frac{1}{3} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{Δx}{2})^{3} + \frac{1}{60} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{Δx}{2})^{5}
+O[(Δx)^{7}]  (C.3)

これを変形すると (i, j) 格子点上における u の x 方向一階微分の式が得られる.

[\frac{∂π}{∂x} ]_{i(u),k} = \frac{π_{i+1,j} − π_{i,j}}{Δx} − \frac{1}{24} [\frac{∂^{3}π}{∂x^{3}} ]_{i(u),k} (Δx)^{2}
−\frac{1}{1920} [\frac{∂^{5}π}{∂x^{5}} ]_{i(u),k} (Δx)^{4} + O[(Δx)^{6}]  (C.4)

上式の (Δx)^{2} 以上の高次項を無視することで, 交互格子を用いた場合の 2 次精度中心差分の式

[\frac{∂u}{∂x}]_{i(u),k} = \frac{u_{i+1,j} − u_{i,j}}{Δx}  (C.5)

が得られる. このときの誤差の大きさは

|\frac{1}{24} [\frac{∂^{3}u}{∂x^{3}} ]_{i(u),k} (Δx)^{2}|  (C.6)

となる.

C.2 4 次精度中心差分

2 次精度中心差分の式を求める際に用いた (C.1), (C.2) に加え, (i, j) から x 方向に ±Δ3x/2 だけずれたフラックス格子点での u (u_{i+1(u),j}, u_{i−2(u),j}) の値を u_{i,j} のテーラー展開として求める.

u_{i+1(u),j} = u_{i,j} + [\frac{∂u}{∂x}]_{i,k} \frac{3Δx}{2} + \frac{1}{2!} [\frac{∂^{2}u}{∂x^{2}} ]_{i,k} (\frac{3Δx}{2})^{2} + \frac{1}{3!} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{3Δx}{2})^{3}
+\frac{1}{4!} [\frac{∂^{4}u}{∂x^{4}} ]_{i,k} (\frac{3Δx}{2})^{4} + \frac{1}{5!} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{3Δx}{2})^{5} + · · · .  (C.7)
u_{i−2(u),j} = u_{i,j} − [\frac{∂u}{∂x}]_{i,k} \frac{3Δx}{2} + \frac{1}{2!} [\frac{∂^{2}u}{∂x^{2}} ]_{i,k} (\frac{3Δx}{2})^{2} − \frac{1}{3!} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{3Δx}{2})^{3}
+\frac{1}{4!} [\frac{∂^{4}u}{∂x^{4}} ]_{i,k} (\frac{3Δx}{2})^{4} − \frac{1}{5!} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{3Δx}{2})^{5} + · · · .  (C.8)

(C.7) − (C.8) より,

u_{i+1(u),j} − u_{i−2(u),j} = [\frac{∂u}{∂x}]_{i,k} 3Δx + \frac{1}{3} [\frac{∂^{3}u}{∂x^{3}} ]_{i,k} (\frac{3Δx}{2})^{3} + \frac{1}{60} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{3Δx}{2})^{5}
+O[(Δx)^{7}]  (C.9)

\underline{2 次元非静力学モデルの離散化付録 C 差分式の導出と誤差 36} (C.3)×27−(C.9) を行い (Δx)^{3} の項を消去すると,

27(u_{i(u),j} − u_{i−1(u),j}) − (u_{i+1(u),j} − u_{i−2(u),j}) = [\frac{∂u}{∂x}]_{i,k} 24Δx − \frac{216}{60} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (\frac{Δx}{2})^{5}
+O[(Δx)^{7}]  (C.10)

これを変形して (i, j) 格子点上における u の x 方向一階微分の式が得られる.

[\frac{∂u}{∂x}]_{i,k} = \frac{9}{8} (\frac{u_{i,j} − u_{i−1,j}}{Δx} ) − \frac{1}{24} (\frac{u_{i+1,j} − u_{i−2,j}}{Δx} )
+\frac{3}{640} [\frac{∂^{5}u}{∂x^{5}} ]_{i,k} (Δx)^{4} + O[(Δx)^{6}]  (C.11)

上式の (Δx)^{4} 以上の高次項を無視することで, 交互格子を用いた場合の 4 次精度中心差分式

[\frac{∂u}{∂x}]_{i,k} = \frac{9}{8} (\frac{u_{i,j} − u_{i−1,j}}{Δx} ) − \frac{1}{24} (\frac{u_{i+1,j} − u_{i−2,j}}{Δx} )  (C.12)

が得られる. このときの誤差の大きさは

|\frac{3}{640} [\frac{∂^{5}u}{∂x^{5}} ]_{i(u),k} (Δx)^{4}|  (C.13)

となる.

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化付録 D 化学物性値の計算法 37}

付録D 化学物性値の計算法

D.1 相平衡条件

D.1.1 飽和蒸気圧

飽和蒸気圧は, Antoine の式より求める.

ln p^{∗} = (A − B/(C + T − T_{0})) ∗ log(10.0) + ln(133.322)  (D.1)

ここで, p^{∗} は飽和蒸気圧, T は温度, T_{0} = 273.15 である. A, B, C は Antoine 係数である. それらの値は化学便覧改訂 4 版から得る. 化学便覧改訂 4 版では, 圧力の単位が mmHg, 温度の単位が^{∘}C であるので, 単位の換算項が付加されている.

表 D.1: 水, アンモニアの Antoine 係数

A  B  C
H_{2}O(l) 7.9186968 1636.909 224.92
H_{2}O(s) 8.184254 1791.3 238.1
NH_{3}(s) 9.96382 1617.907 272.55

任意の温度が与えられた場合, 凝縮量は飽和蒸気圧と分圧の差として見積もることができる.

D.1.2 圧平衡定数

硫化アンモニウムの生成反応

NH_{3} + H_{2}S → NH_{4}SH  (D.2)

の圧平衡定数は,

K_{p} = ln(p_{NH_{3}} · p_{H_{2}S}) = 61.781 − \frac{10834}{T} − ln 10^{2}  (D.3)

である. 圧平衡定数を用いることで, 任意の温度に対するアンモニアと硫化水素のモル比の積を求めることができる.

D.2 生成のエンタルピー変化

D.2.1 潜熱

飽和蒸気圧と潜熱はクラウジウス・クラペイロンの式,

\frac{dp_{v}}{dT} = \frac{p_{v}L_{v}}{RT^{2}}  (D.4)

で関係づけられる. この式を L_{v} の式としてまとめなおすことで, 潜熱は以下のように与えられる.

L_{v} = \frac{d ln p_{v}}{dT} R_{v}T^{2}  (D.5)

但し R_{v} は凝縮成分に対する気体定数である. Antoine の式を代入すると,

L_{v} = \{\frac{B ln(10.0)}{(C + T − T_{0})^{2}} \} R_{v}T^{2}  (D.6)

である.

D.2.2 反応熱

硫化アンモニウムの生成反応

NH_{3} + H_{2}S → NH_{4}SH  (D.7)

において, NH_{4}SH のエントロピーと NH_{3} と H_{2}S のエントロピーの差が, 反応に伴うエントロピー変化に対応する.

NH_{4}SH のモルエントロピーは,

s_{NH_{4}SH} = −\frac{∂μ_{NH_{4}SH}}{∂T}
= −\frac{∂}{∂T} (μ_{NH_{3}}^{∘} + μ_{H_{2}S}^{∘} + RT K_{p} − RT ln p_{0}^{2})
= s_{NH_{3}}^{∘} + s_{H_{2}S}^{∘} − RT \frac{∂K_{p}}{∂T} − RK_{p} + R ln p_{0}^{2}  (D.8)

である. ここで μ_{NH_{3}}^{∘}, μ_{H_{2}S}^{∘} は NH_{3} と N_{2}S の標準化学ポテンシャル, s_{NH_{3}}^{∘}, s_{H_{2}S}^{∘} はそれに対応するエントロピー, K_{p} は (D.7) の反応式の圧平衡定数である. NH_{3} と H_{2}S のモルエントロピーの和は,

s_{NH_{3}} + s_{H_{2}S} = s_{NH_{3}}^{∘} + s_{H_{2}S}^{∘} − R ln(p_{NH_{3}} · p_{H_{2}S}) + R ln p_{0}^{2}
= s_{NH_{3}}^{∘} + s_{H_{2}S}^{∘} − R ln K_{p} + R ln p_{0}^{2}  (D.9)

(D.8) と (D.9) の差

Δs = RT \frac{∂K_{p}}{∂T}  (D.10)

が反応のエントロピー変化に相当する. モル当たりの反応熱は,

L_{NH_{4}SH} = T Δs = RT^{2} \frac{∂K_{p}}{∂T}  (D.11)

である. NH_{4}SH 生成反応の圧平衡定数を代入すると,

L_{NH_{4}SH} = \frac{10834}{T^{2}}RT^{2} = 10834R  (D.12)

である.

\underline{2 次元非静力学モデルの離散化付録 E 単位の換算等の計算 40}

付録E 単位の換算等の計算

E.1 混合比とモル比 (分圧) との換算式

混合気体を考える. 乾燥成分の添え字を d とし, 湿潤成分の添え字を v とする. ここでは任意の湿潤成分気体 (添え字 x) の混合比, 分圧を求める. 混合比 q_{x} を分圧 p_{x} で表現すると,

q_{x} = \frac{ρ_{x}}{ρ_{d}}
= (\frac{R_{d}T}{p} \frac{1/M_{d} + ∑_{v} q_{v}/M_{v}}{1/M_{d}} ) \frac{p_{x}}{R_{x}T}
= \frac{p_{x}}{p} M_{x}(1/M_{d} + ∑_{v} q_{v}/M_{v})
q_{x} (1 − \frac{p_{x}}{p} ) = \frac{p_{x}}{p} M_{x}(1/M_{d} + ∑_{v≠x} q_{v}/M_{v})
q_{x} = \frac{p_{x}}{p} \frac{M_{x}(1/M_{d} + ∑_{v≠x} q_{v}/M_{v})}{(1 − \frac{p_{x}}{p} )}
q_{x} = \frac{p_{x}M_{x}(1/M_{d} + ∑_{v≠x} q_{v}/M_{v})}{p − p_{x}}
q_{x} = \frac{p_{x}M_{x}(1/M_{d} + ∑_{v≠x} q_{v}/M_{v})}{p_{0}π^{c_{p}_{d}/R_{d}} − p_{x}}  (E.1)

となり, モル比 x_{x} で表現すると,

q_{x} = \frac{x_{x}M_{x}(1/M_{d} + ∑_{v≠x} q_{v}/M_{v})}{1 − x_{x}}  (E.2)

となる. 一方で, 分圧 p_{x} を混合比 q_{x} で表現すると,

p_{x} = \frac{q_{x}p}{M_{x}(1/M_{d} + ∑_{v} q_{v}/M_{v})}
p_{x} = \frac{q_{x}p_{0}π^{c_{p}_{d}/R_{d}}}{M_{x}(1/M_{d} + ∑_{v} q_{v}/M_{v})}  (E.3)

\overline{main.tex 2010 年 3 月 17 日 (山下達也, 杉山耕一朗, 石渡正樹)}

\underline{2 次元非静力学モデルの離散化付録 E 単位の換算等の計算 41} であり, モル比 x_{x} で表現すると,

x_{x} = \frac{q_{x}}{M_{x}(1/M_{d} + ∑_{v} q_{v}/M_{v})}  (E.4)

である.

\underline{2 次元非静力学モデルの離散化謝辞 42}

謝辞

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

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

において公開されているものである. \begin{matrix}c \\ ◯ \end{matrix}山下達也, 高橋こう子, 杉山耕一朗, 小高正嗣, 中島健介, 林祥介 (T. Yamashita, K. Takahashi, K. Sugiyama, M. Odaka, K. Nakajima and Y.-Y. Hayashi) 2003. 本資源は, 著作者の諸権利に抵触しない (迷惑をかけない) 限りにおいて自由に利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめることをお願いする (無保証無責任原則).

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

dcstaff@gfd-dennou.org

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