%
% An Introduction to Atmospheric Gravity Waves (2002) 
% 1st edition
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Style  Setting             %%%%%%%%
% フォント: 12point (最大), 片面印刷
\documentclass[a4j,12pt,oneside]{ltjreport}

\input{preemble}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Text Start                 %%%%%%%%
\begin{document}
%%\appendix
%\setcounter{chapter}{0}
%\setcounter{section}{0}
%\setcounter{equation}{0}
%%\def\thesection{付録\Alph{section}}
%%\def\theequation{\Alph{section}.\arabic{equation}}
\chapter{モデルについて}
\section{基礎方程式系}

対流圏における対流の発生, その対流からの内部重力波の発生, 
クリティカルレベルにおける内部重力波の吸収,
クリティカルレベルからの二次的な内部重力波の発生, 
その波の高高度への伝播と砕波をシミュレーションするため, 
非線形・非静力学・圧縮性の鉛直2次元モデルを作成した. 
シミュレーションを行う内部重力波の振動数は
局所的な慣性振動数よりも大きいことから,
地球の自転の効果は無視する.
基礎方程式系は\cite{goya1999non}
で用いられたものを使用した.
基礎方程式系は次のとおりである:
\begin{subequations}\label{main_equations.eq}
\begin{align}
\pdv{\delta\tilde{u}}{t}
=&-\pdv{\tilde{u}u}{x}-\pdv{\tilde{u}w}{z}
-\pdv{p}{x}
-K_H\pdv[4]{\delta\tilde{u}}{x}
+K_V\pdv[2]{\delta\tilde{u}}{z}
-R_f\delta\tilde{u},
\label{1a.eq}
\\
\pdv{\delta\tilde{w}}{t}
=&-\pdv{\tilde{w}u}{x}-\pdv{\tilde{w}w}{z}
-\pdv{p}{z}-\rho g
-K_H\pdv[4]{\delta\tilde{w}}{x}
+K_V\pdv[2]{\delta\tilde{w}}{z}
-R_f\delta\tilde{w},
\label{1b.eq}
\\
\pdv{\delta\tilde{T}}{t}
=&-\pdv{\tilde{u}T}{x}-\pdv{\tilde{w}T}{z}
-\frac{R\tilde{T}}{c_v}
\left(\pdv{u}{x}+\pdv{w}{z}\right)
-K_H\pdv[4]{\delta\tilde{T}}{x}
+K_V\pdv[2]{\delta\tilde{T}}{z}
-\alpha\delta\tilde{T},
\label{1c.eq}
\\
\pdv{\rho}{t}+&\pdv{\tilde{u}}{x}+\pdv{\tilde{w}}{z}=0,
\label{1d.eq}
\\
p=&R\tilde{T}.
\label{1e.eq}
\end{align}	
\end{subequations}
記号の定義は以下のとおりである:
\begin{align*}
u &:\text{水平速度},\\
w &:\text{鉛直速度},\\
p &:\text{圧力},\\
\rho &:\text{密度},\\
T &:\text{温度},\\
g &:\text{重力加速度},\\
R &:\text{乾燥大気の気体定数},\\
c_v &:\text{定積比熱},\\
K_H &:\text{水平方向の渦拡散係数},\\
K_V &:\text{鉛直方向の渦拡散係数},\\
R_f &:\text{レイリー摩擦係数},\\
\alpha &:\text{ニュートン冷却係数},\\
\end{align*}
力学変数[$ u,w,T $]と$ \rho $の積を $ \tilde{} $ つきの量で表す:
例えば
\begin{equation}
\tilde{u}\equiv u(x,z,t)\cdot\rho(x,z,t)
\end{equation}
である.
さらに初期の基本場からの変化を表す場合は$ \delta $つきの量で表す.
すなわち
\begin{subequations}
\begin{align}
&\delta\tilde{u}=\tilde{u}(x,z,t)-\tilde{u}_B(z),
\;\;\;\;\;
\tilde{u}_B\equiv\rho_B(z)u_B(z),
\label{tilde_u.eq}
\\
&\delta\tilde{w}(x,z,t)=\tilde{w}(x,z,t),
\label{tilde_w.eq}
\\
&\delta\tilde{T}=\tilde{T}(x,z,t)-\tilde{T}_B(z),
\;\;\;\;\;
\tilde{T}_B\equiv\rho_B(z)T_B(z).
\label{tilde_T.eq}
\end{align}
\end{subequations}
添字の$ B $は基本場の量を表す.


