% 表題  2 次元 2 重周期境界領域でのβ面順圧流体モデルの定式化
%
% 履歴  2005/10/08 竹広 真一  新規作成
%       2005/10/14 竹広 真一  粘性解追加, 平行流修正
%       2005/10/21 竹広 真一  モドン解追加
%       
%
\documentclass[12pt,a4j,notitlepage]{jarticle}
\usepackage{Dennou6}
\usepackage{psfrag}
\usepackage{ascmac}
\usepackage{times}
\Dtitle[2 重周期領域順圧モデル]
       {2 次元 2 重周期境界領域での\\順圧流体モデルの定式化}
\Dauthor[竹広 真一]{竹広 真一, SPMODEL 開発グループ}
\Ddate{平成 20 年 10 月 10 日} 
\Dparskip

\pagestyle{Dheadings}

\begin{document}

\maketitle

\section{はじめに}

  この文章では, 2 次元 2 重周期境界領域でのβ面順圧流体モデルの定式化と
  テスト実験の問題設定について述べる. 

\section{支配方程式系}

  \subsection{支配方程式系}

    $x$, $y$ 方向に各々 $L_x$, $L_y$ の拡がりを持つ
    矩形領域での 2 次元非圧縮流体を考える. 
    支配方程式は, 2 次元β面上での運動方程式と連続の式である:
    %
    \begin{eqnarray}
      \Deqlab{運動方程式 x 成分}
      & & \DP{u}{t} + u\DP{u}{x}+ v\DP{u}{y} - f(y) v 
          = -\Dinv{\rho}\DP{p}{x} + F^{visc}_x, \\
      \Deqlab{運動方程式 y 成分}
      & & \DP{v}{t} + u\DP{v}{x}+ v\DP{v}{y} + f(y) u
          = -\Dinv{\rho}\DP{p}{y} + F^{visc}_y, \\
      \Deqlab{連続の式}
      & & \DP{u}{x}+ \DP{y}{y} = 0.
    \end{eqnarray}
    %
    ここで $u,v$ は速度の $x,y$ 成分, $p$ は圧力, 
    $\rho$ は密度であり定数である. 
    $f(y)$ がコリオリパラメター, 
    $F^{visc}_x,F^{visc}_y$ は粘性散逸項であり, 
    それぞれ具体的な形はここでは定めず一般的に定式化する. 

    境界条件は $x$, $y$ 両方向に周期的であるとする. 
    \begin{equation}
      u(0,y) = u(L_x,y), \ v(0,y) = v(L_x,y), \qquad
      u(x,0) = u(x,L_y), \ v(x,0) = v(x,L_y). 
    \end{equation}

  \subsection{渦度方程式(ポテンシャル渦度保存則)}

    \Deqref{運動方程式 x 成分}に $y$ 微分, 
    \Deqref{運動方程式 y 成分}に $x$ 微分を施して
    差引くことにより, 圧力を消去できて渦度の時間変化の式が得られる. 
    %
    \begin{eqnarray*}
      & &
        -\DP{}{y}\left[\DP{u}{t} + u\DP{u}{x}+ v\DP{u}{y} - f(y) v \right]
        +\DP{}{x}\left[\DP{v}{t} + u\DP{v}{x}+ v\DP{v}{y} + f(y) u \right]
      \\ & & \qquad
          = - \DP{F^{visc}_x}{y}+ \DP{F^{visc}_y}{x}, 
      \\& &
        \DP{}{t}\left(\DP{v}{x}-\DP{u}{y}\right)
          + u\DP{}{x}\left(\DP{v}{x}-\DP{u}{y}\right)
          + v\DP{}{y}\left(\DP{v}{x}-\DP{u}{y}\right) \\
       & &\quad
          - \DP{u}{y}\DP{u}{x}-\DP{v}{y}\DP{u}{y}
          + \DP{u}{x}\DP{v}{x}+\DP{v}{x}\DP{v}{y}
          - f(y)\left(\DP{u}{x}+\DP{v}{y}\right) +\DD{f(y)}{y}v 
      \\ & & \qquad
          = - \DP{F^{visc}_x}{y}+ \DP{F^{visc}_y}{x}, 
      \\& &
        \left(\DP{}{t}+u\DP{}{x}+v\DP{}{y}\right)
                  \left(\DP{v}{x}-\DP{u}{y}\right)
          + \left(\DP{v}{x}-\DP{u}{y}\right)
               \left(\DP{u}{x}+\DP{v}{y}\right)
      \\ & & \qquad
          = - \DP{F^{visc}_x}{y}+ \DP{F^{visc}_y}{x}, 
      \\& &
        \left(\DP{}{t}+u\DP{}{x}+v\DP{}{y}\right)
                  \left(\DP{v}{x}-\DP{u}{y}\right)
          +\DD{f(y)}{y}v 
            = - \DP{F^{visc}_x}{y}+ \DP{F^{visc}_y}{x}, 
    \end{eqnarray*}
    %
    途中, 連続の式\Deqref{連続の式}を用いた. 
    したがって, 渦度 $\Ddsty \zeta(x,y,t) = \DP{v}{x}-\DP{u}{y}$ を
    定義すると, 
    %
    \begin{equation}
      \Deqlab{渦度方程式}
        \left(\DP{}{t}+u\DP{}{x}+v\DP{}{y}\right)\zeta(x,y,t)
          +\DD{f(y)}{y}v 
            = F^{visc}_\zeta,
    \end{equation}
    %
    ここで $\Ddsty F^{visc}_\zeta =- \DP{F^{visc}_x}{y}+ \DP{F^{visc}_y}{x}$
    は渦度の散逸項である. 
    あるいは左辺をまとめて, 
    %
    \begin{equation}
      \Deqlab{ポテンシャル渦度保存則}
        \left(\DP{}{t}+u\DP{}{x}+v\DP{}{y}\right)[\zeta(x,y,t) + f(y)]
            = F^{visc}_\zeta.
    \end{equation}
    %
    散逸がない場合に $\zeta + f(y)$ は保存量となる. 
    $q(x,y,t) \equiv \zeta(x,y,t) + f(y)$ は系の絶対渦度あるいは
    ポテンシャル渦度である. これに対して $\zeta(x,y,t)$ は相対渦度, 
    $f(y)$ は座標系のもつ渦度
    \footnote{気象学では惑星渦度}と呼ばれる. 

  \subsection{流線関数を用いた表現}

    連続の式\Deqref{連続の式}から速度場を流線関数
    $\psi(x,y,t)$ を用いて表すことができる. 
    すなわち, 
    %
    \begin{equation}
      \Deqlab{速度と流線関数}
      u = -\DP{\psi}{y}, \quad v = \DP{\psi}{x}
    \end{equation}
    %
    とすると, 相対渦度は
    %
    \begin{equation}
      \zeta = \DP{v}{x}-\DP{u}{y} =\DP[2]{\psi}{x}+\DP[2]{\psi}{y}
             =\Dlapla{\psi}.
    \end{equation}
    %
    ここで $\Ddsty \Dlapla=\DP[2]{}{x}+\DP[2]{}{y}$ は 2 次元の
    ラプラシアン演算子である. 
    %
    流線関数を用いて渦度方程式\Deqref{渦度方程式}を表すと
    %
    \begin{equation}
        \DP{\zeta}{t}-\DP{\psi}{y}\DP{\zeta}{x}+\DP{\psi}{x}\DP{\zeta}{y}
          +\DD{f(y)}{y}\DP{\psi}{x}
            = F^{visc}_\zeta.
    \end{equation}
    %
    すなわち, 
    %
    \begin{displaymath}
      \Deqlab{渦度方程式(流線関数)}
        \DP{\zeta}{t} + J(\psi,\zeta) +\beta(y)\DP{\psi}{x}
            = F^{visc}_\zeta, 
       \quad \zeta(x,y,t) = \Dlapla\psi(x,y,t).
    \end{displaymath}
    %
    ここで $\Ddsty J(f,g)\equiv \DP{f}{x}\DP{g}{y}-\DP{f}{y}\DP{g}{x}$ は
    ヤコビアン演算子, 
    $\Ddsty \beta(y)\equiv \DD{f(y)}{y}$ は
    コリオリパラメターの $y$ 微分である. 
    ポテンシャル渦度保存則\Deqref{ポテンシャル渦度保存則}は
    %
    \begin{eqnarray}
      \Deqlab{ポテンシャル渦度保存則(流線関数)}
        \DP{q}{t}+ J(\psi,q)
           & =& F^{visc}_\zeta, \\
        q(x,y,t) &=& \zeta(x,y,t) + f(y) = \Dlapla\psi(x,y,t)+f(y).
    \end{eqnarray}

    境界条件は $x$, $y$ 両方向に周期的であるので, 
    %
    \begin{eqnarray}
      & & \zeta(0,y) = \zeta(L_x,y), \quad \zeta(x,0) = \zeta(x,L_y), \\
      & & \psi(0,y) = \psi(L_x,y), \quad \psi(x,0) = \psi(x,L_y). 
    \end{eqnarray}
    %
    あたえるコリオリパラメターも $y$ 方向に周期的でなければならない. 
    しかしながら\Deqref{渦度方程式(流線関数)}には $f(y)$ でなく
    その微分 $\beta(y)$ しか方程式にあらわれないので
    %
    \begin{equation}
      \beta(0) = \beta(L_y),
    \end{equation}
    %
    であれば解くことができる
    \footnote{$f(y)$ が周期的でない場合には圧力を求めようとすると破綻する. 
      この場合ポテンシャル渦度保存則を用いて計算を行なえない. }. 

  \subsection{粘性散逸の表現}

    粘性散逸項として, 超粘性型のものを考えることにする. 
    %
    \begin{equation}
      F^{visc}_x = (-1)^{p+1}\nu_{2p}\Dlapla[2p] u, \quad
      F^{visc}_x = (-1)^{p+1}\nu_{2p}\Dlapla[2p] v,
    \end{equation}
    %
    $p$ は超粘性の次数を表す整数パラメターである. 
    このとき渦度に対する散逸は, 
    %
    \begin{equation}
      F^{visc}_\zeta = (-1)^{p+1}\nu_{2p}\Dlapla[2p] \zeta.
    \end{equation}

  \subsection{支配方程式(まとめ)}

    以上, 支配方程式をまとめると, 
    %
    \begin{eqnarray}
      \Deqlab{流線関数渦度方程式}
        & & \DP{\zeta}{t} + J(\psi,\zeta) +\beta(y)\DP{\psi}{x}
              = (-1)^{p+1}\nu_{2p}\Dlapla[2p] \zeta, \\
        & &  \zeta(x,y,t) = \Dlapla\psi(x,y,t),
    \end{eqnarray}
    %
    境界条件は
    %
    \begin{eqnarray}
      \Deqlab{渦度流線関数境界条件}
        & & \zeta(0,y) = \zeta(L_x,y), \quad \zeta(x,0) = \zeta(x,L_y), \\
        & & \psi(0,y) = \psi(L_x,y), \quad \psi(x,0) = \psi(x,L_y). 
    \end{eqnarray}
    %
    ただし, 与える物理パラメター $\beta(y)$ は
    %
    \begin{equation}
      \Deqlab{ベータパラメター境界条件}
      \beta(0) = \beta(L_y),
    \end{equation}
    %
    を満たしていなければならない. 
    $\beta$ が一定の場合がいわゆる $\beta$ 面モデルと呼ばれるものである. 

