移流方程式

移流方程式(advection equation)

はじめに

移流現象とは，物質や運動量などの物理量が流れに従って運ばれる現象である（運動量や運動量に関係する物理量の場合は非線形となる．前に説明したBurgers方程式の左辺は非線形な移流現象を表している）．

非常に単純な方程式であるが，流体の運動を支配するNavier-Stokes方程式を数値的に解く際に最も重要な（線形であったとしても面倒な）部分となる．ここでは，流れの速度が一定の $u$ で運ばれる物質の濃度が $C$ の場合（並進方程式）を考える．一次元の場合の模式図を熱拡散の場合と同様に以下に示す．

図(a)の概念図から方程式を誘導する（差分法の解説であるが，基礎式の誘導はコントロールボリュームを用いた有限体積法的である．支配方程式が保存則がベースであるので．）．まず対象となる微小領域の濃度変化を考える．単位時間内にこの領域内に流入するフラックスは $u C (x) \cdot Δ y \cdot 1$ であり，流出するフラックスは $u C (x + Δ x) \cdot Δ y \cdot 1$ ，その結果として微小領域内濃度の時間変化は $\frac{\partial C}{\partial t} \cdot Δ x \cdot Δ y \cdot 1$ となる．テイラー展開を行い， $Δ x \to 0$ の極限を取ると以下の式が得られる．

\begin{matrix} (1) & \frac{\partial C}{\partial t} = - \frac{\partial u C}{\partial x} \end{matrix}

この式の形式を保存形と呼ぶ．これは境界面での物理量の収支を計算しているので，物質の総量が保存される（ある計算セルから出て行った量は必ず隣の計算セルに入る．式自身が「物質がなくならないこと」（連続式）を保証している）．一方，連続式が満たされる $\frac{\partial u}{\partial x} = 0$ と仮定して式を変形することもできる．

\begin{matrix} (2) & \frac{\partial C}{\partial t} = - \frac{\partial u C}{\partial x} = - u \frac{\partial C}{\partial x} - C \frac{\partial u}{\partial x} = - u \frac{\partial C}{\partial x} \end{matrix}

この式は非保存型と呼ばれ，式自身は連続式を保証していない．今回の場合は移流速度 $u$ が一定であるため，保存形と非保存形は一致する．また，多くの場合，全て左辺に移項させて以下のように表すことが多い．

\begin{matrix} (3) & \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} = 0 \end{matrix}

FTCS

熱拡散の時にも述べたが，方程式( $2$ )を近似する様々な差分が考えられるが，まずは熱拡散と同様にFTCS法を用いた離散化を試みる．差分法の所で学んだ時間に関する前進差分，空間に関する２階の中心差分を適用すると，以下の式が得られる．

\begin{matrix} (4) & \frac{C_{i}^{n + 1} - C_{i}^{n}}{Δ t} = - u \frac{C_{i + 1}^{n} - C_{i - 1}^{n}}{2 Δ x} \end{matrix}

整理すると

\begin{matrix} (5) & C_{i}^{n + 1} = C_{i}^{n} - \frac{u Δ t}{Δ x} \frac{C_{i + 1}^{n} - C_{i - 1}^{n}}{2} \end{matrix}

次に，より物理的に理解しやすい方法（有限体積法的な離散化方法）を説明する．はじめに示した図(b)に見られるような離散化を行う(時間の添字がついていない変数は全て既知の時間ステップ $n$ の値とする) ．空間の中心差分は計算セル（微小領域）の境界面における濃度を両側の計算点における値の平均値で表すことに対応する．

\begin{matrix} (6) & \frac{C_{i}^{n + 1} - C_{i}^{n}}{Δ t} \cdot Δ x \cdot Δ y \cdot 1 = u \frac{C_{i - 1}^{n} + C_{i}^{n}}{2} \cdot Δ y \cdot 1 - u \frac{C_{i}^{n} + C_{i + 1}^{n}}{2} \cdot Δ y \cdot 1 \end{matrix}

整理整頓すると式( $5$ )と同じ式が得られる．

\begin{matrix} (7) & C_{i}^{n + 1} = C_{i}^{n} - \frac{u Δ t}{Δ x} \frac{C_{i + 1}^{n} - C_{i - 1}^{n}}{2} \end{matrix}

ここで，プログラムの作成・実行前にVon Neumannの安定解析を行う．

\begin{matrix} (8) & | \frac{r^{n + 1}}{r^{n}} | = | 1 - i \frac{u Δ t}{Δ x} \sin (k Δ x) | \end{matrix}

複素数であるので，絶対値を計算すると常に1よりも大きくどのような条件でも安定にならないことがわかる( $u = 0$ でそのような条件になるが，それでは移流しない事になる)．このようなスキームは無条件不安定なスキームと呼ばれている．

\begin{matrix} (9) & | \frac{r^{n + 1}}{r^{n}} | = \sqrt{1 + {(\frac{u Δ t}{Δ x})}^{2} \sin^{2} (k Δ x)} > 1 \end{matrix}

