.. sectionauthor::   清水康行 <yasu@i-ric.org>

.. _2d_advection_eq:

=====================================
2次元移流方程式
=====================================

基礎式
================

2次元移流方程式を次式で表します.

.. math::
    :label: 2d_advection

    {{\partial f}\over{\partial t}} +
    u{{\partial f}\over{\partial x}}+
    v{{\partial f}\over{\partial y}} = 0

ここで, :math:`t` は時間, :math:`x, y` は直交座標軸, :math:`f` は移流される変数, 
:math:`u,v` は :math:`x, y` 方向の移流速度です.

1次元の移流方程式と同様，この方程式の理論解は :math:`f` はその分布形を保ったまま，
:math:`x` 方向には :math:`u` の速度で, :math:`ｙ` 方向には :math:`v` の速度で移動
するという結果となります． 

.. _2dkazakami:

風上差分
=================

移流項の風上差分式は下記となります．

.. math::
    :label: 2d_ad_1

    u>0 \mbox{のとき}; \; \; 
    u\cfrac{\partial f}{\partial x}=u\cfrac
    {f(i,j)-f(i-1,j)}{\Delta x}    \\
    u<0 \mbox{のとき}; \; \; 
    u\cfrac{\partial f}{\partial x}=u\cfrac
    {f(i+1,j)-f(i,j)}{\Delta x}    

.. math::
    :label: 2d_ad_2

    v>0 \mbox{のとき}; \; \; 
    v\cfrac{\partial f}{\partial x}=v\cfrac
    {f(i,j)-f(i,j-1)}{\Delta y}    \\
    v<0 \mbox{のとき}; \; \; 
    v\cfrac{\partial f}{\partial x}=v\cfrac
    {f(i,j+1)-f(i,j)}{\Delta y}    

ここで, :math:`\Delta x, \Delta y, \Delta t` はそれぞれ
:math:`x, y, t` 方向の計算刻み幅で, :math:`(i,j)` はそれぞれ
:math:`(x,y)` 方向の格子番号です．
したがって, 例えば, :math:`u>0, \; v>0` の場合, :eq:`2d_advection` 式の
差分式は下記のようになります．

.. math::
    :label: 2d_diff_1

    \cfrac{f_n(i,j)-f(i,j)}{\Delta t}+
    u\cfrac{f(i,j)-f(i-1,j)}{\Delta x}+
    v\cfrac{f(i,j)-f(i,j-1)}{\Delta y}=0

ただし, :math:`f_n(i,j)` は :math:`t=t+\Delta t` における :math:`f(i,j)` の値です．
:eq:`2d_diff_1` 式より :math:`f_n(i,j)` は次式を:math:`\Delta t` 毎に
繰り返し計算することにより求めることが出来ます. なお，この場合は :math:`x` 方向にも
:math:`y` 方向にも後進差分ということになります．

.. math::
    :label: 2d_diff_2

    f_n(i,j)=f(i,j)-\Delta t \left\{
    u\cfrac{f(i,j)-f(i-1,j)}{\Delta x}+
    v\cfrac{f(i,j)-f(i,j-1)}{\Delta y}\right\}

後進差分による計算例
---------------------

:math:`u>0, \; v>0` の場合の計算例を示します．:math:`f` の初期分布は :numref:`3dcone` 図に示すような
円錐形の分布を与えます．


        .. _3dcone:

.. figure:: images/04/3dcone.png
    :width: 70%

    : 2次元移流計算の初期条件

ちなみに，円錐形分布は次式で与えます．

.. math::
    :label: cone_profile

    
    r<R \; \mbox{のとき;} &  f(x,y)=\cfrac{R-r}{R} f_p \\
    r>R \; \mbox{のとき;} &  f(x,y)=0

ただし, 

.. math::
    :label: teimen

    r=\sqrt{(X_p-x)^2+(Y_p-y)^2}
    
また，:math:`X_p, Y_p` ピークの位置, :math:`f_p` はピークの値, :math:`R` は底面の半径です．

