next up previous
: 3 乱流モデル : Two dimensional anelastic model : 1 モデル離散化の概要


2 大気モデル

2.1 運動方程式

第 1 部に示した式(1)$\sim$(3)は以下のように変形した後に離散化する.

\begin{eqnarray*}
\DP{u}{t} &=& -\DP{\hat{P}}{x} + \alpha , \\
\DP{v}{t} &=& \beta , \\
\DP{w}{t} &=& -\DP{\hat{P}}{z}
\end{eqnarray*}

ただし,

\begin{eqnarray*}
\hat{P} &\equiv& c_{p}\Theta _{0}\pi - \gamma , \\
\alpha &...
... - w\DP{w}{z} + g\frac{\theta}{\Theta _{0}}
+ D(w)\right) \Dd z
\end{eqnarray*}

である. このように変形するのはモデル上下端において $\pi $ の境界条件を時 間変化させないようにするためである.

これらを以下のように離散化する. 移流項 $\mbox{D[UVW]ADV}$ は保存型と移流 型の混合型で計算する. 摩擦項 $[\mbox{D[UVW]VIS}]_{i+\frac{1}{2},j}^{N}$, $[\mbox{D[UVW]NLV}]_{i+\frac{1}{2},j}^{N}$ の時間積分は前進差分, その他の 項の積分は leap frog スキームと前進差分の組み合わせを用いて行う. 圧力項 $\hat{P}$ の導出について第2.3節を参照.

$\displaystyle u_{i+\frac{1}{2},j}^{n+1}$ $\textstyle =$ $\displaystyle u_{i+\frac{1}{2},j}^{N} + dt \left\{
\frac{\hat{P}_{i+1,j}-\hat{P}_{i,j}}{\Delta x}
+ \alpha _{i+\frac{1}{2},j} \right\},$ (1)
$\displaystyle v_{i+\frac{1}{2},j}^{n+1}$ $\textstyle =$ $\displaystyle v_{i+\frac{1}{2},j}^{N}
+ dt \beta _{i+\frac{1}{2},j}^{n},$ (2)
$\displaystyle w_{i,j+\frac{1}{2}}^{n+1}$ $\textstyle =$ $\displaystyle w_{i,j+\frac{1}{2}}^{N} + dt
\frac{\hat{P}_{i,j+1}-\hat{P}_{i,j}}{\Delta z_{j+\frac{1}{2}}}.$ (3)