実際にプログラムを作成して計算してみる．計算は，領域 $x = 0 \sim 1$ mである物質（スカラー量）の濃度 $C$ が初期（ $t = 0$ ）に $x = 0.2 \sim 0.5$ mで $C = 1$ ，それ以外の区間で $C = 0$ である分布（矩形波）が $x$ の正の方向に速度 $u = 1.0$ m/sで移流する現象を扱う（この移流問題は，単純なようで実は難しい）．

厳密解は，分布の形を変えず速度1.0m/sで $x$ の正の方向へ移動するものとなる．時間ステップを $Δ t = 0.01$ s，空間方向の刻み $Δ x$ =0.05m（領域を20分割）とする．また，この問題では境界の影響を受けない部分に着目した計算を行うので，合理的な境界条件は設定しない（というかゼロで固定）．

 program advection
  implicit none
  real(8) :: dx=0.05d0, dt=0.01d0, u=1.d0, s, x
  real(8), dimension(20) :: c=0.d0, c_old
  integer :: i, n, i_end=20, n_end=30, write_interval=2
  character(128) fname
  ! parameter
  s = u*dt/dx
  ! writing initial condition
  write (fname, '("advection_result", i3.3, ".txt")') 0 
  open(88, file=fname, status='replace')
  x = -0.5*dx
  do i=1, i_end
    x=x+dx
    if(x>0.2 .and. x<0.5) then
      c(i)=1.0
    end if
    write(88,*) x, c(i)
  end do
  close(88)
  c_old = c
  ! main loop
  do n=1, n_end
    do i=2, i_end-1 ! exclude boundaries
      c(i)=c_old(i)-0.5*s*(c_old(i+1)-c_old(i-1))
    end do
    ! writing result
    if(mod(n,write_interval).eq.0) then
      write (fname, '("advection_result", i3.3, ".txt")') n 
      open(99, file=fname, status='replace')
      x = -0.5*dx
      do i=1, i_end
        x=x+dx
        write(99,*) x, c(i)
      end do
      close(99)
    endif
    c_old = c ! stock data
  end do  
end program advection

結果のアニメーションを作成するGnuplotプログラムを以下に示す．

 reset
dt = 0.01
set term gif animate optimize size 640,480
set output 'advection_FTCS.gif'
set title "Advection by FTCS"
set xlabel "x (m)"
set ylabel "C"
set yrange [-0.5:1.5]
do for [j = 0:30:2 ] {
plot sprintf("./advection_result%03d.txt", j) with linespoints lw 3 title sprintf("time = %f s", j*dt)}

下のアニメーションに示すように，分布は本来あるべき初期の形で移動（並進）することなく，徐々に発散して行くことがわかる．

風上差分

FTCSスキームでは計算が発散することがわかった．そこで，他の差分スキームの適用を試みる．時間差分は前進差分のままに固定して，空間の離散化スキームを変更してみよう．高次スキームは後ほど議論することにして，まずは，低次の前進差分と後退差分を検証することにする．前進差分による式( $2$ )の離散化式，FTFSスキームは以下の通りとなる．

\begin{matrix} (10) & C_{i}^{n + 1} = C_{i}^{n} - \frac{u Δ t}{Δ x} (C_{i + 1}^{n} - C_{i}^{n}) \end{matrix}

このFTFSスキームによる結果アニメーションを以下に示す．アニメーションから確認できる様にFTCSスキームよりひどく，全く並進している様子がない．これまでと同様なVon Neumannの安定条件を計算すると，

\begin{matrix} (11) & | \frac{r^{n + 1}}{r^{n}} | = \sqrt{1 + 2 c (1 + c) {1 - \cos (k Δ x)}} \end{matrix}

ここで $c = \frac{u Δ t}{Δ x}$ であり，クーラン数と呼ばれている．この式中の $1 - \cos (k Δ x) \geq 0$ であるため， $| \frac{r^{n + 1}}{r^{n}} | \leq 1$ を満足するためには， $2 c (1 + c) \leq 0$ である必要がある．このことから， $- 1 \leq c \leq 0$ でなければならず，今回の問題の様に速度が正の場合には安定にはなりえず，不安定化することがわかる．

今度は空間スキームに後退差分を利用するFTBSスキームを考え，離散化式は以下の様になる．

\begin{matrix} (12) & C_{i}^{n + 1} = C_{i}^{n} - \frac{u Δ t}{Δ x} (C_{i}^{n} - C_{i - 1}^{n}) \end{matrix}

下に示すアニメーションから，FTBSでは形は変化してしまっている（矩形から拡散した様ななだらかな形状への変化）が，並進移動に関して最もよく表現できていることがわかる．FTBSスキームの安定条件を考える．

\begin{matrix} (13) & | \frac{r^{n + 1}}{r^{n}} | = \sqrt{1 - 2 c (1 - c) {1 - \cos (k Δ x)}} \end{matrix}

