%
% 表題   SPMODEL サンプルプログラムドキュメント
%        kh.tex
%
% 履歴   2004/06/01 小高 正嗣  
%        2004/06/08 小高 正嗣  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4j,12pt,dvipdfmx]{jarticle}
%
% プリアンブル
%
\usepackage{Dennou6}                      % dennou-sty-6 を使用

\Dtitle[{\small 2 次元ブシネスクモデル: KH 不安定}] % ヘッダ部タイトル, モデル名と同じ
{ {\large SPMODEL サンプルプログラム} \\  % ここは変更しない 
  水路領域 \\ 2 次元ブシネスク方程式モデル:\\           % モデル名
  ケルビン・ヘルムホルツ不安定の計算 \\
  {\large kh1.f90}                   % プログラム名
}
%
%
\Dauthor[小高 正嗣, 竹広 真一]{小高 正嗣, 竹広 真一}            % 著者
% \Ddate{2004 年 6 月 8 日}                % 日付
\Ddate{2026 年 7 月 8 日}                % 日付
\Dnoparindent                             % 段落の字下げをしない
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 本文開始
%
% 基本的な章節構成は以下の通り. 
%
%   1. はじめに 
%   2. 支配方程式系
%   3. 離散化
%   4. 利用モジュールとその他の設定
%   5. 数値実験
%   6. 参考文献
%   -. 謝辞
% 
% 必要に応じて構成を変更すること.
%
\begin{document}

\maketitle                                % タイトルの作成
\tableofcontents                          % 目次の作成
\pagebreak                                

%----------------------------------------------------------------------
\Dparskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section{概要}

SPMODEL サンプルプログラム『kh1.f90』に用いられている基礎方程式と境界条
件, および, このプログラムを用いた数値実験の方法について解説する. 基礎方
程式は 2 次元のブシネスク方程式系である. 計算はスペクトル法を用いて行い, 
展開関数は水平方向にはフーリエ級数, 鉛直方法にはチェビシェフ多項式を用いる.
スペクトル変換と逆変換および微分演算には, SPMODEL ライブラリ(spml) の 
et\_module を用いている. 数値実験では水路中心で向きが反転する平行流を与
えた場合に生じるケルビン・ヘルムホルツ不安定の計算を行う.

{\bf プログラム名} \\                  % プログラム名
{\footnotesize
kh1.f90
}

{\bf プログラム取得元}\\               % プログラム取得先
{\footnotesize
http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 2d-channel-et\slash boussinesq\slash kh-instability\slash SIGEN.htm
}


{\bf SPMODEL サンプルプログラム目次}\\ % サンプルプログラム目次(変更しない)
{\footnotesize
 http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash sample.htm 
}

