% 表題  3 次元 3 重周期境界領域での 3 次元非圧縮流体の強制問題の定式化
%
% 履歴  2016/07/29 竹広 真一 
%
\documentclass[12pt,a4j,notitlepage]{jarticle}
\usepackage{Dennou6}
\usepackage{psfrag}
\usepackage{ascmac}
\usepackage{times}
\usepackage[dvipdfmx]{hyperref}

\Dtitle[3 重周期境界強制乱流]
       {3 重周期境界領域での\\3 次元非圧縮流体の強制問題の定式化}
\Dauthor[竹広 真一]{竹広 真一, SPMODEL 開発グループ}

\Ddate{平成 28 年 7 月 29 日} 

\Dparskip

\pagestyle{Dheadings}

\begin{document}

\maketitle

\section{はじめに}

  直線直交座標 $(x,y,z)$ で張られた
  $0 \le x \le 2\pi, 0 \le y \le 2\pi, 0 \le z \le 2\pi$ の
  立方体の領域に詰めこまれている密度一定の非圧縮流体の運動を考える. 
  境界は $x,y,z$ それぞれの方向に対して周期的である. 


\section{支配方程式}

  支配方程式は, 非圧縮流体の拡張されたナビエストークス方程式である. 
  %
  \begin{equation}
    \Deqlab{ナビエストークス}
    \DP{\Dvect{u}}{t} + (\Dvect{u}\cdot\Dgrad)\Dvect{u} 
    = - \frac{1}{\rho}\Dgrad p + (-1)^{p+1}\nu_{2p}\nabla^{2p}\Dvect{u}, 
    \quad
    \Ddiv\Dvect{u} = 0. 
  \end{equation}
  %
  ここで $\Dvect{u}$ は速度, $p$ は圧力, $\rho$ は密度(一定), 
  $\nu_{2p}$ は高階粘性の係数であり $p=1$ で通常の粘性および
  ナビエストークス方程式となる. 

  この方程式の回転をとると, 渦度方程式が得られる. 
  %
  \begin{equation}
    \Deqlab{渦度方程式}
    \DP{\Dvect{\omega}}{t} 
    = \Drot(\Dvect{u}\times\Dvect{\omega})
    + (-1)^{p+1}\nu_{2p}\nabla^{2p}\Dvect{\omega}.
  \end{equation}
  %
  $\Ddiv\Dvect{\omega}=0$ であるからこの方程式の 2 成分が独立な式である. 
  渦度方程式の適当な 2 成分の時間発展を計算し, 
  残りの渦度成分は渦度非発散条件から求められる. 
  渦度全ての成分が求まれば, 渦度を再度回転をとった式
  %
  \begin{equation}
    \Drot\Dvect{\omega} = - \Dlapla\Dvect{u}.
  \end{equation}
  %
  から速度場が求められる. 


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

  \subsection{フーリエ展開}

    スペクトル法による数値計算を行うために, 
    各物理量をフーリエ変換する. すなわち, 
    %
    \begin{eqnarray}
      \omega_i(x,y,z,t)
          &=&\sum_{n=-N}^N\sum_{m=-M}^M\sum_{m=-L}^L
              \tilde{\omega}_i(l,m,n,t) e^{i(lx+my+nz)}, \\
      \tilde{\omega}_i(l,m,n,t) 
          &=& \frac{1}{(2\pi)^3}\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi}
                 \omega_i(x,y,z,t) e^{-i(lx+my+nz)} dx dy dz.
    \end{eqnarray}
    %
    ここで, $\omega_i$ は渦度の $i$ 成分であり, 実数である. 
    したがってそのフーリエ展開係数 $\tilde{\omega}_i$ は
    %
    \begin{equation}
      \tilde{\omega}_i(-l,-m,-n) = \tilde{\omega}_i^*(l,m,n)
    \end{equation}
    %
    の関係を満たしている. ただし ${}^*$ は複素共役を表す. 

    逆変換の関係 $\omega_i$ を\Deqref{渦度方程式}に代入し, 
    $ e^{-i(lx+my+nz)}$ をかけて全領域で積分することにより, 
    直交関係から以下の式が得られる. 
    %
    \begin{equation}
      \Deqlab{スペクトル渦度方程式}
      \DD{}{t}\tilde{\omega}_i(l,m,n,t)
      = \widetilde{[\Drot(\Dvect{u}\times\Dvect{\omega})]}_i(l,m,n,t)
      + \nu^*(l,m,n,p)\tilde{\omega}_i(l,m,n,t).
    \end{equation}
    %
    ここで, 波数毎の超粘性係数 $\nu^*(n,p)$ を
    %
    \begin{equation}
        \nu^*(l,m,n,p)\equiv (-1)^p\nu_{2p}(l^2+m^2+n^2)^p
    \end{equation}
    %
    と定義した($\nu^*(n,p)$ は 0 以上であることに注意). 


  \subsection{時間積分}

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

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



\section{実験設定}

  \subsection{強制設定}

    小スケールにおける強制は, スペクトル空間 $(l,m,n)$ 上における
    球殻領域 $ K_{min} \le \sqrt{l^2+m^2+n^2} \le K_{max}$ 内に
    存在するスペクトル成分からなる運動エネルギーを 1 に固定する. 
    ここで $K_{min},K_{max}$ は強制を入力する最小点最大全波数である. 

  \subsection{初期条件}

    初期条件は静止状態である. 

  \subsection{パラメター}

    未定. 


\begin{Dreference}
  \item SPMODEL ライブラリレファレンスマニュアル「eee\_mpi\_module.f90」,\\
    \url{http://www.gfd-dennou.org/library/spmodel/spml_current/html/index.html}
  \item ISPACK ver.1 解説文書「P3PACK 使用の手引」,\\
    \url{http://www.gfd-dennou.org/library/ispack/}
\end{Dreference}


\end{document}

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