Pythonコードと計算例
^^^^^^^^^^^^^^^^^^^^^^^^^^

:eq:`2d_diff_2` を使ったPythonコードを
https://github.com/YasuShimizu/Python-Code/blob/main/2d_backward_p3d.py
に示します．また，このコードを用いて計算された計算結果を下図に示します．
なお，この例では :math:`u=v=0.5` (m/s), :math:`f_p=0.5, R=20` (m)としてます．

        .. _3dcone_backward:

.. figure:: images/04/3dcone_backward.gif
    :width: 95%

    : 2次元移流方程式の計算結果(後進差分)

１次元の時と同様，この計算結果は数値拡散の影響でピークの値が徐々に下がりながら，円錐の形も
潰れながら移動していおり，けっして正確な計算結果とは言えません．

移流方向が変化する場合の風上差分法
-------------------------------------

:eq:`2d_ad_1` 式, :eq:`2d_ad_2` 式の移流項を一つの式で表した，:math:`f` の更新
式は下記となります．

.. Math::

    f_n(i,j)=f(i,j)
    
.. Math::

    -\cfrac{\Delta t}{2 \Delta x}
    \left[(u+|u|)\left\{f(i,j)-f(i-1,j)\right\}+(u-|u|)\left\{f(i+1,j)-f(i,j)\right\} \right. 

.. Math::
    :label: kazakami2d  

    +\left.(v+|v|)\left\{f(i,j)-f(i,j-1)\right\}+(v-|v|)\left\{f(i,j+1)-f(i,j)\right\} \right]

:eq:`kazakami2d` 式は :math:`u, v` の大きさや符号に関わらず適用可能な風上差分法です．

では，この2次元風上差分法を :math:`u, v` が変化する状況で試してみましょう． いま，:math:`u,v` が次式で
振動する状況を想定します．

.. Math::
    :label: circular_uv

    u = & -r_t \omega \sin \theta \\
    v = &  r_t \omega \cos \theta 

この式で表される運動は円運動で, :math:`r_t` は円運動の曲率半径， :math:`\omega` は角速度, 
:math:`\theta` は位相角です． もし, :numref:`3dcone` 図に示した円錐形のピークが，:eq:`circular_uv` 
式に従って運動するとした場合の結果は下図( :numref:`3d_cone_exact_solution` )のようになります．
なお, この例では, :math:`T` を円運動の1周の周期として :math:`T=` 100sec, 
:math:`r_t` = 30m, `f_p` =0.5 とした計算です．
また，:math:`T` と :math:`\omega` の関係は次式で表されます．

.. Math::
    :label: shuuuki

    \omega=\cfrac{2\pi}{T}

.. _3d_cone_exact_solution:

.. figure:: images/04/CircularMovementCone.gif
    :width: 110%

    : 円錐形状の円運動

:numref:`3d_cone_exact_solution` は単に円錐形の分布を :eq:`circular_uv` 式の速度で
動かしただけで，:eq:`2d_advection` 式の差分計算結果ではありませんが，ある意味これが
これから行う数値計算の目標となる厳密解です． なお参考までに，:numref:`3d_cone_exact_solution` 図のアニメーションを
作成するPythonコードを下記に示します．

https://github.com/YasuShimizu/Python-Code/blob/main/3dcone_circle_m3d.py

前置きが長くなってしまいましたが，本節の目標である2次元移流方程式の風上差分解法を
やってみましょう．

Pythonコードと計算例
^^^^^^^^^^^^^^^^^^^^^^^^^^^

:eq:`2d_advection` 式を用いたPythonコードを下記に示します．

https://github.com/YasuShimizu/Python-Code/blob/main/2d_upwind_circle.py

このコードを使って，:eq:`cone_profile` 式の初期条件で, :math:`u,v` を :eq:`circular_uv` 式
で変化させながら計算した結果は下図となります．

.. _3d_upwind_circular:

.. figure:: images/04/2d_upwind_circle.gif

    : 2次元風上差分法による円錐形分布の円運動の計算