{\bf SPMODEL の使い方}\\               % SPMODEL の使い方(変更しない)
{\footnotesize
http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\/
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section{支配方程式系}

ここでは支配方程式系と境界条件を記す.

\subsection{支配方程式系}

支配方程式系は 2 次元のブシネスク方程式系である(Chandrasekhar, 1961).
\begin{eqnarray}
&& \DP{u}{t}+(U+u)\DP{u}{x}+v\DP{(U+u)}{y} = -\frac{1}{\rho _{0}}\DP{p}{x} 
                                     + \nu \Dgrad ^{2p} u, \\
&& \DP{v}{t}+(U+u)\DP{v}{x}+v\DP{v}{y} = -\frac{1}{\rho _{0}}\DP{p}{y} 
                                     - \frac{\rho}{\rho _{0}}g
                                     + \nu \Dgrad ^{2p} v, \\
&& \DP{u}{x}+\DP{v}{y} = 0,\\
&& \DP{\rho}{t}+(U+u)\DP{\rho}{x}+v\DP{\rho}{y} = \kappa \Dgrad^{2p}\rho.
\end{eqnarray}
各記号の定義は\Dtabref{変数, 物理定数の定義}に表す. 基本場の速度 $U$ は
$y$ にのみ依存すると仮定する.

渦度 $\zeta$ と流線関数 $\psi$ 
\begin{eqnarray}
  &&\zeta =\DP{v}{x}-\DP{u}{y}, \\
  &&v=\DP{\psi}{x}, \quad u=-\DP{\psi}{y},
\end{eqnarray}
を導入し, 渦度と密度の予報方程式として支配方程式系を表すと以下のようになる
\footnote{渦度方程式に表れる $-\DD{U}{z}\DP{u}{x}$ 項は無視する.}
\begin{eqnarray}
&& \DP{\zeta }{t} 
   + U \DP{\zeta}{x} - \DD[2]{U}{y}\DP{\psi}{x}
   + J(\psi, \zeta) + \frac{g}{\rho _{0}}\DP{\rho }{x}
   = \nu \Dgrad ^{2p} \zeta + F_{\zeta}, \Deqlab{渦度方程式}\\
&& \DP{\rho}{t} + 
   U\DP{\rho }{x} + 
   J(\psi, \rho) 
   = \kappa \Dgrad ^{2p} \rho. \Deqlab{密度の予報方程式}
\end{eqnarray}
ここで $J(A, B)$ はヤコビアン
\[
  J(A,B) = \DP{A}{x}\DP{B}{y} - \DP{A}{y}\DP{B}{x}
\]
である.

\begin{table}[h]
 \begin{center}
  \begin{tabular}{c l l  } \hline
   記号    &\qquad\qquad&  変数/物理定数   \\ \hline 
   $x$      && 水平座標        \\
   $y$      && 鉛直座標        \\
   $t$      && 時間            \\ 
   $u$      && $x$ 方向速度    \\ 
   $U(y)$   && $x$ 方向速度(基本場)\\ 
   $v$      && $y$ 方向速度    \\ 
   $\rho_{0}$ && 平均密度      \\
   $\rho $   && 密度           \\
   $\psi$   && 流線関数        \\ 
   $\zeta$  && 渦度            \\ 
   $g$      && 重力加速度      \\ 
   $p     $ && 超粘性の次数    \\
   $\nu$    && 動粘性係数      \\ 
   $\kappa$ && 熱拡散係数      \\ \hline
 \end{tabular} 
\caption{変数, 物理定数の定義}
\Dtablab{変数, 物理定数の定義}
\end{center}
\end{table}

\subsection{境界条件}

境界条件は水平に周期境界条件, 鉛直方向には $y=0, y_{m}$ に置いた剛体壁
面で $v=0$, 応力なし, 熱フラックスなしとする. 
水平方向の境界条件は計算領域をそれぞれ $x_{m}$ とすると
\begin{equation}
  \zeta(x+x_{m},y) = \zeta(x,y)
\end{equation}
などと表される. 鉛直方向の境界条件は $y=0, y_{m}$ において,
\begin{eqnarray}
 &&\psi = \mbox{Const.}, \\
 &&\zeta = 0, \\
 &&\DP{\rho}{y} = 0
\end{eqnarray}
である. 境界で与える $\psi $ の値は簡単のため $0$ とした.


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \newpage                                  % 改ページ
% \section{離散化}
% この節では方程式系の空間離散化および使用した時間積分法について説明し, プロ
% グラム内で実際に用いられている方程式を記述する. 

% \subsection{空間離散化}
% 支配方程式\Deqref{渦度方程式}, \Deqref{密度の予報方程式}の離散表現は以下
% のようになる.
% \begin{eqnarray}
% && \DP{\zeta _{i,j}}{t} + 
%    U_{i,j}\DP{\zeta _{i,j}}{x} -
%    \DD[2]{U_{i,j}}{z}\DP{\psi_{i,j}}{x} + 
%    \DP{\psi_{i,j}}{x}\DP{\zeta _{i,j}}{y} - 
%    \DP{\psi_{i,j}}{y}\DP{\zeta _{i,j}}{x}  
%    + \frac{g}{\rho _{0}}\DP{\rho _{i,j}}{x} \\ \nonumber
% && \quad  
%    = \nu \Dgrad ^{2p} \zeta_{i,j} , \Deqlab{離散化された渦度方程式}\\
% && \DP{\rho_{i,j}}{t} + 
%    U_{i,j}\DP{\rho _{i,j}}{x} + 
%    \DP{\psi_{i,j}}{x}\DP{\rho _{i,j}}{y} - 
%    \DP{\psi_{i,j}}{y}\DP{\rho _{i,j}}{x}  
%    = \kappa \Dgrad ^{2p} \rho_{i,j}. \Deqlab{離散化された密度の予報方程式}
% \end{eqnarray}
% ここで添字 $i,j$ は格子点 $(x_{i}, y_{j})$ 上の値であることを表す.


% \subsection{空間方向のスペクトル表現}
% 空間離散化した支配方程式\Deqref{離散化された渦度方程式}, \Deqref{離散
%   化された密度の予報方程式}をスペクトル法を用いて表現する.  スペクトル
% 展開は水平方向にフーリエ級数, 鉛直方法には境界条件にあわせてフーリエ正
% 弦級数またはフーリエ余弦級数を用いて行う. 非線形項を扱う場合は, 先に格
% 子点上での非線形項の値を計算し, その値のスペクトルを求める方法(変換法)
% を用いる. 浮力項についても同様に扱う. 以下では $k,l$  をそれぞれ$x,y$ 
% 方向波数, $K,L$ を切断波数, $I,J$ を格子点数とする.

% $\zeta_{i,j}, \psi _{i,j}, \rho _{i,j}$ はスペクトル逆変換によって以下
% のように展開される.
% \begin{eqnarray}
%  \zeta_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                  \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right)
%                  \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \zeta_{k,l}, \\
%  \psi_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                 \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right)
%                 \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \psi_{k,l}, \\
%  \rho_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                 \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right)
%                 \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)\hat \rho_{k,l}.
% \end{eqnarray}
% スペクトル係数 $\hat \zeta_{i,j}, \hat \psi _{i,j}, \hat \rho _{i,j}$ は
% スペクトル変換によって以下のように与えられる.
% \begin{eqnarray}
%  \hat \zeta_{k,l} &=& \frac{1}{I}\frac{1}{J}
%                       \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%                       \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%                       \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\zeta_{i,j}, \\
%  \hat \psi_{k,l} &=& \frac{1}{I}\frac{1}{J}
%                      \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%                      \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%                        \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)\psi_{i,j}, \\
%  \hat \rho_{k,l} &=& \frac{1}{I}\frac{1}{J}
%                      \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%                      \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%                      \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)\rho_{i,j}.
% \end{eqnarray}