\section{モデルの設定と計算方法}

\begin{figure}[htb]
	\begin{center}
		\includegraphics[width=10cm]{./mthesis_fig/grid.jpg}
		\begin{quote}
			\begin{footnotesize}
				\caption{\label{grid.fig}
				格子点の配置.
				空間変数を$ (x,z)=(i\Delta x,j\Delta z ) $
				と差分化した時, 格子点$ (i,j) $に
				熱力学変数を配置し, そこから$ 1/2 $格子点だけ
				$ x,z $方向にずれた点でそれぞれ$ u,w $
				を配置する.}
			\end{footnotesize}
		\end{quote}
	\end{center}
\end{figure}

\eqref{main_equations.eq} の数値計算には図\ref{grid.fig} のように
物理変数を配置したスタッガード格子を用いる.
%\chapter{Goya \& Miyahara (1999)の方程式の差分化}
%\setcounter{equation}{0}
差分化した基礎方程式系は以下のとおりである.
\begin{subequations}

\subsection*{運動量方程式}

\begin{align}
\qty(\pdv{\delta\tilde{u}}{t})_{i+\frac{1}{2},j}
=
&-\frac{\tilde{u}_{i+1,j}u_{i+1,j}
	-\tilde{u}_{i,j}u_{i,j}}
{\Delta x}
\notag\\
&-\frac{\tilde{u}_{i+\frac{1}{2},j+\frac{1}{2}}
	w_{i+\frac{1}{2},j+\frac{1}{2}}
	-\tilde{u}_{i+\frac{1}{2},j-\frac{1}{2}}
	w_{i+\frac{1}{2},j-\frac{1}{2}}}
{\Delta z}
\notag\\
&-\frac{p_{i+1,j}-p_{i,j}}{\Delta x}
\notag\\
&-K_H\frac{
	\delta\tilde{u}_{i+\frac{5}{2},j}
	-4\delta\tilde{u}_{i+\frac{3}{2},j}
	+6\delta\tilde{u}_{i+\frac{1}{2},j}
	-4\delta\tilde{u}_{i-\frac{1}{2},j}
	+\delta\tilde{u}_{i-\frac{3}{2},j}	
}
{(\Delta x)^4}
\notag\\
&+K_V\frac{
	\delta\tilde{u}_{i+\frac{1}{2},j+1}
	-2\delta\tilde{u}_{i+\frac{1}{2},j}
	+\delta\tilde{u}_{i+\frac{1}{2},j-1}
}
{(\Delta z)^2}
\notag\\
&-R_{fj}\delta\tilde{u}_{i+\frac{1}{2},j}
\end{align}

\begin{align}
\qty(\pdv{\delta\tilde{w}}{t})_{i,j+\frac{1}{2}}
=
&-\frac{\tilde{w}_{i+\frac{1}{2},j+\frac{1}{2}}
	u_{i+\frac{1}{2},j+\frac{1}{2}}
	-\tilde{w}_{i-\frac{1}{2},j+\frac{1}{2}}
	u_{i-\frac{1}{2},j+\frac{1}{2}}}
{\Delta x}
\notag\\
&-\frac{\tilde{w}_{i+1,j}w_{i+1,j}
	-\tilde{w}_{i,j}w_{i,j}}
{\Delta z}
\notag\\
&-\frac{\delta p_{i,j+1}-\delta p_{i,j}}{\Delta z}
-g\frac{\delta\rho_{i,j+1}+\delta\rho_{i,j}}{2}
\notag\\
&-K_H\frac{
	\delta\tilde{w}_{i+2,j+\frac{1}{2}}
	-4\delta\tilde{w}_{i+1,j+\frac{1}{2}}
	+6\delta\tilde{u}_{i,j+\frac{1}{2}}
	-4\delta\tilde{u}_{i-1,j+\frac{1}{2}}
	+\delta\tilde{u}_{i-2,j+\frac{1}{2}}	
}
{(\Delta x)^4}
\notag\\
&+K_V\frac{
	\delta\tilde{w}_{i,j++\frac{3}{2}}
	-2\delta\tilde{u}_{i,j+\frac{1}{2}}
	+\delta\tilde{u}_{i,j-\frac{1}{2}}
}
{(\Delta z)^2}
\notag\\
&-R_{fj}\delta\tilde{w}_{i,j+\frac{1}{2}}
\end{align}
ここで基本場について静水圧平衡が成り立つことから
\begin{equation*}
\frac{(p_B)_{i,j+1}-(p_B)_{i,j}}{\Delta z}
+g\frac{(\rho_B)_{i,j+1}+(\rho_B)_{i,j}}{2}
=0
\end{equation*}
であることを用いた.