\begin{displaymath}
N = \left\{
\begin{array}{lcl}
n-1 & \mbox{for} & \mbox{l...
... \Delta t & \mbox{for} & \mbox{forward}.
\end{array} \right.
\end{displaymath} (4)


$\displaystyle \alpha _{i+\frac{1}{2},j}$ $\textstyle =$ $\displaystyle - ( \gamma _{i+1,j} - \gamma _{i,j})
+ [\mbox{DUADV}]_{i+\frac{1}...
...
+ [\mbox{DUPRS}]_{i+\frac{1}{2},j}^{n}
+ [\mbox{DUCOLI}]_{i+\frac{1}{2},j}^{n}$  
    $\displaystyle + [\mbox{DUVIS}]_{i+\frac{1}{2},j}^{N}
+ [\mbox{DUNLV}]_{i+\frac{1}{2},j}^{N} ,$ (5)
$\displaystyle \beta _{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle [\mbox{DVADV}]_{i+\frac{1}{2},j}^{n}
+ [\mbox{DVCOLI}]_{i+\frac{1...
...
+ [\mbox{DVVIS}]_{i+\frac{1}{2},j}^{N}
+ [\mbox{DVNLV}]_{i+\frac{1}{2},j}^{N},$ (6)
$\displaystyle \gamma _{i,j}$ $\textstyle =$ $\displaystyle \sum _{j'=0}^{j}
\left(
[\mbox{DWADV}]_{i,j'-\frac{1}{2}}^{n}
+ [...
...{DWVIS}]_{i,j'-\frac{1}{2}}^{N}
+ [\mbox{DWNLV}]_{i,j'-\frac{1}{2}}^{N}
\right.$  
    $\displaystyle \left.
+ [\mbox{BUOY}]_{i,j'-\frac{1}{2}}^{n}
\right)\Delta z_{j'-\frac{1}{2}}$ (7)


$\displaystyle \mbox{DUADV}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle - \frac{1}{\Delta x}\left\{
\left(u_{i+1,j}^{n}u_{i+1,j}^{n} - u_{i-1,j}^{n}u_{i-1,j}^{n}\right)
\right.$  
    $\displaystyle +
\left.
\left(
\rho _{0,j+\frac{1}{2}}
u_{i+\frac{1}{2},j+\frac{...
...\frac{1}{2},j-\frac{1}{2}}^{n}
\right)
\Delta x\right/(\rho _{0,j}\Delta z_{j})$  
    $\displaystyle \left.
- u_{i+\frac{1}{2},j}^{n}
\left(\Ddiv \rho _{0}\Dvect{v}/\rho _{0}\right)
_{i+\frac{1}{2},j}^{n}
\right\},$ (8)
$\displaystyle \mbox{DVADV}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle - \left\{
\left(v_{i+1,j}^{n}u_{i+1,j}^{n} - v_{i-1,j}^{n}u_{i-1,j}^{n}\right)
\right.$  
    $\displaystyle +
\left.
\left(
\rho _{0,j+\frac{1}{2}}
v_{i+\frac{1}{2},j+\frac{...
...\frac{1}{2},j-\frac{1}{2}}^{n}
\right)
\Delta x\right/(\rho _{0,j}\Delta z_{j})$  
    $\displaystyle \left.
\left.
- v_{i+\frac{1}{2},j}^{n}
\left(\Ddiv \rho _{0}\Dvect{v}/\rho _{0}\right)_{i+\frac{1}{2},j}^{n}
\right\}
\right/ \Delta x,$ (9)
$\displaystyle \mbox{DWADV}_{i,j+\frac{1}{2}}^{n}$ $\textstyle =$ $\displaystyle - \frac{1}{\Delta x}\left\{
\left(w_{i+\frac{1}{2},j+\frac{1}{2}}...
...ac{1}{2},j+\frac{1}{2}}^{n}
u_{i-\frac{1}{2},j+\frac{1}{2}}^{n}
\right)
\right.$  
    $\displaystyle + \left.
\left( \rho _{0,j+1}w_{i,j+1}^{n}w_{i,j+1}^{n}
- \rho _{...
...^{n}
\right)
\Delta x \right/
(\rho _{0,j+\frac{1}{2}}\Delta z_{j+\frac{1}{2}})$  
    $\displaystyle \left.
- w_{i,j+\frac{1}{2}}^{n}
\left(\Ddiv \rho _{0}\Dvect{v}/\rho _{0,j}\right)
_{i,j+\frac{1}{2}}^{n}
\right\},$ (10)
       
$\displaystyle \mbox{DUVIS}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle \frac{1}{(\Delta x)^{2}}\left\{
\left[ K_{i+1,j}^{n}
\left(u_{i+\...
...n}
\left(u_{i+\frac{1}{2},j}^{n}-u_{i-\frac{1}{2},j}^{n}\right)
\right]
\right.$  
    $\displaystyle + \frac{(\Delta x)^{2}}{\rho _{0,j}\Delta z_{j}}
\left[ \left.
\r...
...1}^{n}-u_{i+\frac{1}{2},j}^{n}\right)
\right/\Delta z_{j+\frac{1}{2}} -
\right.$  
    $\displaystyle \left.
\left.
\left. \rho _{0,j-\frac{1}{2}}
K_{i+\frac{1}{2},j-\...
...{i+\frac{1}{2},j-1}^{n}\right)
\right/\Delta z_{j-\frac{1}{2}}
\right]\right\},$ (11)
$\displaystyle \mbox{DVVIS}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle \frac{1}{(\Delta x)^{2}}
\left\{
\left[ K_{i+1,j}^{n}
\left(v_{i+...
...n}
\left(v_{i+\frac{1}{2},j}^{n}-v_{i-\frac{1}{2},j}^{n}\right)
\right]
\right.$  
    $\displaystyle + \frac{(\Delta x)^{2}}{\rho _{0,j}\Delta z_{j}}
\left[ \left.
\r...
...1}^{n}-u_{i+\frac{1}{2},j}^{n}\right)
\right/\Delta z_{j+\frac{1}{2}} -
\right.$  
    $\displaystyle \left.
\left.
\left. \rho _{0,j-\frac{1}{2}}
K_{i+\frac{1}{2},j-\...
...{i+\frac{1}{2},j-1}^{n}\right)
\right/\Delta z_{j-\frac{1}{2}}
\right]\right\},$ (12)
$\displaystyle \mbox{DWVIS}_{i,j+\frac{1}{2}}^{n}$ $\textstyle =$ $\displaystyle \frac{1}{(\Delta x)^{2}}\left\{
\left[ K_{i+\frac{1}{2},j+\frac{1...
...
\left(w_{i,j+\frac{1}{2}}^{n}-w_{i-1,j+\frac{1}{2}}^{n}\right)
\right]
\right.$  
    $\displaystyle + \frac{(\Delta x)^{2}}{\rho _{0,j}\Delta z_{j}}
\left[ \left.
\r...
...j+1}^{n}
\left(w_{i,j+1}^{n}-w_{i,j}^{n}\right)
\right/\Delta z_{j+1} -
\right.$  
    $\displaystyle \left.
\left.
\left. \rho _{0,j}K_{i,j}^{n}
\left(w_{i,j+\frac{1}{2}}^{n}-w_{i,j-\frac{1}{2}}^{n}\right)
\right/\Delta z_{j-1}
\right] \right\}.$ (13)


$\displaystyle \mbox{DUNLV}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle \left\{
\left[
\left(u_{i+\frac{3}{2},j}^{n}-u_{i+\frac{1}{2},j}^...
...t(u_{i+\frac{1}{2},j+1}^{n}-u_{i+\frac{1}{2},j}^{n}
\right)^{3} \right. \right.$  
    $\displaystyle \left.
\left.
\left. - \left(u_{i+\frac{1}{2},j}^{n}-u_{i+\frac{1...
...]
\right\}
\right/
(16.0 \cdot 10^{3} \cdot \rho _{0,j}\Delta z_{j}/\Delta x ),$ (14)
$\displaystyle \mbox{DVNLV}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle \left\{
\left[
\left(v_{i+\frac{3}{2},j}^{n}-v_{i+\frac{1}{2},j}^...
...t(v_{i+\frac{1}{2},j+1}^{n}-v_{i+\frac{1}{2},j}^{n}
\right)^{3} \right.
\right.$  
    $\displaystyle \left.
\left.
\left. - \left(v_{i+\frac{1}{2},j}^{n}-v_{i+\frac{1...
...]
\right\}
\right/
(16.0 \cdot 10^{3} \cdot \rho _{0,j}\Delta z_{j}/\Delta x ),$ (15)
$\displaystyle \mbox{DWNLV}_{i,j+\frac{1}{2}}^{n}$ $\textstyle =$ $\displaystyle \left\{
\left[
\left(w_{i+1,j+\frac{1}{2}}^{n}-w_{i,j+\frac{1}{2}...
...ft(w_{i,j+\frac{1}{2}}^{n}-w_{i-1,j+\frac{1}{2}}^{n}\right)^{3}
\right]
\right.$  
    $\displaystyle \left.
\left.
+ 0.1
\left[ \left(w_{i,j++\frac{3}{2}}^{n}-w_{i,j+...
...,j+\frac{1}{2}}^{n}-w_{i,j-\frac{1}{2}}^{n}\right)^{3}
\right]
\right\}
\right/$  
    $\displaystyle (16.0 \cdot 10^{3} \cdot \rho _{0,j+\frac{1}{2}}
\Delta z_{j+\frac{1}{2}}/\Delta x ),$ (16)
$\displaystyle \mbox{BUOY}_{i,j+\frac{1}{2}}^{n}$ $\textstyle =$ $\displaystyle \frac{g}{\Theta _{0,j+\frac{1}{2}}}
\theta _{i,j+\frac{1}{2}}^{n},$ (17)
$\displaystyle \mbox{DUCOLI}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle - f v_{i+\frac{1}{2},j}^{n},$ (18)
$\displaystyle \mbox{DVCOLI}_{i+\frac{1}{2},j}^{n}$ $\textstyle =$ $\displaystyle + f u_{i+\frac{1}{2},j}^{n}.$ (19)

\begin{eqnarray*}
&& u_{i+\frac{1}{2},j+\frac{1}{2}}^{n} =
0.5\left(u_{i+\frac...
...w_{i,j}^{n}}
{\rho _{0,j+\frac{1}{2}}\Delta z_{j+\frac{1}{2}}}.
\end{eqnarray*}

2.2 熱力学の式

移流項 $[\mbox{DTADV}]_{i,j}^{n}, [\mbox{DTAD0}]_{i,j}^{n}$ は 4 次中央 差分で空間離散化する. 摩擦項 $[\mbox{DTDIF}]_{i,j}^{N}$, $[\mbox{DTDI0}]_{i,j}^{N}$ と放射加熱項 $Q_{rad,i,j}^{N}$ と散逸加熱項 $Q_{dis,i,j}^{N}$ の時間積分は常に前進差分を用いて行う. 放射加熱項の計算方法は第5節を参照.

$\displaystyle \theta _{i,j}^{n+1}$ $\textstyle =$ $\displaystyle \theta _{i,j}^{N} +
dt \left\{
\frac{\Theta _{0,j}}{T_{0,j}}(Q_{rad,i,j}^{N} + Q_{dis,i,j}^{N})
\right.$  
    $\displaystyle + \left.
[\mbox{DTADV}]_{i,j}^{n} +
[\mbox{DTAD0}]_{i,j}^{n} +
[\mbox{DTDIF}]_{i,j}^{N} +
[\mbox{DTDI0}]_{i,j}^{N} \right\}$ (20)


\begin{displaymath}
Q_{dis,i,j}^{N} = \frac{C_{\epsilon}}{lc_{p}}
(\varepsilon _{i,j}^{N})^{\frac{3}{2}},
\end{displaymath} (21)


    $\displaystyle \mbox{DTADV}_{i,j}^{n}
= - \left\{
\frac{1}{\rho _{0,j}\Delta x}\...
...c{1}{2},j)}^{n}
+ \frac{1}{24}F\theta _{x(i-\frac{3}{2},j)}^{n} \right]
\right.$  
    $\displaystyle + \left.
\frac{1}{\rho _{0,j}\Delta z_{j}}\left[
- \frac{1}{24}F\...
...c{1}{2})}^{n}
+ \frac{1}{24}F\theta _{z(i,j-\frac{3}{2})}^{n} \right] \right\},$ (22)
    $\displaystyle F\theta _{x(i+\frac{1}{2},j)}^{n}
=
\rho _{0,j}u_{i+\frac{1}{2},j...
...^{n}
+ \frac{9}{16}\theta _{i,j}^{n}
- \frac{1}{16}\theta _{i-1,j}^{n} \right),$  
    $\displaystyle F\theta _{z(i,j+\frac{1}{2})}^{n}
=
\rho _{0,j+\frac{1}{2}}w_{i,j...
...^{n}
+ \frac{9}{16}\theta _{i,j}^{n}
- \frac{1}{16}\theta _{i,j-1}^{n} \right).$  


    $\displaystyle \mbox{DTAD0}_{i,j}^{n}
= - \left\{
\frac{1}{\rho _{0,j}\Delta x}\...
...(i-\frac{1}{2},j)}^{n}
+ \frac{1}{24}F_{x(i-\frac{3}{2},j)}^{n} \right]
\right.$  
    $\displaystyle + \left.
\frac{1}{\rho _{0,j}\Delta z_{j}}\left[
- \frac{1}{24}F\...
...c{1}{2})}^{n}
+ \frac{1}{24}F\Theta _{0,z(j-\frac{3}{2})}^{n} \right] \right\},$ (23)
    $\displaystyle F_{x(i+\frac{1}{2},j)}^{n}
=
\rho _{0,j}u_{i+\frac{1}{2},j}^{n},$  
    $\displaystyle F\Theta _{0,z(j+\frac{1}{2})}^{n}
=
\rho _{0,j+\frac{1}{2}}w_{i,j...
...eta _{0,j+1}
+ \frac{9}{16}\Theta _{0,j}
- \frac{1}{16}\Theta _{0,j-1} \right).$  


$\displaystyle \mbox{DTDIF}_{i,j}^{n}$ $\textstyle =$ $\displaystyle \frac{1}{(\Delta x)^{2}}
\left[
\tilde{K}_{i+\frac{1}{2},j}^{n}(\...
...ilde{K}_{i-\frac{1}{2},j}^{n}(\theta _{i,j}^{n} -\theta _{i-1,j}^{n})
\right] +$  
    $\displaystyle \frac{1}{\rho _{0,j}\Delta z_{j}}
\left[ \rho _{0,j+\frac{1}{2}}
...
...c{
(\theta _{i,j}^{n} -\theta _{i,j-1}^{n})}{\Delta z_{j-\frac{1}{2}}}
\right],$ (24)


\begin{displaymath}
\mbox{DTDI0}_{i,j}^{n}
=
\frac{1}{\rho _{0,j}\Delta z_{j}...
...{n} -\Theta _{0,j-1}^{n})}{\Delta z_{j-\frac{1}{2}}}
\right],
\end{displaymath} (25)


\begin{displaymath}
\tilde{K}_{i+\frac{1}{2},j}^{n}
= 0.5(\tilde{K}_{i+1,j}^{...
...{2}}^{n}
= 0.5(\tilde{K}_{i,j+1}^{n} +\tilde{K}_{i,j}^{n} ).
\end{displaymath}

2.3 圧力診断式

圧力の診断式([*])を第2.1節で行った変形 にあわせて以下のように変形する.

\begin{displaymath}
\left(\DP[2]{}{x} +\frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{...
...\rho _{0}}\Ddiv \rho _{0}\Dvect{v}\right)
+ \DP{}{x}\alpha ,
\end{displaymath}

これを dimension-reduction 法を用いて解く. 上の式を適当に離散化すると,
\begin{displaymath}
\Dvect{D_{x}P} + \Dvect{D_{z}P} = \Dvect{S}
\end{displaymath} (26)

と書くことができる. ただし $\Dvect{P}, \Dvect{D_{x}}, \Dvect{D_{z}},
\Dvect{S}$ はそれぞれ,

\begin{displaymath}
\hat{P}, \quad
\DP[2]{}{x}, \quad
\frac{1}{\rho _{0}}\DP{...
...frac{1}{\rho _{0}}\Ddiv \Dvect{v}\right)
+ \DP{}{x}\alpha ,
\end{displaymath}

を離散化した行列である. 行列 $\Dvect{D_{z}}$ の固有値を $\lambda_{i}$, 固有ベクトルを $V_{i}(z)$ とし, 固有列ベクトルを並べた行列を $\Dvect{V}$, 対角成分が $\lambda_{i}$ である行列を $\Dvect{\Lambda }$ とすると, $\Dvect{D_{z}V} = \Dvect{V\Lambda }$ となる. $\Dvect{P} = \Dvect{V}\cdot
\Dvect{H}$ と展開すると,

\begin{displaymath}
\Dvect{V}\Dvect{D_{x}H} + \Dvect{V\Lambda H} = \Dvect{S}.
\end{displaymath}

したがって,
\begin{displaymath}
\left(\Dvect{D_{x}} + \Dvect{\Lambda }\right) \Dvect{H}
= \Dvect{V}^{-1}\Dvect{S},
\end{displaymath} (27)

となる.

固有値を $\lambda_{i}$ と固有ベクトルを $V_{i}(z)$ を求めるために必要な $z$ 方向の係数行列 $\Dvect{D}_{z}$

\begin{displaymath}
\frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{}{z} \pi
\end{displaymath}

を差分化して与える. 最初の微分は 2 次中央差分, 次の微分は 4 次中央差分で 解く. これは連続の式は 4 次中央差分, 圧力の微分は 2 次中央差分であること による. よって係数行列は 5重対角行列になる. 係数行列の成分を $A_{ij}$ と すると以下のようになる.
$\displaystyle A_{i,i+2}$ $\textstyle =$ $\displaystyle -\frac{1}{\rho _{0}\Delta z_{i}}\left(
\frac{1}{24\Delta z_{i+\frac{3}{2}}}\frac{(\rho _
{0,i+2}+\rho _{0,i+1})}{2}\right),$ (28)
$\displaystyle A_{i,i+1}$ $\textstyle =$ $\displaystyle \frac{1}{\rho _{0}\Delta z_{i}}\left(
\frac{1}{24\Delta z_{i+\fra...
...rac{9}{8\Delta m_{i+\frac{1}{2}}}\frac{(\rho _
{0,i+1}+\rho _{0,i})}{2}\right),$ (29)
$\displaystyle A_{i,i}$ $\textstyle =$ $\displaystyle - \frac{1}{\rho _{0}\Delta z_{i}}\left(
\frac{9}{8\Delta zm_{i+\f...
...rac{9}{8\Delta z_{i-\frac{1}{2}}}\frac{(\rho _
{0,i}+\rho _{0,i-1})}{2}\right),$ (30)
$\displaystyle A_{i+1,i}$ $\textstyle =$ $\displaystyle \frac{1}{\rho _{0}\Delta z_{i}}\left(
\frac{1}{24\Delta z_{i-\fra...
...rac{9}{8\Delta z_{i-\frac{1}{2}}}\frac{(\rho _
{0,i}+\rho _{0,i-1})}{2}\right),$ (31)
$\displaystyle A_{i+2,i}$ $\textstyle =$ $\displaystyle -\frac{1}{\rho _{0}\Delta z_{i}}\left(
\frac{1}{24\Delta z_{i-\frac{3}{2}}}\frac{(\rho _
{0,i-1}+\rho _{0,i-2})}{2} \right).$ (32)

境界条件は上下壁で $\DP{}{z}=0$ となるようにする.

$x$ 方向についても適当な固有関数に展開して各モードに対する展開係数 を求める. ここでは三角関数で展開する.

$\displaystyle \Dvect{H}$ $\textstyle =$ $\displaystyle \sum_{k=1}^{NX/2-1}
\left[\Dvect{H}\right]_{k_{x}},$ (33)
$\displaystyle \Dvect{V}^{-1}\Dvect{S}$ $\textstyle =$ $\displaystyle \sum_{k=1}^{NX/2-1}
\left[\Dvect{V}^{-1}\Dvect{S}\right]_{k_{x}},$ (34)
$\displaystyle \left(\Dvect{D_{x}} + \Dvect{\Lambda }\right)$ $\textstyle =$ $\displaystyle \sum_{k_{x}=1}^{NX/2-1}
\left[\Dvect{D_{x}} + \Dvect{\Lambda }\right]_{k_{x}}$ (35)


\begin{displaymath}
\left[\Dvect{H}\right]_{k_{x}} =
\left[\Dvect{D_{x}} + \D...
...ht]_{k_{x}}^{-1}
\left[\Dvect{V}^{-1}\Dvect{S}\right]_{k_{x}}
\end{displaymath} (36)

2.4 基本場

温度分布 $T_{j}$ を与え, 静水圧平衡の式と状態方程式を用いて $P_{0,j}$$\rho _{0,j}$ を計算する.

    $\displaystyle \ln P_{0,j} = \ln P_{00}
- \sum _{j=1}^{j}\frac{g}{RT_{0,j}}\Delta z_{j},$ (37)
    $\displaystyle \rho _{0,j} = \frac{P_{0,j}}{RT_{0,j}}.$ (38)

$P_{0,j}, \rho _{0,j}$ を求めた後 $\Pi _{0,j}, \Theta _{0,j}$ を計算する.
    $\displaystyle \Pi _{0,j} = \left(\frac{P_{0,j}}{P_00}\right)^{\kappa },$ (39)
    $\displaystyle \Theta _{0,j} = \frac{T_{0,j}}{\Pi _{0,j}}.$ (40)


next up previous
: 3 乱流モデル : Two dimensional anelastic model : 1 モデル離散化の概要
Odaka Masatsugu 平成19年4月26日