背景
流体に限りませんが，数値計算の分野では，作成したプログラムの検証に困ることがよくあります．なぜなら，実験では観察できないような流れの詳細を捉える目的で数値計算が行われることが多く，そのような場合には数値計算結果と比較する実験結果が存在しないからです．ドキュメンタリーやドラマ，マンガ，アニメなどでは，すごさを示すために「スーパーコンピュータによるシミュレーションを･･･」という文言が用いられることがありますが，あれはあくまで演出です．スーパーコンピュータを使ったからといって正しい結果が得られる保証はありません．
プログラムの正しさを評価するために，テスト問題を解いて正しく解けているかを確認することが行われます．厳密解が分かっている問題を解いて，得られた結果と厳密解を比較し，正しくプログラムを実装できているか，実装した計算方法の特性はどうなのかを確認します．つまり，本当に計算したい問題に近いテスト問題を見つけてそれを計算し，テスト問題が正しく解けているからプログラムは正しそうだと判断するのです．
本記事では，悪名高いテスト問題の一つであるGresho渦をとりあげ，Gresho渦を計算するための初期値を導出します．
Gresho渦
Gresho渦は，2次元非粘性流体中で回転する単一渦です．遠心力と圧力勾配が釣り合うことで時間的に変化しない定在渦となります．流れの非圧縮性を仮定すると，Gresho渦の周方向速度$u_\theta$は次式で表されます．
u_\theta = \left\{
\begin{array}{ll}
5r & 0 \leq r< 0.2\\
2-5r & 0.2 \leq r < 0.4\\
0 & r \geq 0.4
\end{array}
\right.
半径方向速度$u_r$は$0$です．つまり，Gresho渦は軸対称流れです．非圧縮非粘性流れの支配方程式を極座標系で表せば計算は簡単なのですが，この軸対称流れをいわゆる直交座標系でどれだけ正しく解けるか（軸対称をどれだけ保持できるかや周方向運動量保存性）が評価のポイントです．
速度以外の分布も求められており，圧力$p$は
p = \left\{
\begin{array}{ll}
5+\frac{25}{2}r^2 & 0 \leq r< 0.2\\
9 – 4 \ln 0.2 + \frac{25}{2} r^2 – 20 r + 4 \ln r & 0.2 \leq r < 0.4\\
3 + 4 \ln 2 & r \geq 0.4
\end{array}
\right.,
渦度$\omega$は
\omega = \left\{
\begin{array}{ll}
10 & 0 \leq r< 0.2\\
\frac{2}{r}-10 & 0.2 \leq r < 0.4\\
0 & r \geq 0.4
\end{array}
\right.
と与えられます．
グラフからもわかりますが，渦度場は$r=0.2, r=0.4$で大きな不連続が存在します．
流れ関数
流れ関数とは，2次元非圧縮性流れにおいてのみ存在する量で，とある2点の流れ関数値の差がその2点間を通過する流量に相当します．
2次元非圧縮性流れの計算には，速度と圧力を用いる方法だけでなく，渦度と流れ関数を用いる方法があります．速度と圧力を用いる方法では，Gresho渦を計算する上で速度分布も圧力分布も分かっているので初期値を正しく与えられますが，渦度と流れ関数を使う方法では，渦度分布はわかっているものの，流れ関数分布が不明です．
本記事では，Gresho渦を表す流れ関数$\psi$の分布を導出します．
仮定
流れ関数を導出するにあたり，二つの仮定を導入します．
- 周方向に変化しない（$\psi = \psi(r)$）と仮定します．流れ関数が周方向に変化すると，半径方向速度$u_r$が$0$でなくなってしまうためです．
- 流れ関数は区間の両端で不連続ではない．今回導出する$\psi$の初期値は，数値計算に使うことを目的にしているので，できるかぎり取り扱いを簡便にします．
流れ関数の導出
流れ関数$\psi$と周方向速度$u_\theta$および渦度$\omega$の関係はそれぞれ次式で表されます．
u_\theta = -\frac{1}{r}\frac{\partial}{\partial r}\left(r\psi\right)
\omega = -\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial \psi}{\partial r}\right)
各区間の渦度を積分し，$\psi$を導出します．
$0 \leq r < 0.2$の区間における渦度の式を積分すると
\psi = -\frac{5}{3}r^2+\frac{C_1}{r}+D_1
が得られます．$C_1, D_1$は積分定数です．$\psi$を$u_\theta$の式に代入すると
u_\theta = 5r-\frac{D_1}{r}
となり，$D_1=0$であることがわかります．しかし，$C_1$は$u_\theta$の式に代入しても$\omega$の式に代入しても計算過程で微分されて消えるため，値が特定できません．
次は$0.2 \leq r < 0.4$の区間です．
渦度の式を積分すると
\psi = \frac{5}{3}r^2-r+\frac{C_2}{r}+D_2
となって，同様に$\psi$を$u_\theta$の式に代入すると，
u_\theta = -5r+2-\frac{D_2}{r}
より$D_2=0$であることがわかりますが，もう一つの積分定数$C_2$は計算過程で微分されて消えるため，これも特定できません．
$r\geq 0.4$では
\psi = \frac{C_3}{r}+D_3
を得ます．もう予想できると思いますが，$u_\theta$の式に$\psi$を代入した結果は
u_\theta = -\frac{D_3}{r}
となり，$D_3=0$と決まります．$C_3$は相変わらず分かりません．
とりあえずここまでに得られた$\psi$をまとめます．
\psi = \left\{
\begin{array}{ll}
-\frac{5}{3}r^2+\frac{C_1}{r} & 0 \leq r< 0.2\\
\frac{5}{3}r^2-r+\frac{C_2}{r} & 0.2 \leq r < 0.4\\
\frac{C_3}{r} & r \geq 0.4
\end{array}
\right.
積分定数の決定
上で求めた流れ関数の式は，それが定められた区間でだけ成立すればよいため，積分定数は任意にとることができます．しかしながら，流れ関数を初期値として用いることを考えると，積分定数を任意にとることで生じる区間の端での不連続は，計算の破綻につながります．そのような不連続を避けるように格子を切ったり，不連続があるところでif文を用いて対処することは可能かも知れませんが，それはGresho渦に限られ，計算手法の特性を正しく評価したことにはなりません．
そこで，積分定数が任意であることを利用し，区間の端で分布が不連続にならないように$C_1,C_2,C_3$を定めます1．そのために仮定の2を設けていたのです．$r=0.2, 0.4$で$\psi$は不連続的に変化しないので，区分関数のどちらの式を使っても値が同じとします．
\begin{array}{ll}
-\frac{5}{3}r^2+\frac{C_1}{r} = \frac{5}{3}r^2-r+\frac{C_2}{r} & \mathrm{at}\quad r= 0.2\\
\frac{5}{3}r^2-r+\frac{C_2}{r}=\frac{C_3}{r} & \mathrm{at}\quad r = 0.4
\end{array}
を満たす$C_1, C_2, C_3$を求めます．
$C_1$については，$C_1\neq 0$の時に$\psi$が$\infty$あるいは$-\infty$に発散するので，$C_1=0$を採用します．そうすると，順次$C_2,C_3$が求まります．
\psi = \left\{
\begin{array}{ll}
-\frac{5}{3}r^2 & 0 \leq r< 0.2\\
\frac{5}{3}r^2-r+\frac{1}{75r} & 0.2 \leq r < 0.4\\
-\frac{1}{25r} & r \geq 0.4
\end{array}
\right.
これで流れ関数の分布が定まりました．求められた分布を図示すると，確かに値の不連続がなく，なめらかに分布しているように見えます．
速度および渦度の離散的な導出
流れ関数の分布が求められたので，$\psi$を離散化し，差分法によって速度$u_\theta$および渦度$\omega$の初期値を計算します．$u_\theta$と$\psi$の関係式を見ると，$r$について2次式の$\psi$に$r$がかけられており，被微分関数としては$r$の3次式であることがわかります．$\omega$と$\psi$の関係式では，$\psi$が$r$で微分された後に$r^2$がかけられており，これも$r$の3次式であることがわかります．つまり，3次精度以上の差分式を用いて離散化することで，正確な値が得られることが分かります2．
\begin{eqnarray}
u_{\theta}&=& -\frac{1}{r}\frac{\partial}{\partial r}\left(r\psi\right)\\
u_{\theta i}&\simeq& -\frac{1}{r_i}\frac{-r_{i+2}\psi_{i+2}+8r_{i+1}\psi_{i+1}-8r_{i-1}\psi_{i-1}+r_{i-2}\psi_{i-2}}{12\varDelta r}
\end{eqnarray}
$r_i$は座標値，下付き添字の$i$は格子点の番号です．
$\omega$は4次精度の中心差分を2回適用します．
\begin{eqnarray}
\omega &=& -\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial \psi}{\partial r}\right)
\\
\varDelta \psi _i &=& r^2_i\frac{-\psi_{i+2}+8\psi_{i+1}-8\psi_{i-1}+\psi_{i-2}}{12\varDelta r}\\
\omega_i&\simeq& -\frac{1}{r^2_i}\frac{-\varDelta\psi_{i+2}+8\varDelta\psi_{i+1}-8\varDelta\psi_{i-1}+\varDelta\psi_{i-2}}{12\varDelta r}
\end{eqnarray}
これはよくない方法ですが，手軽さを重視して使っています．計算結果の質を高めたい場合にはこの方法は使うべきではありません．
数値実験
実際に初期値を計算してみましょう．計算環境は次のとおりです．
- Anaconda3 4.2.0.0 64bit
- Python 3.5.2
- Jupyter 4.2.0
# coding:utf-8
import numpy as np
from matplotlib import pyplot as plt
R=0.5 #計算領域
N=51 #格子点数
dr = R/(N-1) #格子間隔
r = np.linspace(0,R,N) #格子点の座標値
# 流れ関数
def ψ_(r):
if 0<=r<0.2:
return -5./3.*r**2
if 0.2<=r<0.4:
return 5./3.*r**2-r+1./(75.*r)
if r>=0.4:
return -1./(25.*r)
return 0.
ψ = np.vectorize(ψ_)(r)
# 4次精度中心差分を用いた周方向速度の計算
# 格子点が使えない両端2点は計算しない
u_numerical = np.zeros_like(r)
u_numerical[2:-2] = -1/r[2:-2]*(-r[4:]*ψ[4:]+8*r[3:-1]*ψ[3:-1]-8*r[1:-3]*ψ[1:-3]+r[:-4]*ψ[:-4])/(12*dr)
# 4次精度中心差分を用いた渦度の計算
# 中心差分を2回実行するので，計算しない点が両端4点に広がる
dψ = np.zeros_like(r)
dψ[2:-2] = r[2:-2]*r[2:-2] * (-ψ[4:]+8*ψ[3:-1]-8*ψ[1:-3]+ψ[:-4])/(12*dr)
ω_numerical = np.zeros_like(r)
ω_numerical[4:-4] = -1/r[4:-4]**2 *(-dψ[6:-2]+8*dψ[5:-3]-8*dψ[3:-5]+dψ[2:-6])/(12*dr)
計算された値をmatplotlibでプロットすると，それぞれ次のようになります．
図には正しい値が赤の実線で併記されています．区間の端で勾配が不連続になるため，厳密な値からは若干離れていますが，それ以外の場所では分布をよく再現しています．
試みに，$C_1,C_2,C_3$全てを$0$にした$\psi$を表示してみると，確かに区間の端で値が不連続に変化しています．
この$\psi$を使って$u_\theta$と$\omega$を計算してみると，区間の端近く以外では厳密な値と一致しているようですが，区間をまたいで格子点を参照する箇所ではあまりにもひどい結果しか得られません．
このような初期値を使って計算をしても，正しい結果が得られるとは到底考えられないので，本記事で示した流れ関数の分布が初期値としてよりよいことが確認できます．
まとめ
本記事では，Gresho渦を表す流れ関数の分布として，以下の式を提示しました．
\psi = \left\{
\begin{array}{ll}
-\frac{5}{3}r^2 & 0 \leq r< 0.2\\
\frac{5}{3}r^2-r+\frac{1}{75r} & 0.2 \leq r < 0.4\\
-\frac{1}{25r} & r \geq 0.4
\end{array}
\right.
この分布を用いることで分布の不連続性を考慮する必要がなくなり，周方向速度および渦度の分布が離散的に求められることが確認できました．