\subsection*{連続の式}

\begin{equation}\label{ContEq_d}
\qty(\pdv{\rho}{t})_{i,j}
=-\frac{\tilde{u}_{i+\frac{1}{2},j}
	-\tilde{u}_{i-\frac{1}{2},j}}
{\Delta x}
-\frac{\tilde{w}_{i,j+\frac{1}{2}}
	-\tilde{u}_{i,j-\frac{1}{2}}}
{\Delta z}
\end{equation}

\subsection*{熱力学の式}
\begin{align}
\qty(\pdv{\delta\tilde{T}}{t})_{i,j}
=
&-\frac{
	\tilde{u}_{i+\frac{1}{2},j}T_{i+\frac{1}{2},j}
	-\tilde{u}_{i-\frac{1}{2},j}T_{i-\frac{1}{2},j}
}{\Delta x}
\notag \\
&-\frac{
	\tilde{w}_{i,j+\frac{1}{2}}T_{i,j+\frac{1}{2}}
	-\tilde{w}_{i,j-\frac{1}{2}}T_{i,j-\frac{1}{2}}
}{\Delta z}
\notag \\
&+\frac{R\tilde{T}_{I,j}}{c_v}
\qty[
\frac{u_{i+\frac{1}{2},j}-u_{i-\frac{1}{2},j}}{\Delta x}
+\frac{w_{i,j+\frac{1}{2}}-w_{i,j-\frac{1}{2}}}{\Delta z}
]
\notag \\
&-K_H\frac{
	\delta\tilde{T}_{i+2,j}
	-4\delta\tilde{T}_{i+1,j}
	+6\delta\tilde{T}_{i,j}
	-4\delta\tilde{T}_{i-1,j}
	+\delta\tilde{T}_{i-2,j}}
{(\Delta x)^4}
\notag \\
&+K_V\frac{
	\delta\tilde{T}_{i,j+1}
	-2\delta\tilde{T}_{i,j}
	+\delta\tilde{T}_{i,j-1}}
{(\Delta z)^2}
\notag \\
&-\alpha_{i,j}\delta\tilde{T}_{i,j}
\end{align}
\end{subequations}

モデルの領域は水平方向に$ \SI{30}{km} $,
鉛直方向に$ \SI{150}{km} $である.
時間積分の計算には leap-frog 法と Asselin フィルタを用いた.


\section{基本状態}

\begin{figure}[htbp]
\begin{tabular}{ll}
	\begin{minipage}{0.49\hsize}
		\begin{center}
		\includegraphics[height=8cm]{./mthesis_fig/U_b.png}
		\begin{quote} 
			\begin{footnotesize} 
				\caption{\label{U_b.fig}
					水平風の基本場の鉛直プロファイル
					$ U_B(z)$ $[\si{m.s^{-1}]} $.} 
			\end{footnotesize}
		\end{quote} 
		\end{center}
	\end{minipage}
% % % % % % % % % % % % % % % % % % % % % %
	\begin{minipage}{0.49\hsize}
		\begin{center}
		\includegraphics[height=8.3cm]{./mthesis_fig/T_b.png}
		\begin{quote}
		\begin{footnotesize}
			\caption{
				温度の基本場の鉛直プロファイル
				$ T_B(z)\ [\si{K}] $.
				\label{T_b.fig}} 
		\end{footnotesize}		
		\end{quote}
		\end{center}
	\end{minipage}
	\\
	\begin{minipage}{0.49\hsize}
		\begin{center}
		\includegraphics[height=8cm]{./mthesis_fig/K_v.png}
		\begin{quote}
		\begin{footnotesize}
			\caption{\label{K_v.fig}
				鉛直渦粘性係数の鉛直プロファイル
				$ K_V(z)$ $[\si{m^2.s^{-1}}] $.}
		\end{footnotesize}
		\end{quote}
		\end{center}	
	\end{minipage}
	\begin{minipage}{0.49\hsize}
		\begin{center}
		\includegraphics[height=8cm]{./mthesis_fig/Rf.png}
		\begin{quote}
		\begin{footnotesize}
		\caption{\label{Rf.fig}
			レイリー摩擦係数の鉛直プロファイル
			$ R_f(z)\ [\si{s^{-1}}] $.}
		\end{footnotesize}
		\end{quote}
		\end{center}
	\end{minipage}