予想どおりではありますが，風上差分法で初期形状が徐々に拡散し，ピークの値も減衰してしまってます．
計算結果は安定していますが，1次元の時と同様，風上差分法では正確な計算結果が得られないことが示されました．

演習問題(4)
^^^^^^^^^^^^^^

.. list-table:: 

    * - 初期条件を円錐形分布とする変数 :math:`f` が四角形の形状の辺に沿って連続的に移動する場合の :math:`f` の時間変化を計算して下さい．
    
高精度の差分法
=================

ここで，1次元の場合と同様に2次元の高精度差分法として2次元版のCIP法について述べます．

2次元CIP法による移流項の計算( :math:`u<0, v<0` の場合)
--------------------------------------------------------------

説明は簡単のためにまず :math:`u<0` , :math:`v<0` の場合について行います.

    .. _2duv:

.. figure:: images/04/2duv.png
   :width: 90%

   :  計算平面上での :math:`f` の推定

:numref:`2duv` 示すように :math:`x, y` 方向の計算点番号を :math:`i,j` として
時刻 :math:`t` で :math:`i, i+1, j, j+1` の範囲にある点(図中の白丸)が
:math:`\Delta t` 時間で図中の黒丸の点に移動すると考えます.

1次元の計算で説明したのとまったく同じようにこの問題は :numref:`2duv` の
黒丸の点における :math:`f` の分布形をいかに推定するかという問題に帰着します.
即ち, 

.. math::
    :label: 2d_1

    X=-u \Delta t, \; \;  Y=-v \Delta t

のときこの分布形 :math:`F(X,Y)` が推定できれば

.. math::
    :label: 2d_2

    f(i,j,t+\Delta t) = F(X, Y)

により :math:`f` の時間更新ができます. ここで, :math:`f` は　:math:`(i,j)` 格子の離散値, 
:math:`F` は連続関数であることに注意して下さい．

:math:`F` の関数形は1次元問題のときは3次曲線としましたが, 
ここでは3次曲面とします. :math:`F(X,Y)` は :math:`X=0, Y=0` で

.. math::
    :label: 2d_3

    F(X,Y)=f(i,j), \; \; \; F_x(X,Y)=f_x(i,j), \; \; \;F_y(X,Y) = f_y(i,j) 

(ただし, :math:`F_x(X,Y)=\displaystyle{{\partial F}\over{\partial x}}` , 
:math:`F_y(X,Y)=\displaystyle{{\partial F}\over{\partial y}} (X,Y)` ,
:math:`f_x(i,j)=\displaystyle{{\partial f}\over{\partial x}} (i,j)` ,
:math:`f_y(i,j)=\displaystyle{{\partial f}\over{\partial y}} (i,j)` です．)

を満たす必要があるので, :math:`F(X,Y)` を次式のように置きます.

.. math::
    :label: 2d_5

    F(X,Y) =
    \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + f_x(i,j) \right] X +  \\
    \left[ (b_1 Y + d_1 X + f_1) Y + f_y(i,j) \right] Y + f(i,j)


この式は :eq:`2d_3` 式を自動的に満たしています. 他の条件としては,

.. math::
    :label: 2d_6

    F(0,-\Delta y)=f(i,j+1), \; \; \;
    F(-\Delta x, 0)=f(i+1,j) \\
    F(-\Delta x, -\Delta y) = f(i+1,j+1), \; \; \; 
    F_x(0, -\Delta y)=f_x(i,j+1) \\
    F_x(-\Delta x, 0)=f_x(i+1,j), \; \; \;
    F_x(-\Delta x, -\Delta y) = f_x(i+1,j+1) \\
    F_y(0, -\Delta y)=f_y(i,j+1), \; \; \;
    F_y(-\Delta x, 0)=f_y(i+1,j)

などを満たす必要があります. これらの関係を用いることにより, :eq:`2d_5` 式中の係数は以下の
ようになります.

.. math::

    a_1 = {{
    \left[ f_x(i+1,j)+f_x(i,j) \right] \Delta x
    +2 \left[ f(i,j)-f(i+1,j) \right]}\over{\Delta x^3}}