\section{数値計算 : 2 重フーリエ級数展開と時間積分}

  \subsection{2 重フーリエ展開}

    スペクトル法による数値計算を行うために, 
    各物理量を境界条件を満たす函数系である 2 重フーリエ級数展開する. 
    すなわち, 
    %
    \begin{eqnarray}
      \zeta(x,y,t)
         &=& \sum_{k=-K}^{K} \sum_{l=-L}^{L}
               \tilde{\zeta}_{kl}(t)
                 \exp\left(i\frac{2\pi x}{L_x}k\right)
                 \exp\left(i\frac{2\pi y}{L_y}l\right),\\
      \tilde{\zeta}_{kl}(t)
         &=& \Dinv{L_x L_y}\int_0^{L_x}\int_0^{L_y}\zeta(x,y,t)
                 \exp\left(-i\frac{2\pi k}{L_x}x\right)
                 \exp\left(-i\frac{2\pi l}{L_y}y\right) dx \ dy,
    \end{eqnarray}
    %
    $\tilde{\zeta}_{kl}(t)$ は渦度のスペクトルの波数 $k,l$ 成分, 
    $K,L$ はそれぞれ $x,y$ 方向の切断波数である. 
    流線関数 $\psi(x,y,t)$ のスペクトルも同様に 
    $\tilde{\psi}_{kl}(t)$ を定義する. 

    これらを\Deqref{流線関数渦度方程式}に代入し, 
    $\Ddsty \exp\left(-i\frac{2\pi k'}{L_x}x\right)
            \exp\left(-i\frac{2\pi l'}{L_y}y\right)$
    をかけて全領域で積分することにより, 
    基底関数の直交関係から各スペクトル成分の式が得られる. 
    %
    \begin{displaymath}
        \DD{\tilde{\zeta}_{kl}}{t} 
           = - \widetilde{[J(\psi,\zeta)]}_{kl} 
             - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}
             + (-1)^{p+1}\nu_{2p}\widetilde{[\Dlapla[2p]\zeta]}_{kl}.
    \end{displaymath}
    %
    ただし $[\ldots]_{kl}$ はスペクトルの $k,l$ 成分であることを表している. 
    %
    \begin{displaymath}
        \DD{\tilde{\zeta}_{kl}}{t} 
           = - \widetilde{[J(\psi,\zeta)]}_{kl} 
             - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}
             + (-1)^{p+1}\nu_{2p}\widetilde{[\Dlapla[2p]\zeta]}_{kl}.
    \end{displaymath}
    %
    散逸項はラプラシアンが 2 重指数函数系の固有演算子であることから
    %
    \begin{eqnarray*}
      (-1)^{p+1}\nu_{2p}\widetilde{[\Dlapla[2p]\zeta]}_{kl}
      &=& (-1)^{p+1}\nu_{2p}
         \left[-\left(\frac{2\pi k}{L_x}\right)^2
               -\left(\frac{2\pi l}{L_y}\right)^2\right]^p 
          \cdot \tilde{\zeta}_{kl}\\
      &=& -\nu_{2p}\left[ \left(\frac{2\pi k}{L_x}\right)^2
                         +\left(\frac{2\pi l}{L_y}\right)^2\right]^p
           \cdot  \tilde{\zeta}_{kl} 
       = -\nu^*(k,l,p) \ \tilde{\zeta}_{kl}
    \end{eqnarray*}
%
    となる. ここで, 波数毎の超粘性係数 $\nu^*(k,l,p)$ を
    \begin{equation}
        \nu^*(k,l,p)\equiv \nu_{2p}
                            \left[ \left(\frac{2\pi k}{L_x}\right)^2
                                  +\left(\frac{2\pi l}{L_y}\right)^2\right]^p
    \end{equation}
    %
    と定義した($\nu^*(k,l,p)$ は 0 以上であることに注意). 
    したがって, 
    %
    \begin{equation}
      \Deqlab{流線関数渦度方程式スペクトル}
        \DD{\tilde{\zeta}_{kl}}{t} 
           = - \widetilde{[J(\psi,\zeta)]}_{kl} 
             - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}
             - \nu^*(k,l,p)\tilde{\zeta}_{kl}.
    \end{equation}

  \subsection{時間積分}

    以下では時間に関して差分化し, 時間積分を行うための定式化を記す. 
    $\Delta t$ を時間格子間隔, 時刻 $t, t+\Delta t, t-\Delta t$ 
    における $\tilde{\zeta}_{kl}$ の値をそれぞれ
    $\tilde{\zeta}_{kl}^{\tau}$, 
    $\tilde{\zeta}_{kl}^{\tau+1}$, 
    $\tilde{\zeta}_{kl}^{\tau-1}$ 等と表す. 

    \subsubsection{euler スキーム(1 次精度)}

      euler スキームの場合は時間変化項をすべて現在時間で評価して, 
      %
      \begin{eqnarray*}
         \frac{\tilde{\zeta}_{kl}^{\tau+1}-\tilde{\zeta}_{kl}^{\tau}}
              {\Delta t} 
          = - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
            - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
            - \nu^*(k,l,p)\tilde{\zeta}_{kl}^\tau.
      \end{eqnarray*}
      %
      すなわち, 
      % 
      \begin{eqnarray}
         \tilde{\zeta}_{kl}^{\tau+1}
          = \tilde{\zeta}_{kl}^{\tau}
            + (\Delta t) \times \left\{
                 - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
                 - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
                 - \nu^*(k,l,p)\tilde{\zeta}_{kl}^\tau
              \right\}.
      \end{eqnarray}

    \subsubsection{Leap frog and Implicit スキーム}

      散逸項を陰解法(implicit  スキーム), 
      それ以外の項を Leap frog スキームで評価する. 
      まず第一段階として散逸項以外の項に関して積分する. 
      %
      \begin{displaymath}
         \frac{\tilde{\zeta}_{kl}^{\tau+1}-\tilde{\zeta}_{kl}^{\tau-1}}
              {2(\Delta t)} 
          = - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
            - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
            - \nu^*(k,l,p)\tilde{\zeta}_{kl}^{\tau+1}
      \end{displaymath}
      %
      変形すると
      %
      \begin{displaymath}
         [1+2\nu^*(k,l,p)(\Delta t)]\tilde{\zeta}_{kl}^{\tau+1}
          = \tilde{\zeta}_{kl}^{\tau-1}
             + 2(\Delta t) \times \left\{
                  - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
                  - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
               \right\}.
      \end{displaymath}
      %
      したがって
      %
      \begin{equation}
         \tilde{\zeta}_{kl}^{\tau+1}
          = \frac{\Ddsty 
               \tilde{\zeta}_{kl}^{\tau-1}
             + 2(\Delta t) \times \left\{
                  - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
                  - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
               \right\}}
            {1+2\nu^*(k,l,p)(\Delta t)}
      \end{equation}

      Leap frog スキームを用いているので, 
      計算モードの増幅を抑制するために Asselin (1972) の
      時間フィルターをとあるタイムステップ間隔で適用する. 
      フィルターの定義式は
      %
      \begin{equation}
        A^{\tau*} = A^{\tau} 
         + \gamma (A^{\tau+1}-2A^{\tau}+A^{\tau-1})
      \end{equation}
      %
      ここで $A^{\tau*}$ が修正された $\tau$ での $A$ の値, 
      $\gamma $ はフィルターの係数である. 

           
    \subsubsection{Adams-Bashforth and Crank-Nicolson スキーム(2 次精度)}

      散逸項に Clank-Nicolson スキーム, それ以外の項に
      2 次の Adams-Bashforth スキームを適用する. 
      %
      \begin{eqnarray*}
        \frac{\tilde{\zeta}_{kl}^{\tau+1}-\tilde{\zeta}_{kl}^\tau}{\Delta t}
         &=&  \frac{3}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^\tau
             -\frac{1}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^{\tau-1}
          - \nu^*(k,l,p)
                \frac{\tilde{\zeta}_{kl}^{\tau+1}+\tilde{\zeta}_{kl}^{*\tau+1}}
                     {2},
     \end{eqnarray*}
     %
     ただし $\Ddsty\left(\DD{\zeta}{t}\right)_{nd}$ は 
     時間変化の非散逸項を表している:
     %
     \begin{eqnarray}
       \left(\DD{\zeta}{t}\right)_{nd}
        &\equiv&
            - J(\psi,\zeta) - \beta(y)\DP{\psi}{x}, \\
        \left[\widetilde{
              \left(\DD{\zeta}{t}\right)_{nd}}\right]_{kl}^\tau 
        &=&
            - \widetilde{[J(\psi,\zeta)]}_{kl}^\tau
            - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}^\tau
     \end{eqnarray}
     %
     これを解くと, 
     %
     \begin{eqnarray*}
       [1 + \nu^*(k,l,p)(\Delta t)/ 2]\tilde{\zeta}_{kl}^{\tau+1} 
       &=& 
       [1 - \nu^*(k,l,p)(\Delta t)/2]\tilde{\zeta}_{kl}^\tau \\
       & &   + \Delta t \times \left\{
              \frac{3}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^\tau
             -\frac{1}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^{\tau-1}
              \right\}, \\
     \end{eqnarray*}
     %
     すなわち, 
     %
     \begin{eqnarray}
       \tilde{\zeta}_{kl}^{\tau+1} 
         &=& \frac{\Ddsty
              [1 - \nu^*(k,l,p)(\Delta t)/2]\tilde{\zeta}_{kl}^\tau 
             + \Delta t \times \left\{
              \frac{3}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^\tau
             -\frac{1}{2}\left[
                         \widetilde{\left(\DD{\zeta}{t}\right)}_{nd}
                         \right]_{kl}^{\tau-1}
              \right\}}
             {1 + \nu^*(k,l,p)(\Delta t)/ 2}.\nonumber\\
     \end{eqnarray}

    \subsubsection{演算子分割処理法(粘性項)と 4 次精度 Runge-Kutta スキーム}

      変数変換することによって粘性項を解析的に評価し, 
      その他の項を高次の数値積分スキームを適用する
      演算子分割法をとる場合を以下に示す. 
      渦度のスペクトル成分について
      %
      \begin{equation}
        \tilde{\zeta}_{kl}(t) = \hat{\zeta}_{kl}(t) e^{-\nu^*(k,l,p)t}
      \end{equation}
      %
      と変数変換すれば, 
      %
      \begin{displaymath}
        \DD{\tilde{\zeta}_{kl}}{t} 
         =   \DD{\hat{\zeta}_{kl}(t)}{t} e^{-\nu^*(k,l,p)t} 
           - \nu^*(k,l,p)\hat{\zeta}_{kl}(t) e^{\nu^*(k,l,p)t}
         =   \DD{\hat{\zeta}_{kl}(t)}{t} e^{-\nu^*(k,l,p)t} 
           - \nu^*(k,l,p)\tilde{\zeta}_{kl}(t)
      \end{displaymath}
      %
      であるから, \Deqref{流線関数渦度方程式スペクトル}は
      %
      \begin{equation}
        \DD{\hat{\zeta}(t)_{kl}}{t} e^{-\nu^*(k,l,p)t} 
        =  - \widetilde{[J(\psi,\zeta)]}_{kl} 
           - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}.
      \end{equation}
      %
      \begin{equation}
        \DD{\hat{\zeta}(t)_{kl}}{t} 
        = e^{\nu^*(k,l,p)t}\left[ 
           - \widetilde{[J(\psi,\zeta)]}_{kl} 
           - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}
          \right].
      \end{equation}
      %
      この左辺の時間積分を Runge-Kutta スキームで計算する. 
      すなわち, 
      \begin{eqnarray*}
        \tilde{\zeta}_{kl}^{\tau+1} 
          &=& \hat{\zeta}_{kl}^{\tau+1} e^{-\nu^*(k,l,p)\Delta t},\\
        \hat{\zeta}_{kl}^{\tau+1} 
          &=& \tilde{\zeta}_{kl}^\tau
           + (\Delta t) \times \left(  \frac{k_1}{6} +\frac{k_2}{3} 
                                   +\frac{k_3}{3} +\frac{k_4}{6} \right), \\
        k_1 &=& f[\zeta^\tau,\psi^\tau,0]_{kl},\\
        k_2 &=& f[\zeta^{\tau+1/2},\psi^{\tau+1/2},\Delta t/2]_{kl}\\
        k_3 &=& f[\zeta^{*\tau+1/2},\psi^{*\tau+1/2},\Delta t/2]_{kl}\\
        k_4 &=& f[\zeta^{*\tau+1},\psi^{*\tau+1},\Delta t]_{kl}\\
        f[\zeta,\psi,t]_{kl} &=& 
          e^{\nu^*(k,l,p)t}\left[ 
           - \widetilde{[J(\psi,\zeta)]}_{kl} 
           - \widetilde{\left[\beta(y)\DP{\psi}{x}\right]}_{kl}
          \right],\\
        \tilde{\zeta}_{kl}^{\tau+1/2} &=& e^{-\nu^*(k,l,p)\Delta t/2}
             [\tilde{\zeta}_{kl}^\tau + k_1 \Delta t/2],\\
        \tilde{\zeta}_{kl}^{*\tau+1/2} &=& e^{-\nu^*(k,l,p)\Delta t/2}
            [\tilde{\zeta}_{kl}^\tau + k_2 \Delta t/2],\\
        \tilde{\zeta}_{kl}^{*\tau+1} &=& e^{-\nu^*(k,l,p)\Delta t}
            [\tilde{\zeta}_{kl}^\tau + k_3 \Delta t].
      \end{eqnarray*}
      %
      この定式化において 
      $\tilde{\zeta}_{kl}^\tau = \hat{\zeta}_{kl}^\tau$ であることを
      用いていることに注意されたい. 

   \subsubsection{演算子分割処理法(線形項)と 4 次精度 Runge-Kutta スキーム}

      $\beta$ が一定の場合($\beta$ 平面)には, 粘性項だけでなく
      $\beta$ 項も線形となるので演算子分割処理を施すことができる. 
      
      $\beta$ が一定の場合には, 
      $\beta$ 項を渦度を用いて書き直すことができて, 
      %
      \begin{displaymath}
        \widetilde{\beta\left[\DP{\psi}{x}\right]}_{kl}
        = \beta\widetilde{\left[\DP{\psi}{x}\right]}_{kl}
        = -i \beta\left(\frac{2\pi k}{L_x}\right)
                \left[ \left(\frac{2\pi k}{L_x}\right)^2
                       +\left(\frac{2\pi l}{L_y}\right)^2\right]^{-1}
          \cdot\tilde{\zeta}_{kl}
        = i \omega_R(\beta,k,l)\cdot\tilde{\zeta}_{kl}
      \end{displaymath}
      %
      ここで $\omega_R(\beta,k,l)$ はロスビー波の位相速度
      \begin{equation}
        \omega_R(\beta,k,l) 
           \equiv -\beta\left(\frac{2\pi k}{L_x}\right)
                     \left[ \left(\frac{2\pi k}{L_x}\right)^2
                       +\left(\frac{2\pi l}{L_y}\right)^2\right]^{-1}
      \end{equation}
      %
      である. 
      \Deqref{流線関数渦度方程式スペクトル}は
      %
      \begin{equation}
        \Deqlab{流線関数渦度方程式スペクトル}
          \DD{\tilde{\zeta}_{kl}}{t} 
             = - \widetilde{[J(\psi,\zeta)]}_{kl} 
               - i \omega_R(\beta,k,l)\tilde{\zeta}_{kl}
               - \nu^*(k,l,p)\tilde{\zeta}_{kl}.
      \end{equation}
      %
      線形項を消去すべく 
      $ \tilde{\zeta}_{kl}(t) 
         = Z_{kl}(t)\cdot e^{-[i\omega_R(\beta,n,m)+\nu^*(k,l,p)]t}$
      と変換すれば, 
      %
      \begin{eqnarray*}
        \DD{\tilde{\zeta}_{kl}}{t} 
         &=& \DD{Z_{kl}}{t} e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}
               -[i\omega_R(\beta,k,l)+\nu^*(k,l,p)]Z_{kl}(t) 
                    e^{-[i\omega_R(\beta,k,l)+\nu^*(k,l,p)]t}\\
         &=&  \DD{Z_{kl}}{t} e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}
                -[i\omega_R(\beta,k,l) + \nu^*(k,l,p)]\tilde{\zeta}_{kl},
      \end{eqnarray*}
      したがって
      \begin{displaymath}
        \DD{Z_{kl}}{t}
        = e^{[i\omega_R(\beta,k,l)+\nu^*(k,l,p)]t}G[\psi,\zeta]_{kl}
      \end{displaymath}
      %
      ここで $G[\psi,\zeta]_{kl}$ は非線形項
      \begin{displaymath}
        G[\psi,\zeta]_{kl} = -\widetilde{[J(\psi,\zeta)]}_{kl}
      \end{displaymath}
      %
      を表す. 
      %
      ISPACK のスペクトル成分の格納のし方にあわせて
      配列を実部と虚部にわけて表すと, 
      $k > 0 $に対して
      $Re[Z_{kl}]=\xi_{kl} ならば Im[Z_{kl}]=\xi_{-k-l}$ 
      であるから
      %
      \begin{eqnarray}
         \DD{\xi_{kl}}{t}
         &=&  e^{\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t] 
                Re[G[\psi,\zeta]_{kl}]\nonumber\\
         & & - e^{\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t]
                Im[G[\psi,\zeta]_{kl}]\nonumber\\
         &=&  e^{\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t]
                  G[\psi,\zeta]_{kl}
             - e^{\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t]
                  G[\psi,\zeta]_{-k-l}, \\
         \DD{\xi_{-k-l}}{t}
         &=&   e^{\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t] 
                Im[G[\psi,\zeta]_{kl}] \nonumber\\
         & & + e^{\nu^*(k,l,p)]t}\sin[\omega_R(\beta,k,l) t] 
                Re[G[\psi,\zeta]_{kl}] \nonumber\\
         &=&   e^{\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t] 
                  G[\psi,\zeta]_{-k-l}
             + e^{\nu^*(k,l,p)]t}\sin[\omega_R(\beta,k,l) t]
                  G[\psi,\zeta]_{kl}
      \end{eqnarray}
      %
      時間積分のための $\xi_{kl}$, $\xi_n^{-m}$ の各積分ステップでの初期値は
      $\xi_{kl}$ と $Z_{kl}$ の関係において $t=0$ を適用すればよく, 
      $\xi_{kl} = \tilde{\zeta}_{kl}$, 
      $\xi_{-k-l} = \tilde{\zeta}_{-k-l}$
      となる. 上の式を用いて数値積分し, あたらしい時間の
      $\xi_{kl}$, $\xi_{-k-l}$ が求まれば, 
      対応して $\tilde{\zeta}_{kl}$ がつぎのように求められる. 
      %
      \begin{eqnarray}
        \tilde{\zeta}_{kl} 
            &=& Re[Z_{kl} e^{-[i\omega_R(\beta,k,l)+\nu^*(k,l,p)]t}]\nonumber\\
            &=&   Re[Z_{kl}] e^{-\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t]
                + Im[Z_{kl}] e^{-\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t]
                \nonumber\\
            &=&   \xi_{kl} e^{-\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t]
                + \xi_{-k-l} e^{-\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t],\\
        \tilde{\zeta}_{-k-l}
            &=& Im[Z_{kl} e^{-[i\omega_R(\beta,k,l)+\nu^*(k,l,p)]t}]\nonumber\\
            &=&   Im[Z_{kl}] e^{-\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t]
                - Re[Z_{kl}] e^{-\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t]
                \nonumber\\
            &=&   \xi_{-k-l} e^{-\nu^*(k,l,p)t}\cos[\omega_R(\beta,k,l) t]
                - \xi_{kl} e^{-\nu^*(k,l,p)t}\sin[\omega_R(\beta,k,l) t].
      \end{eqnarray}


 \section{テスト実験}


   \subsection{テスト 1 : 線形項のテスト}

       単一モード波は非線形項が 0 となるので, 
       線形方程式の分散関係にしたがう解析解となる(ロスビー波解). 
       %
       \begin{equation}
         \psi=A_0\exp(i\tilde{k}x + i\tilde{l}y - \omega t),
       \end{equation}
       %
       を考える. ここで
       \begin{eqnarray*}
         \tilde{k} = \frac{2\pi k}{L_x}, \quad 
         \tilde{l} = \frac{2\pi l}{L_y}
       \end{eqnarray*}
       %
       は水平波数であり, $k, l$ は整数である. 
       %
       \begin{displaymath}
         \zeta = \Dlapla\psi=-(\tilde{k}^2 + \tilde{l}^2)
                               A_0\exp(i\tilde{k}x + i\tilde{l}y - \omega t)
                            =-(\tilde{k}^2 + \tilde{l}^2)\psi,
       \end{displaymath}
       %
       であるから, 非線形項が
       %
       \begin{displaymath}
         J(\psi,\zeta)=J(\psi,A\psi) = 0,
       \end{displaymath}
       %
       となる. よって, 残りの線形項に解を代入して分散関係が得られ, 
       %
       \begin{displaymath}
         \omega(\beta,k,l) 
              = -\frac{\beta \tilde{k}}{\tilde{k}^2+\tilde{l}^2}
                -i\nu_{2p}(\tilde{k}^2+\tilde{l}^2)^p              
       \end{displaymath}
       %
       となる. 即ち解は
       %
       \begin{eqnarray*}
         \psi &=& A_0\exp\{i[\tilde{k}x + \tilde{l}y - \omega(\beta,k,l)] t\}
                  \exp[-\nu_{2p}(\tilde{k}^2+\tilde{l}^2)^p t]\\
              &=& \psi(t=0)\exp[- i\omega(\beta,k,l) t]
                       \exp[-\nu_{2p}(\tilde{k}^2+\tilde{l}^2)^p t],
       \end{eqnarray*}
       %
       である. ただし $\psi(t=0)=A_0\exp[i(\tilde{k}x + \tilde{l}y)]$ は
       初期値流線関数である. 
       %
       実数関数で書くならば, 
       \begin{eqnarray*}
         \psi =  \exp[-\nu_{2p}(\tilde{k}^2+\tilde{l}^2)^p t]
                &\cdot&
                       [ Re[\psi(t=0)]\cos \omega(\beta,k,l)t
                       + Im[\psi(t=0)]\sin\omega(\beta,k,l)t ] \\
             & &   + i [ Im[\psi(t=0)]\cos\omega(\beta,k,l)t
                       - Re[\psi(t=0)]\sin\omega(\beta,k,l)t ] \}
       \end{eqnarray*}
       %
       振幅 $A_0$ は任意であるが, 
       粘性がない場合に振幅が大きくなると
       ロスビー波解が不安定となることが知られている
       ({\it e.g.} Hoskins and Hollingworth, 1972). 

   \subsection{テスト 2 : 移流項のテスト}

      まずは $\beta=0$ の場合に, 渦が形を変えずに移流される場合の解を
      正しく計算できるかチェックする. 
      
      そのための解析解として, 
      「定常な渦を平行流あるいは剛体回転流に重ねた場」を
      考える. このような場では定常な渦がその形を変えず
      平行流あるいは剛体回転流によって流されて行くだけである. 
      以下そのことを示しておく. 

      今, 定常な渦の流線関数を $\psi'$ とすると, 定常解であることから
      %
      \begin{displaymath}
        J(\psi',\Dlapla\psi') = 0, 
      \end{displaymath}
      %
      である. 一方, 基本流の流線関数を $\psi_0$ とすると, 
      平行流あるいは剛体回転流なので $\Dlapla\psi_0=$const. である. 
      したがって, これらを重ねた場の時間発展は
      %
      \begin{eqnarray*}
        \DP{}{t}(\psi_0+\psi') &=& J[\psi_0+\psi',\Dlapla(\psi_0+\psi')]\\
         &=& J(\psi_0,\Dlapla\psi_0) + J(\psi_0,\Dlapla\psi')
            + J(\psi',\Dlapla\psi_0) + J(psi',\Dlapla\psi'),\\
        \DP{}{t}(\psi_0+\psi') &=& J(\psi_0,\Dlapla\psi')
      \end{eqnarray*}
      %
      左辺は $\psi'$ の渦の流れ $\psi_0$ による移流を表している. 
      基本流が平行流あるいは剛体回転ならば
      $\psi_0$ の流れによって移流されても渦の形が変わらないので, 
      上式がつねに成立ち, 渦が形を変えずに流されていくだけとなる. 

      局所的な渦で定常なものは軸対称な流線および渦度分布を持つものである. 
      そのような渦の流線と渦度は渦の中心からの距離 $r$ だけで表すことが
      できて
      %
      \begin{displaymath}
        \psi'= f(r), \Dlapla\psi'=g(r).
      \end{displaymath}
      %
      これに対する移流項のヤコビアンは
      \begin{eqnarray*}
        J(\psi',\Dlapla\psi') &=&   \DP{\psi'}{x}\DP{(\Dlapla\psi')}{y}
                                - \DP{\psi'}{x}\DP{(\Dlapla\psi')}{y}\\
            &=&   \DP[][y]{\ r}{x}\DD{f}{r}\DP[][x]{\ r}{y}\DD{g}{r}
              - \DP[][y]{\ r}{y}\DD{f}{r}\DP[][y]{\ r}{x}\DD{g}{r}
            = 0.
      \end{eqnarray*}
      %
      したがって定常である. 

     \subsubsection{平行流とガウス型の渦}
     
       平行流の基本流に局所的な軸対称のガウス分布の渦を初期値として与える. 
       しかしながら, 
       今適用しているスペクトル法では, 渦度および流線関数の
       $(k,l)=(0,0)$ スペクトル成分は 0 としなければならない. 
       そこで平行流の代わりに領域の半分づつが互いに反平行な流れとなっている
       状況を与えることにする. 

       すなわち
       \begin{displaymath}
          U_x = \left\{
          \begin{array}{lc}
            \Ddsty
             U_0 \tanh \left(\frac{y-L_y/4}{\sigma}\right)  
                       &\Ddsty\left( 0\ge y\ge \frac{L_y}{2} \right),\\
            \Ddsty
            -U_0 \tanh \left(\frac{y-3L_y/4}{\sigma}\right) 
                       &\Ddsty \left( \frac{L_y}{2}\ge y\ge L_y \right),
          \end{array}
          \right.
          \qquad U_y = 0.
       \end{displaymath}
       %
       $U_0$ が流れの振幅, $\sigma$ が反平行な流れの境界層の厚さを表す. 
       対応する流線関数は,
       \begin{displaymath}
          \psi_0(x,y) = \left\{
          \begin{array}{lc}
            \Ddsty
             -U_0 \sigma \left[
                \log\left|\cosh \left(\frac{y-L_y/4}{\sigma}\right)\right|  
                 -\log\left|\cosh \left(\frac{L_y/4}{\sigma}\right)\right|
                \right]
                       & \Ddsty\left( 0\ge y\ge \frac{L_y}{2} \right),\\
            \Ddsty
             U_0 \sigma \left[
                 \log\left|\cosh \left(\frac{y-3L_y/4}{\sigma}\right)\right|
                 -\log\left|\cosh \left(\frac{L_y/4}{\sigma}\right)\right|
                 \right]
                       & \Ddsty\left( \frac{L_y}{2}\ge y\ge L_y \right), 
          \end{array}
          \right.
       \end{displaymath}
       %
       一方, 重ねる局所渦は
       %
%        \begin{equation}
%           \psi'(x,y,0) = A_0[ \exp
%                          \left(-\frac{(x-x_0)^2+(y-y_0)^2}{a^2}\right)
%                          - \pi a^2]
%        \end{equation}
       \begin{equation}
           \psi'(x,y,0) = A_0[1- (x-x_0)^2+(y-y_0)^2)] \exp
                         \left(-\frac{(x-x_0)^2+(y-y_0)^2}{a^2}\right),
       \end{equation}
       %
       $A_0$ が渦の振幅, $a$ が渦の特徴的な大きさを表すパラメターである.
       第 2 項目は一様渦度成分を 0 とするための定数である. 
       $y_0$ を中心付近にとり, 渦の大きさを十分に小さくすれば, 
       渦は一様流 $U_0$ によって流される. このときの解は
       %
        \begin{equation}
           \psi'(x,y,0) = A_0[1- (x-x_0-U_0 t)^2+(y-y_0)^2)] \exp
                          \left(-\frac{(x-x_0-U_0 t)^2+(y-y_0)^2}{a^2}\right),
        \end{equation}
       %
       となる. 

       注意: 基本流がそもそも不安定(順圧 or K-H 不安定)な状況ので, 
       基本場が大きくくずれるまでの間だけ渦が形を変えずに移流される
       様子が観察できれば OK とする. 

%      \subsubsection{剛体回転流ガウス型の渦}

%        剛体回転流の基本流に局所的な軸対称のガウス分布の渦を初期値として与える. 
%        すなわち
%        %
%        \begin{equation}
%           \psi_0(x,y)  = -\Dinv{2}\omega(x^2 + y^2), 
%           \psi'(x,y,0) = A_0 \exp
%                          \left(-\frac{(x-x_0)^2+(y-y_0)^2}{a^2}\right),
%        \end{equation}
%        %
%        このとき, 渦は角速度 $\omega$ の剛体回転流で流されて行くので, 
%        その中心の座標 $(x_c(t),y_c(t))$ が
%        %
%        \begin{equation}
%          \left(
%            \begin{array}{c}
%              x_c(t)\\y_c(t) 
%            \end{array}
%          \right) 
%          = 
%          \left(
%            \begin{array}{cc}
%              \cos\omega t & -\sin\omega t \\
%              \sin\omega t &  \cos\omega t 
%            \end{array}
%          \right) 
%          \left(
%            \begin{array}{c}
%              x_0\\y_0
%            \end{array}
%          \right)
%          =
%          \left(
%            \begin{array}{c}
%              x_0\cos\omega t - y_0\sin\omega t \\ 
%              x_0\sin\omega t + y_0\cos\omega t 
%            \end{array}
%          \right).
%        \end{equation}
%        %
%        となり, 流線関数の時間発展は
%        %
%        \begin{equation}
%           \psi(x,y,t)  = \psi_0(x,y) +
%                A_0 \exp\left(-\frac{[x-x_c(t)]^2+[y-y_c(t)]^2}
%                                    {a^2}\right),
%        \end{equation}
%        %
%        となる. 

   \subsection{テスト 3 : モドン解}

     $\beta$ 一定の無限領域ベータ平面では孤立渦の解析解である
     モドンが存在する. その形は渦の中心を原点とする平面極座標
     $(r,\theta)$ を用いて次のように表される. 
     %
     \begin{equation}
       \zeta(r,\theta) = \Dlapla\psi 
          = - \frac{c q^2}{a}\sin\theta \times \left\{
       \begin{array}{lc}
         \Ddsty \frac{{\rm J}_1(kr/a)}{{\rm J}_1(k)} & (r < a) \\
         \Ddsty \frac{{\rm K}_1(qr/a)}{{\rm q}_1(q)} & (r \ge a)
       \end{array}
       \right.
     \end{equation}
     %
     ここで $a$ $c$ はそれぞれ孤立渦の半径と伝播速度, 
     パラメター $q = \sqrt{\beta/c}$, $k$ は分散関係に相当する関係式
     %
     \begin{equation}
       -\frac{{\rm J}_2(k)}{k{\rm J}_1(k)} = 
       \frac{{\rm K}_2(q)}{q{\rm K}_1(q)} 
     \end{equation}
     %
     から定まるパラメターである. 

     $q$ が実数であることから伝播速度が $\beta$ と同符号でなければならない
     ことがわかる. 
     $\beta>0$ のとき $c>0$ であり, このときは常に東向きに伝播する.

     パラメターのとりうる値の一つとして, 
     $a=c=\beta=1$ のとき $q=1, k=3.9226$ の場合がある 
     (McWilliams {\it et al.}, 1981, 辻村 1993). 


 \begin{Dreference}
   \item Hoskins, B. J., Hollingsworth, A., 1972: 
     On the simplest example of the barotropic instability 
     of Rossby wave motion. 
     {\it J. Atmos. Sci.}, {\bf 30}, 150-153.
   \item McWilliams, J. C., Flierl, G. R., Larichev, V.D., Reznik, G. M., 1981:
     Numerical studies of barotropic Modons.
     {\it Dyn. Atmos. Oceans}, {\bf 5}, 219--238.
   \item 辻村 豊, 1993: モドンの数値実験, 
       気象研究ノート, 179, 2.5 章, 215--242.
 \end{Dreference}
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section*{謝辞}                           % 謝辞. このまま利用する.
\markright{謝辞}                          % ヘッダに表示

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

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

\end{document}

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