この結果から，安定であるためには $2 c (1 - c) \geq 0$ でなければならず，安定のためには， $0 \leq c \leq 1$ である必要がある． $c$ が正とすれば， $c = \frac{u Δ t}{Δ x} \leq 1$ という条件になり，この条件はCFL条件と呼ばれ，多くの移流スキームの安定性を評価する条件となっている．この計算例では $c = 0.2$ であったため，安定に計算ができている．

物理的な意味でなぜFTBSスキームが安定なのかを考えてみる．今回考えた $u > 0$ の移流計算は，全ての情報が $x > 0$ の方向へ流れていくため，境界面の物理量を評価する際に上流側の影響が相対的に大きいと思われるので，上流側の計算点における値に高い重みが付けられるスキームが式の意味に沿ったものになることが予想される．FTBTスキームは，境界面の値を評価する際に上流側のみを利用しており，この意に沿ったスキームとなっている．この様に上流側に重みを大きくした差分スキームは上流差分，もしくは風上差分と呼ばれている．

FTBSスキームは，安定であったが明らかに精度が低下している．差分法の概要説明で述べた様に，時間に関しても空間に関しても１次精度スキームとなっている．より精度の高い差分スキームとしてLax-Wendroff法が提案されている．風上差分は２階微分の項が誤差の主要項になっており，２階微分は拡散的に振る舞うため分布形がなだらかになると考えられる．そこで，時間差分で２次の項まで再現し，空間を中心差分で展開する方法を考える．まず， $c_{i}^{n + 1}$ に関してテイラー展開を利用する．

\begin{matrix} (14) & C_{i}^{n + 1} = C (x_{i}, t + Δ t) = C (x_{i}, t) + \frac{\partial C (x_{i}, t)}{\partial t} Δ t + \frac{1}{2} \frac{\partial^{2} C (x_{i}, t)}{\partial t^{2}} Δ t^{2} + O (Δ t^{3}) \end{matrix}

ここで，支配方程式から $\frac{\partial C}{\partial t} = - u \frac{\partial C}{\partial x}$ ， $\frac{\partial^{2} C}{\partial t^{2}} = u^{2} \frac{\partial^{2} C}{\partial x^{2}}$ が得られるため，上式に代入すると以下の式が得られる．

\begin{array}{rcl} C (x_{i}, t + Δ t) & = & C (x_{i}, t) - u \frac{\partial C (x_{i}, t)}{\partial x} Δ t + \frac{u^{2}}{2} \frac{\partial^{2} C (x_{i}, t)}{\partial x^{2}} Δ t^{2} + O (Δ t^{3}) \\ (15) & \sim & C_{i}^{n} - (\frac{u Δ t}{Δ x}) \frac{C_{i + 1}^{n} - C_{i - 1}^{n}}{2} + {(\frac{u Δ t}{Δ x})}^{2} \frac{C_{i + 1}^{n} - 2 C_{i}^{n} + C_{i - 1}^{n}}{2} \end{array}

Lax-Wendroff法で計算した結果をアニメーションで示すと以下の様になる．結果を見ると，FTCSスキームの様に振幅が徐々に増大するわけではなく，さらに風上差分の様に振幅が急激に減衰することもない．また，安定条件は風上差分と同じ $c \leq 1$ となっている．しかしながら，Lax-Wendroff法では数値振動（高周波の波）は存在し，初期の最大値(1)や最小値(0)の範囲を超えた値が出現（オーバーシュート・アンダーシュート）していることがわかる．これらを抑制しながら高次精度を実現するスキームが提案されているが，この授業の範囲内ではないため，説明は省略する．

	program advection
	implicit none
	real(8) :: dx=0.05d0, dt=0.01d0, u=1.d0, s, x
	real(8), dimension(20) :: c=0.d0, c_old
	integer :: i, n, i_end=20, n_end=30, write_interval=2
	character(128) fname
	! parameter
	s = u*dt/dx
	! writing initial condition
	write (fname, '("advection_result", i3.3, ".txt")') 0
	open(88, file=fname, status='replace')
	x = -0.5*dx
	do i=1, i_end
	x=x+dx
	if(x>0.2 .and. x<0.5) then
	c(i)=1.0
	end if
	write(88,*) x, c(i)
	end do
	close(88)
	c_old = c
	! main loop
	do n=1, n_end
	do i=2, i_end-1 ! exclude boundaries
	c(i)=c_old(i)-0.5s(c_old(i+1)-c_old(i-1))
	end do
	! writing result
	if(mod(n,write_interval).eq.0) then
	write (fname, '("advection_result", i3.3, ".txt")') n
	open(99, file=fname, status='replace')
	x = -0.5*dx
	do i=1, i_end
	x=x+dx
	write(99,*) x, c(i)
	end do
	close(99)
	endif
	c_old = c ! stock data
	end do
	end program advection

	reset
	dt = 0.01
	set term gif animate optimize size 640,480
	set output 'advection_FTCS.gif'
	set title "Advection by FTCS"
	set xlabel "x (m)"
	set ylabel "C"
	set yrange [-0.5:1.5]
	do for [j = 0:30:2 ] {
	plot sprintf("./advection_result%03d.txt", j) with linespoints lw 3 title sprintf("time = %f s", j*dt)}