.. math::

    b_1 = {{
    \left[ f_y(i,j+1)+f_y(i,j) \right] \Delta y
    +2 \left[ f(i,j)-f(i,j+1) \right]}\over{\Delta y^3}}

.. math::

    c_1 = {{-\left\{
    f(i,j)-f(i,j+1)-f(i+1,j)+f(i+1,j+1)\right\}-\left[f_x(i,j+1)-f_x(i,j) \right] \Delta x
    }\over{\Delta x^2 \Delta y}}

.. math::

    d_1 = {{-\left\{f(i,j)-f(i,j+1)-f(i+1,j)+f(i+1,j+1)\right\}-
    \left[ f_y(i+1,j)-f_y(i,j) \right] \Delta y}\over{\Delta x \Delta y^2}}

.. math::

    e_1 = {{
    3 \left[ f(i+1,j)-f(i,j) \right] +
    \left[ f_x(i+1,j)+2 f_x(i,j) \right] \Delta x
    }\over{\Delta x^2}}

.. math::

    f_1 = {{
    3 \left[ f(i,j+1)-f(i,j) \right] +
    \left[ f_y(i,j+1)+2 f_y(i,j) \right] \Delta y
    }\over{\Delta y^2}}

.. math::
    :label: coefs

    g_1 =
    {{-\left\{f_y(i+1,j)-f_y(i,j)\right\}+c_1 \Delta x^2
    }\over{\Delta x}}

演習問題(5)
^^^^^^^^^^^^^^

.. list-table:: 

    * - :eq:`coefs` 式の誘導を自分で行って, :math:`a_1 \sim g_1` の表現が正しいかどうか確かめて下さい.

2次元CIP法による移流項の計算(一般化)
---------------------------------------

ここまでの説明は :math:`u<0` , :math:`v<0` の場合でしたが, 実際には :math:`u, v` の符号の組み合わせに
よって8通りの組み合わせがありますが, これをすべて定式化してプログラミング
するのは非常に煩雑です. そこで1次元でやったように符号関数( :math:`\mbox{sign}` )関数を導入して
表現することにします.


.. math::
    :label: advect2d_1

    i_s = \mbox{sign}(u), \; \; \;  & j_s = \mbox{sign}(v) \\
    i_m = i-i_s, \; \; \;  & j_m = j-j_s 


格子間の変数の分布形状 :math:`F(X,Y)` は :math:`X, Y` を :numref:`2duv_1` のように定義し, 
:eq:`fxnew` 式を2次元に拡張することにより, :eq:`advect2d_2` のように定義します．

    .. _2duv_1:

.. figure:: images/04/2duv_1.png
   :width: 80%

   :  計算平面上での :math:`X, Y` の定義 

.. math::
    :label: advect2d_2

    F(X,Y) = \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + f_x(i,j) \right] X + \\
    \left[ (b_1 Y + d_1 X + f_1) Y + f_y(i,j) \right] Y + f(i,j)


ただし, 係数 :math:`a_1 \sim g_1` は :math:`i_s, j_s, i_m, j_m` 
を用いることにより, :math:`u, v` の符号に関係なく以下のように表現されます.

.. math::

    a_1 = {{
    i_s \left[ f_x(i_m,j)+f_x(i,j) \right] \Delta x
    -2 \left[ f(i,j)-f(i_m,j) \right]}\over{i_s \Delta x^3}}

.. math::

    b_1 = {{
    j_s \left[ f_y(i,j_m)+f_y(i,j) \right] \Delta y
    -2 \left[ f(i,j)-f(i,j_m) \right]}\over{j_s \Delta y^3}}

.. math::

    c_1 = {{
    -\left\{f(i,j)-f(i,j_m)-f(i_m,j)+f(i_m,j_m)\right\}-i_s \left[f_x(i,j_m)-f_x(i,j) \right] \Delta x
    }\over{j_s \Delta x^2 \Delta y}}