\end{tabular}
\end{figure}
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 二列目は仮置き
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

図\ref{U_b.fig}と図\ref{T_b.fig}にそれぞれ基本状態の
$ U_B(z) $と$ T_B(z) $を示す.
$ U_B(z) $については, すべての高度でゼロの場合と,
中緯度における対流圏と下部成層圏の東西風の鉛直プロファイルを想定した
場合の2つの計算を行った.
ただし, $ \SI{24}{km} $より上は内部重力波の伝播を見るために
風速をゼロとした.
$ T_B(z) $は現実的な大気の鉛直プロファイルを用いた.
ただし, 高度$ \SI{5}{km} $より下では対流を発生させるために
温度減率を乾燥断熱減率よりも僅かに大きい
$ \dv*{T_B(z)}{z}=\SI{-10}{K.km^{-1}} $とした.
また, 数値シミュレーションが困難な
熱圏での極端に速い音速を避けるため, $ \SI{110}{km} $よりも
上では$ \SI{300}{K} $で一定と仮定する.

\section{散逸}

水平渦粘性係数は$ K_H=1.0\times10^{7}\ \si{m^4.s^{-1}} $
で定数とする.
図\ref{K_v.fig}と図\ref{Rf.fig}はそれぞれ
鉛直渦粘性係数とレイリー摩擦係数を表す.
ニュートン冷却係数$ \alpha(z) $は$ R_f(z) $と同じであると仮定する.
上端からの擬似的な波の反射を防ぐために$ K_V $と$ R_f $には
高さとともに増加する値を与える.
鉛直渦拡散の項は, 例えば\eqref{1a.eq}ならば
物理的に正しい
$ \pdv{z}(\rho K_V \pdv{u}{z}) $
の形式を採用すべきであるが,
簡単のため$ K_V\pdv*[2]{\delta\tilde{u}}{z} $の形式を用いる.
\cite{goya1999non}によれば
シミュレーションを行う範囲においては2つのスキームの
差が与える影響は小さい.
\section{モデルの解像度}

格子点数は
水平方向に128点($ \Delta x\simeq \SI{234}{m} $),
鉛直方向に500点($ \Delta z = \SI{300}{m} $)とした.
$ \Delta t =\SI{0.1}{s} $, フィルタの係数は$ 0.2 $とした.

\section{初期条件と境界条件}

初期条件は$ \delta \tilde{w} $の初期摂動を除いてゼロとする:
\begin{equation}
\delta\tilde{u}=0,\;\;\delta\tilde{w}=0,\;\;
\delta\tilde{T}=0,\;\;\delta \rho=0,\;\;
\delta p=0.
\end{equation}
最下層の$ x=\SI{15}{km} $の一点において$ \delta\tilde{w} $に
$ 4\times 10^{-4}\ \si{m^{-2}.kg.s^{-1}} $の
振幅を持つ摂動を与える.

また境界条件について, 上端の境界条件をゼロとする:
\begin{equation}
\delta\tilde{u}^N=0,\;\;\delta\tilde{w}^N=0,\;\;
\delta\tilde{T}^N=0.
\end{equation}
下端の境界条件は以下のように設定する.
すなわち$ \delta\tilde{u} $の鉛直微分をゼロ, 
$ \delta\tilde{w} $と$ \delta\tilde{T} $をゼロとする:
\begin{equation}
\delta\tilde{u}^0=\delta\tilde{u}^1,\;\;
\delta\tilde{w}^0=0,\;\;
\delta\tilde{T}^0=0.
\end{equation}
水平方向の境界条件については周期境界とする.


\end{document}