% 以上を用いると支配方程式\Deqref{離散化された渦度方程式}, \Deqref{離散
% 化された密度の予報方程式}のスペクトル表現は以下のようになる.
% \begin{eqnarray}
%  \DP{\hat \zeta _{k,l}(t)}{t} &=& -
%       \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         U_{i,j}
%         \left(\DP{\zeta }{x}\right)_{i,j}  \nonumber \\
% && +  \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DD[2]{U}{y}\right)_{i,j}
%         \left(\DP{\zeta }{x}\right)_{i,j}  \nonumber \\
% && -  \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi}{x}\right)_{i,j}
%         \left(\DP{\zeta}{y}\right)_{i,j}  \nonumber \\
% &&  + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi}{y}\right)_{i,j}
%         \left(\DP{\zeta}{x}\right)_{i,j}  \nonumber \\
% && - \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \frac{g}{\rho _{0}}\left(\DP{\rho}{x}\right)_{i,j} \nonumber \\
% && - \nu \left[
%      \left(\frac{2\pi k}{x_{m}}\right)^{2p}+
%      \left(\frac{\pi l}{y_{m}}\right)^{2p}\right]\hat \zeta _{k,l},
%      \Deqlab{渦度方程式のスペクトル表現}\\
%  \DP{\hat \rho _{k,l}(t)}{t} &=& -
%       \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         U_{i,j}
%         \left(\DP{\rho}{x}\right)_{i,j}  \nonumber \\
% && -      \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi}{x}\right)_{i,j}
%         \left(\DP{\rho}{y}\right)_{i,j}  \nonumber \\
% &&  + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi}{y}\right)_{i,j}
%         \left(\DP{\rho}{x}\right)_{i,j}  \nonumber \\
% && - \kappa \left[
%      \left(\frac{2\pi k}{x_{m}}\right)^{2p}+
%      \left(\frac{\pi l}{y_{m}}\right)^{2p}\right]\hat \rho _{k,l}.
%      \Deqlab{密度の予報方程式のスペクトル表現}
% \end{eqnarray}
% ここで
% \begin{eqnarray}
%  \left(\DP{\zeta}{x}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \sin\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{2\pi i k}{x_{m}}\hat \zeta _{k,l}, \\
%  \left(\DP{\zeta}{y}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \cos\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{\pi l}{y_{m}}\hat \zeta _{k,l}, \\
%  \left(\DP{\psi}{x}\right)_{i,j} &=&  \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \sin\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{2\pi i k}{x_{m}}\hat \psi _{k,l}, \\
%  \left(\DP{\psi}{y}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \cos\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{\pi l}{y_{m}}\hat \psi _{k,l}, \\
%  \left(\DP{\rho}{x}\right)_{i,j} &=& \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \cos\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{2\pi i k}{x_{m}}\hat \rho _{k,l}, \\
%  \left(\DP{\rho}{y}\right)_{i,j} &=& - \sum _{k=-K}^{K}\sum _{l=0}^{L}
%                             \exp\left(\frac{2\pi i k }{x_{m}}\right)
%                             \sin\left(\frac{\pi l }{y_{m}}\right)
%                             \frac{\pi l}{y_{m}}\hat \rho _{k,l}
% \end{eqnarray}
% である.

% \subsection{時間積分}
% ここでは時間積分法について記述し, プログラム内で実際に用いられている方
% 程式を記述する. 以下では $\Delta t$ を時間格子間隔, 時刻 $\tau \Delta t$ 
% における $\hat \zeta_{k,l}$ の値を $\hat \zeta_{k,l}^{\tau}$ 等と表す.

% 時間方向の離散化は Euler スキームを用いて行う. 時空間方向に離散化され
% た方程式は以下のように表される.
% \begin{eqnarray}
%  &&\hat \zeta _{k,l}^{\tau +1} = \hat \zeta _{k,l}^{\tau} + \Delta t \hat F_{k,l}^{\tau},  \\
%  &&\hat \rho _{k,l}^{\tau +1} = \hat \rho _{k,l}^{\tau} + \Delta t \hat G_{k,l}^{\tau}, 
% \end{eqnarray}
% \begin{eqnarray}
% \hat F_{k,l}^{\tau} &=& -
%       \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         U_{i,j}
%         \left(\DP{\zeta }{x}\right)_{i,j}  \nonumber \\
% && +  \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DD[2]{U}{y}\right)_{i,j}
%         \left(\DP{\zeta }{x}\right)_{i,j}^{\tau}  \nonumber \\
% && -  \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi_{i,j}}{x}\right)_{i,j}^{\tau}
%         \left(\DP{\zeta _{i,j}}{y}\right)_{i,j}^{\tau}  \nonumber \\
% &&  + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi_{i,j}}{y}\right)_{i,j}^{\tau}
%         \left(\DP{\zeta _{i,j}}{x}\right)_{i,j}^{\tau}  \nonumber \\
% && - \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \sin\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \frac{g}{\rho _{0}}\left(\DP{\rho_{i,j}}{x}\right)_{i,j}^{\tau} \nonumber \\
% && - \nu \left[
%      \left(\frac{2\pi k}{x_{m}}\right)^{2p}+
%      \left(\frac{\pi l}{y_{m}}\right)^{2p}\right]\hat \zeta _{k,l}^{\tau}, 
%    \Deqlab{渦度の増分}
%    \\
% \hat G_{i,j}^{\tau} &=& -
%       \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         U_{i,j}
%         \left(\DP{\rho}{x}\right)_{i,j}^{\tau}  \nonumber \\
% && -  \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi_{i,j}}{x}\right)_{i,j}^{\tau}
%         \left(\DP{\rho _{i,j}}{y}\right)_{i,j}^{\tau}  \nonumber \\
% &&  + \frac{1}{I}\frac{1}{J} \sum _{i=0}^{I-1}\sum _{j=0}^{J-1}
%         \exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)
%         \cos\left(\frac{\pi l y_{j}}{y_{m}}\right)
%         \left(\DP{\psi_{i,j}}{y}\right)_{i,j}^{\tau}
%         \left(\DP{\rho _{i,j}}{x}\right)_{i,j}^{\tau}  \nonumber \\
% && - \kappa \left[
%      \left(\frac{2\pi k}{x_{m}}\right)^{2p}+
%      \left(\frac{\pi l}{y_{m}}\right)^{2p}\right]\hat \rho _{k,l}^{\tau}.
% \end{eqnarray}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section{使用モジュールとその他の設定}

スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の
et\_module に含まれる関数を用いて行う. spml が下位で使用す
る ISPACK の仕様から, 格子点数 $I,J$ は偶数で, かつ 
$I/2, J/2=2^{a}3^{b}5^{c}$ ($a, b, c$  は 0 または整数) でなければならな
い. 非線形項の計算によって生じるエリアジングを防ぐため, 格子点数 $I,J$ 
と切断波数 $K,L$ は  $I>3K,J>3K/2$ を満たすように与える.
時間方向の離散化は Euler スキームを用いて行う. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section{数値実験}

数値実験では水路中心で向きが反転する平行流を与え, ケルビン・ヘルムホルツ
不安定が生じる様子を計算する. 基本場の速度分布は以下のように与える.
\begin{equation}
  U(y) = U_{0} \tanh\left(\frac{y - y_{m}/2}{A_{0}}\right)
\end{equation}

初期の密度分布は
\begin{equation}
  \rho = \rho _{0} - \frac{\delta \rho}{2}
         \left[\tanh \left(\frac{y -y_{m}/2}{A_{0}}\right) + 1\right]
\end{equation}
とし, 流線関数は微小な振幅を持つ乱数として与える. 
パラメータは\Dtabref{使用したパラメータの値}にまとめた値を用いる.

格子点数 $I, J$ と切断波数 $K, L$ はそれぞれ $I=128, J=64, K=L=42$ とす
る. 時間格子間隔 $\Delta t$ は $0.25\times 10^{-4}$ sec, 計算ステップ数
は 40,000 ステップである.