.. math::
 
    d_1 = {{
    -\left\{f(i,j)-f(i,j_m)-f(i_m,j)+f(i_m,j_m)\right\}-
    j_s \left[ f_y(i_m,j)-f_y(i,j) \right] \Delta y}\over{i_s \Delta x \Delta y^2}}

.. math::

    e_1 = {{
    3 \left[ f(i_m,j)-f(i,j) \right] +
    i_s \left[ f_x(i_m,j)+2 f_x(i,j) \right] \Delta x
    }\over{\Delta x^2}}

.. math::

    f_1 = {{
    3 \left[ f(i,j_m)-f(i,j) \right] +
    j_s \left[ f_y(i,j_m)+2 f_y(i,j) \right] \Delta y
    }\over{\Delta y^2}}

.. math::
    :label: coefs_1
    
    g_1 =
    {{-\left\{f_y(i_m,j)-f_y(i,j)\right\}+c_1 \Delta x^2
    }\over{i_s \Delta x}}

演習問題(6)
^^^^^^^^^^^^^^

.. list-table:: 

    * - :eq:`coefs_1` 式の誘導を自分で行って, :math:`a_1 \sim g_1` の表現が正しいかどうか確かめて下さい.

微分量の更新
-----------------

CIP法では変数 :math:`f` の他にもその微分量 :math:`f_x, f_y` も時々刻々更新して行く必要
があります. そこで :eq:`2d_advection` 式の両辺をまず :math:`x` で偏微分します.

.. math::
    :label: bibun2_1

    {{\partial^2 f}\over{\partial t \partial x}}
    + {{\partial u}\over{\partial x}} {{\partial f}\over{\partial x}}
    + u {{\partial^2 f}\over{\partial x^2}}
    + {{\partial v}\over{\partial x}} {{\partial f}\over{\partial y}}
    + v {{\partial^2 f}\over{\partial x \partial x}} =0


:math:`\displaystyle{{{\partial f}\over{\partial x}}}` を :math:`f_x` とすると,

.. math::
    :label: bibun2_2

    {{\partial f_x}\over{\partial t}}
    + u{{\partial f_x}\over{\partial x}}
    + v{{\partial f_x}\over{\partial y}}
    = -\left( f_x {{\partial u}\over{\partial x}} +
    f_y {{\partial v}\over{\partial x}} \right)

となりまして, 左辺は :`f_x` に関する2次元の移流方程式になっていることが分かります.
従って, 前述の2次元のCIP法をそのまま適用すれば良いことに
なります. また, 右辺に関してはもし :math:`u, v` が定数であればゼロとなりますが,
変数の場合には, 1次元のCIP法で説明した[ :ref:`label1` ]を適用して
分離解法で計算することになります. なお右辺の差分形式は中央差分で行けます.

同様に, :eq:`2d_advection` 式を :math:`y` で偏微分して
:math:`\displaystyle{{{\partial f}\over{\partial y}}}` を :math:`f_y` とすると,

.. math::
    :label: bibun2_3

    {{\partial f_y}\over{\partial t}}
    + u{{\partial f_y}\over{\partial x}}
    + v{{\partial f_y}\over{\partial y}}
    = -\left( f_x {{\partial u}\over{\partial y}} +
    f_y {{\partial v}\over{\partial y}} \right)

が得られ :eq:`bibun2_2` 式と同様な形となり, 左辺の移流項の部分はCIP法で, 右辺の
Source項の部分は中央差分で計算可能となります.

Pythonコードと計算例
-----------------------------

以上で述べた2次元CIP法のPhthonコードを下記に示します．

https://github.com/YasuShimizu/2d_advection_cip

        .. _3dcone_cip:

このコードによる計算結果を下記に示します．計算条件は前記の :ref:`2dkazakami` と同じく円運動をする円錐形分布です．

.. figure:: images/04/2d_cip_circle.gif
    :width: 95%

    : 2次元移流方程式の計算結果(CIP法)

計算結果によれば，円錐形の形状はほとんど形を保ったまま円運動を続け，2次元の場合もCIP法は移流方程式の数値計算法として
高精度の計算結果が得られることが分かります．