○移流方程式の解法

今回は流体シミュレーションにむけて、移流方程式（波動方程式）の数値的解法について述べます。

・移流方程式

移流方程式は、連続の式において速度一定とした場合の式に相当します。

$\displaystyle{\frac{\partial u(x,t) }{\partial t }}+c\displaystyle{\frac{\partial u(x,t)}{\partial x }}=0$

ここで$c$は移流速度（波動速度）に相当します。

初期条件を$u(x,t=0)=f(x)$とすると、この式の解は$u(x,t)=f(x-ct)$と書けます。すなわち、時刻$t$でのプロファイルは、初期$t=0$でのプロファイルの形を変えずに$ct$だけ移動したものになります。

一般的な流体計算において速度は空間的に変動し、また時間発展しますが、ここではまず、速度一定の場合について考えます。

・風上差分法

時間・空間方向ともに1次精度の解法です。$c>0$のとき、時間、空間の変数$t$,$x$に対し、以下のように差分近似します。

$\displaystyle{\frac{\partial u(x,t) }{\partial t }}\rightarrow \displaystyle{\frac{u_j^{n+1}-u_j^n }{\Delta t}}$

$c<0$のときは空間微分を以下のように差分近似します。

$\displaystyle{\frac{\partial u(x,t) }{\partial x }}\rightarrow \displaystyle{\frac{u_{j+1}^n-u_{j}^n }{\Delta x}}$

このように常に流れの上流（風上）側の差分を取ることにより、安定に計算を行う方法です。

このとき、解くべき離散化された移流方程式は、以下のようになります。

$u_j^{n+1}=u_j^n-r(u_{j}^n-u_{j-1}^n)\ \ \ \ \ {\rm if}\ c>0$

$u_j^{n+1}=u_j^n-r(u_{j+1}^n-u_{j}^n) \ \ \ \ \ {\rm if}\ c<0$

$r\equiv\displaystyle{\frac{c\Delta t}{\Delta x}}$

以下に$c>0$, 初期条件

$u(x,0)=1\ \ \ \ \ {\rm if}\ 0.1 < x < 0.3 $

$u(x,0)=0\ \ \ \ \ {\rm otherwise}$

境界条件 $u(0,t)=0$ とした場合のプログラム例を記します。

     1  program main
     2    implicit none
     3  
     4    integer,parameter:: jmax=200
     5    integer,parameter:: nmax=250
     6    integer,parameter:: nint=25
     7  
     8    double precision, dimension(0:jmax):: u
     9    double precision, dimension(1:jmax):: unew
    10    double precision delx,delt,c,r,x
    11  
    12    integer j,n
    13  
    14    character(len=10) fname
    15  
    16    delx=2.0d0/dble(jmax)
    17    delt=delx/2.0d0
    18    c=1.0d0
    19    r=c*delt/delx
    20  
    21    do j=0,jmax
    22       if(j<10 .or. j>30)then
    23          u(j)=0.0d0
    24       else
    25          u(j)=1.0d0
    26       end if
    27    end do
    28  
    29    do n=0,nmax
    30  
    31       if(mod(n,nint)==0) then
    32          write(fname,'(a,i5.5,a)') 't',n,'.dat'
    33          open(10,file=fname)
    34          do j=0,jmax
    35             write(10,'(2f12.5)') dble(j)*delx,u(j)
    36          end do
    37          close(10)
    38       end if
    39  
    40       do j=1,jmax
    41          unew(j)=u(j)*(1.0d0-r)+u(j-1)*r
    42       end do
    43  
    44       do j=1,jmax
    45          u(j)=unew(j)
    46       end do
    47  
    48    end do
    49  
    50  end program main

風上差分法を安定に解くための条件は、$r<1.0$です（クーラン条件）。$r>1.0$のとき、$c\Delta t>\Delta x$となります。情報は移流速度cで伝わりますが、$c\Delta t>\Delta x$のとき時間Δtの間に情報が格子間隔Δxを超えて伝わることになってしまい、計算が破たんします。

・ラックス・ヴェンドロフ法

時間・空間方向ともに2次精度をもつ解法です。

時間・空間に対し、テイラー展開の2次の項まで扱うことにより、以下のような差分近似を行います。

$u_j^{n+1}=u_j^n-c \displaystyle{\frac{\Delta t}{2\Delta x}}(u_{j+1}^n-u_{j-1}^n)+ \displaystyle{\frac{1}{2}}c^2\displaystyle{\biggl(\frac{\Delta t}{\Delta x}\biggr)^2}(u_{j+1}^n-2u_j^n+u_{j-1}^n)$

上述の風上差分法のプログラムにおいて、40-42行目を以下のように変更すると、ラックス・ヴェンドロフ法により移流方程式を解くことができます。

do j=1,jmax-1
   unew(j)= u(j+1)*0.5d0*r*(r-1.0d0)& 
             &+u(j)*(1.0d0-r*r)+u(j-1)*0.5d0*r*(1.0d0+r)
end do

ラックス・ヴェンドロフ法を安定に解くための条件も$r<1$です。

ラックス・ヴェンドロフ法は、流体方程式の解法として、よく使われてきた手法の１つです。

○課題○ 風上差分法とラックス・ヴェンドロフ法で安定性条件を満たす場合の計算結果を比較し、それぞれ矩形波が数値誤差のために時間とともにどのように変形するかを確認せよ。

○課題○ 初期条件がガウシアン(例えば$f(x)=\exp(-((x-0.3)/0.1)^2)$)のときの計算結果を風上差分法とラックス・ヴェンドロフ法の場合で比較し、精度の違いを確認せよ。

○課題○ 風上差分法やラックス・ヴェンドロフ法でΔt、Δxの値を変え、安定性条件を満たさない場合に計算が不安定になり、破たんすることを確認せよ。

○課題○ 空間微分を1次精度の中間差分で差分近似した場合、すなわち

$u_j^{n+1}=u_j^n-c \displaystyle{\frac{\Delta t}{2\Delta x}}(u_{j+1}^n-u_{j-1}^n)$

とした場合、この解法は任意の$r=c\Delta t/\Delta x$に対して不安定であることを確認せよ。