\Dfigref{密度分布の時間変化}に計算された密度分布の時間変化の様子を示す.


\begin{table}[h]
 \begin{center}
  \begin{tabular}{c l l  } \hline
   パラメータ &\qquad\qquad&  数値   \\ \hline 
   $\rho_{0}$    && 1000 kgm${}^{-3}$  \\
   $\delta \rho$ && 50   kgm${}^{-2}$ \\
   $g$           && 9.8  msec${}^{-2}$   \\ 
   $p$           && 5       \\ 
   $\nu$         && 10${}^{-16}$ m${}^{2}$sec${}^{-1}$  \\ 
   $\kappa$      && 10${}^{-16}$ m${}^{2}$sec${}^{-1}$  \\ 
   $x_{m}$       && 0.18 m       \\ 
   $y_{m}$       && 0.06 m       \\ 
   $U_{0}$       && 0.3 msec${}^{-1}$ \\
   $A_{0}$       && 0.002 m       \\ \hline
 \end{tabular} 
\caption{使用したパラメータの値}
\Dtablab{使用したパラメータの値}
\end{center}
\end{table}

\begin{figure}[p]
\begin{center}
\Depsf[15cm][]{./figs/kh1.ps}
\vspace*{-1.0cm}
\caption{密度分布の時間変化. 
         左上から右下へと順に $t=$0.82, 0.84, 0.86, 0.88, 0.9, 0.92,
         0.94, 0.96, 0.98, 1.0 sec の結果. 寒色から暖色へと色が変化す
         るにつれて密度が高いことを示している.}  \Dfiglab{密度分布の時間
         変化}
\end{center}
\end{figure}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\newpage
\section{参考文献}
\begin{description}
 \item  Chandrasekhar, S., 1961: 
        Hydrodynamic and Hdromagnetic stability. Oxford University Press.

 \item GFD-online (酒井 敏, 飯澤 功, 荒巻 英治), 1997:
       実験室の中の空と海, 内部重力波, 
       http:\slash \slash www.gfd-dennou.org\slash arch\slash gfd-exp\slash gfd\_exp\slash exp\_j\slash doc\slash iw\slash guide01.htm.	     

 \item 竹広真一, 石岡圭一, 森川靖大, 小高正嗣, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2004: 
       階層的地球流体力学スペクトルモデル集 (SPMODEL),
       http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash ,
       地球流体電脳倶楽部. 
\end{description}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section*{謝辞}                           % 謝辞. このまま利用する.
\markright{謝辞}                          % ヘッダに表示
\addcontentsline{toc}{section}{謝辞}      % 目次に表示

本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用
の実験の一環として
\begin{center}
 http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 
\end{center}
において公開されているものである 
(\copyright 地球流体電脳倶楽部スペクトルモデルプロジェクト 
spmodel@gfd-dennou.org 2002. ). 
本資源は, 著作者の諸権利に抵触しない(迷惑をかけない)限りにおいて自由に
利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめ
ることをお願いする(無保証無責任原則).

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}


%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
