概要
文献12を参考に,2次元非圧縮性流体の数値解析プログラムをFortranで作成した.2次元周期流を仮定し,渦度輸送方程式をスペクトル法で解いた.離散フーリエ変換にはIntel MKL3を用いた.
$1024^2$自由度の解析.Paraview4で可視化した.
記事内へのリンク
支配方程式の導出
初期条件の与え方
2次元DFT
非線形項の評価(変換法)
時間積分法
計算例
ソースコード
更新履歴
2019.12.01 : 公開
支配方程式の導出
2次元の非圧縮性流体を仮定し,スペクトル法で解く支配方程式を導出する.156
渦度輸送方程式の導出
基礎方程式は次に示す連続の式(質量保存則)とNavier-Stokes方程式(運動量保存則)である.56
連続の式
\begin{align}
\nabla\cdot\mathbf{u}
=
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}
=
0
\tag{1}
\end{align}
Navier-Stokes方程式
$$
\begin{align}
\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}
=-\nabla p +\nu\nabla^2\mathbf{u}
\tag{2}\
\end{align}
$$
諸量は代表長さ・速さで無次元化している.ここで,$p=\rm{静圧/密度}$,$\nu=1/\rm{Re}=const.$である.$x$方向速度$u$, $y$方向速度$v$についてそれぞれ書き下せば,次のようになる.
$$
\begin{align}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}
=-\frac{\partial p}{\partial x}
+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)
\tag{3}
\end{align}
$$
$$
\begin{align}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}
=-\frac{\partial p}{\partial y}
+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)
\tag{4}
\end{align}
$$
渦度輸送方程式
渦度$\zeta$を導入し,式(2)を扱いやすい形に変形する.ここで$\zeta$の定義は,
$$
\begin{align}
\zeta=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}
\end{align}
$$
であり,3次元で定義される渦度$\omega=\nabla\times\mathbf{u}$の第3成分と一致する.
式(4)を$x$で偏微分した式から,式(3)を$y$で偏微分した式を引くと,
\begin{align}
L.H.S.
&=
\frac{\partial}{\partial x}\left(
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}
\right)
-
\frac{\partial}{\partial y}\left(
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}
\right)\\
&=
\frac{\partial \zeta}{\partial t}
+u\frac{\partial \zeta}{\partial x}+v\frac{\partial \zeta}{\partial y}\\
R.H.S.
&=
\frac{\partial}{\partial x}\left(-\frac{\partial p}{\partial y}
+\nu\nabla^2v
\right)
-
\frac{\partial}{\partial y}\left(-\frac{\partial p}{\partial x}
+\nu\nabla^2u
\right)\\
&=\nu\nabla^2\zeta
\end{align}
よって,渦度$\zeta$に対する方程式(渦度輸送方程式)が得られる.
$$
\begin{align}
\frac{D \zeta}{D t}=\frac{\partial \zeta}{\partial t}
+u\frac{\partial \zeta}{\partial x}+v\frac{\partial \zeta}{\partial y}
=\nu\nabla^2\zeta
\tag{5}
\end{align}
$$
ここで2次元のラプラシアンと実質微分は次のように定義される.
$$
\begin{align}
\nabla^2=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial x^2}
,;;
\frac{D}{Dt}=\frac{\partial }{\partial t}
+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}.
\end{align}
$$
流れ関数の導入
式(5)に流れ関数$\psi$を導入する.$\psi$を文献1に従い,次のように定義する.
$$
\begin{align}
u=-\frac{\partial \psi}{\partial y},;;v=\frac{\partial \psi}{\partial x}.
\end{align}
$$
一般的な定義5と符号が逆であるが,このように定義することで渦度$\zeta$との関係式が次のようになる.
$$
\begin{align}
\zeta=\nabla^2\psi
\tag{6}
\end{align}
$$
よって,渦度輸送方程式(式(5))は$\psi$のみに関する偏微分方程式となる.
$$
\begin{align}
\frac{D (\nabla^2\psi)}{D t}
=\nu\nabla^2(\nabla^2\psi)
\tag{7}
\end{align}
$$
フーリエ空間での支配方程式
境界条件
計算領域を$0\leq x,y \leq 2\pi$として,周期境界条件をもちいる.すなわち任意の時間$t$における境界の$\psi$の値は,
$$
\begin{align}
\psi(0,y,t)=\psi(2\pi,y,t),;;\psi(x,0,t)=\psi(x,2\pi,t)\
\end{align}
$$
を満たす.
関数展開
$x,y$の両方向に周期的なので,展開関数としてフーリエ級数を選択する.$\psi(x,y,t)$は次のように定義される.
$$
\begin{align}
\psi(x,y,t)=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
a_{kl}(t)e^{ikx}e^{ily}
\tag{8}
\end{align}
$$
展開係数$a_{kl}(t)$の自由度は$\psi$が実数であることから,
$$
\begin{align}
a_{(-k)(-l)}=a^*_{kl}
\end{align}
$$
の制約がかかり,$N^2$となる.ここで,$^*$は複素共役を意味する.
\begin{align}
\because \;\;{\rm Im}\left[a_{kl}e^{i(kx+ly)}-a_{(-k)(-l)}e^{-i(kx+ly)}\right]=0
\end{align}
常微分方程式の導出
展開係数$a_{kl}$に対する常微分方程式を導く.式(7)の残差の式
$$
R(x,y,t)=\frac{D (\nabla^2\psi)}{D t}
-\nu\nabla^2(\nabla^2\psi)
$$
に$\psi$の展開式(式(8))を代入し,$e^{-ikx}e^{-ily}$をかけて,計算領域全域で積分したものが$0$になるようにすれば良い.すなわち,次式を要請する.
$$
\int_{0}^{2\pi}\int_{0}^{2\pi}R(x,y,t)e^{-ikx}e^{-ily}dxdy=0
$$
時間微分項
\begin{align}
&\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{\partial(\nabla^2\psi)}{\partial t}
e^{-ikx}e^{-ily}dxdy\\
&=
\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{\partial}{\partial t}\left[\nabla^2\left(\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
a_{kl}(t)e^{ikx}e^{ily}\right)\right]e^{-ikx}e^{-ily}dxdy\\
&=
\sum_k\sum_l\left\{-(k^2+l^2)\frac{d a_{kl}}{d t}\right\}
\int_{0}^{2\pi}\int_{0}^{2\pi}dxdy\\
&=
-(2\pi)^2\sum_k\sum_l(k^2+l^2)\frac{d a_{kl}}{d t}
\end{align}
線形項
\begin{align}
&\int_{0}^{2\pi}\int_{0}^{2\pi}\nu\nabla^2(\nabla^2\psi)
e^{-ikx}e^{-ily}dxdy\\
&=-(2\pi)^2\sum_k\sum_l(k^2+l^2)^2a_{kl}
\end{align}
非線形項
変換法で評価するので,$\psi$に展開式を代入して変形したりはしない.また,上述の2項で付随した$(2\pi)^2$で割ると,
$$
\begin{align}
F_{kl}=\frac{1}{(2\pi)^2}\int_0^{2\pi}\int_0^{2\pi}
\left(
u\frac{\partial \zeta}{\partial x}+v\frac{\partial \zeta}{\partial y}
\right)e^{-ikx}e^{-ily}dxdy
\tag{9}
\end{align}
$$
となり,$F_{kl}$は非線形項の展開係数とみなせる.(かつ2次元のFFTが適用できそうな式形になる)
常微分方程式
以上を踏まえて,流れ関数$\psi$の展開係数$a_{kl}$に対する常微分方程式は,次のようになる.
$$
\begin{align}
\frac{d a_{kl}}{d t}=\frac{F_{kl}}{(k^2+l^2)}-\nu(k^2+l^2)a_{kl}
\tag{10}
\end{align}
$$
物理空間における格子点
物理空間の自由度も同様に$N^2$として,$0\leq x,y \leq 2\pi$の計算領域に格子点$(x_m,y_n)$を置く.
$$
\begin{align}
x_m=\frac{2\pi}{N}m;,;;;y_n=\frac{2\pi}{N}n;;(0\leq m,n\leq N-1)
\tag{11}
\end{align}
$$
初期条件
式(10)は流れ関数$\psi$の展開係数$a_{kl}$に関する常微分方程式である.しかし,流れ場は$\psi$よりも渦度$\zeta$で定義する方がわかりやすい(と思う).よって,$\zeta(x,y,0)$で与えて,$\psi(x,y,0)$に変換する.$\zeta$の展開係数を$b_{kl}$と置くと,式(8)と同様に
$$
\begin{align}
\zeta(x,y,0)=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
b_{kl}(0)e^{ikx}e^{ily}
\tag{12}
\end{align}
$$
と展開できる.後に述べる2次元FFTで$\zeta\rightarrow b_{kl}$の変換ができる.
$$
\begin{align}
b_{kl}=\frac{1}{2\pi}\int_0^{2\pi}\int_0^{2\pi}\zeta e^{-ikx}e^{-ily}dxdy
\tag{13}
\end{align}
$$
流れ関数への変換
式(6)で示された$\zeta$と$\psi$の関係は,2次元のPoisson方程式である.
$$
\begin{align}
\zeta(x,y,0)=\nabla^2\psi(x,y,0)=\frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2}
\tag{14}
\end{align}
$$
スペクトル法では,Poisson方程式は反復法を用いずに解くことができる.上式から残差を定義し,重み$e^{-ikx}e^{-ily}$をかけて計算領域で積分する.式(8),(11),(12)を用いれば,
$$
\begin{align}
a_{kl}=-\frac{b_{kl}}{k^2+l^2}
\tag{15}
\end{align}
$$
が得られる.よって$b_{kl}$から$a_{kl}$を求められる.
2次元DFT
定義
まず,順方向・逆方向のDFTを定義する.格子点$(x_m,y_n)$を定義して,
$$
\begin{align}
x_m=\frac{2\pi}{M}m;;(0\leq m\leq M-1),;;;
y_n=\frac{2\pi}{N}n;;(0\leq n\leq N-1)
\end{align}
$$
関数$g(x,y)$をその格子点上で定義する.
$$
g_{mn}=g(x_m,y_n)
$$
正変換
\begin{align}
G_{kl}
=\frac{1}{MN}
\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}
g_{mn}
\exp\left(-2\pi i \frac{m\cdot k}{M}\right)
\exp\left(-2\pi i \frac{n\cdot l}{N}\right)\tag{16}\\
(0\leq k\leq M-1,\;\; 0\leq l\leq N-1)
\end{align}
逆変換
\begin{align}
g_{mn}=\sum_{k=0}^{M-1}\sum_{l=0}^{N-1}G_{kl}
\exp\left(2\pi i \frac{m\cdot k}{M}\right)
\exp\left(2\pi i \frac{n\cdot l}{N}\right)\tag{17}\\
(0\leq m\leq M-1,\;\; 0\leq n\leq N-1)
\end{align}
実DFT
$g$を実数とする,2次元実DFTは,1次元複素DFTと1次元実DFTを組み合わせとみなせる(多次元FFTのアルゴリズムは必ずそうとは限らないが).例えば逆変換(式(17))は,
- 複素IDFT.中間変数$g'_{kn}$は複素数.$(0\leq k\leq M-1,;; 0\leq n\leq N-1)$
\begin{align}
g'_{kn}&=\sum_{l=0}^{N-1}G_{kl}
\exp\left(2\pi i \frac{n\cdot l}{N}\right)\\
\end{align}
- 実IDFT.$g_{mn}$は実数.$(0\leq m\leq M-1)$
$$
\begin{align}
g_{mn}&=\sum_{k=0}^{N-1}g'_{kn}
\exp\left(2\pi i \frac{m\cdot k}{M}\right)\
\end{align}
$$
級数展開式の変形
例えば,式(8)で示された$a_{kl}\rightarrow \psi$の変換に式(17)の逆変換を適用することを考える.
(再掲)
$$
\begin{align}
\psi(x,y,t)=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
a_{kl}(t)e^{ikx}e^{ily}
\tag{8}
\end{align}
$$
$x,y$も式(11)のように離散化する.
$$
\begin{align}
\psi(x_m,y_n)=\psi_{mn}=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
a_{kl}(t)
\exp\left(2\pi i \frac{m\cdot k}{N}\right)
\exp\left(2\pi i \frac{n\cdot l}{N}\right)\tag{18}
\
(m,n=0,..,N-1)
\end{align}
$$
問題になるのは,$k,l$のとる範囲である.
まず,内側の$\sum_l$ : $y$方向の級数展開について,
\begin{align}
\sum_{l=-N/2}^{-1}a_{kl}\exp\left(2\pi i \frac{l\cdot n}{N}\right)
&=\sum_{l=N/2}^{N-1}a_{k(l-N)}\exp\left(2\pi i \frac{(l-N)\cdot n}{N}\right)\\
&=\sum_{l=N/2}^{N-1}a_{k(l-N)}\exp\left(2\pi i \frac{l\cdot n}{N}\right)\\
\end{align}
より,中間変数$\psi'_{kn}$(複素数)は
\begin{align}
\psi'_{kn}
&=\sum_{l=-N/2}^{N/2-1}a_{kl}\exp\left(2\pi i \frac{l\cdot n}{N}\right)\\
&=\sum_{l=0}^{N/2-1}a_{kl}\exp\left(2\pi i \frac{l\cdot n}{N}\right)
+\sum_{l=N/2}^{N-1}a_{k(l-N)}\exp\left(2\pi i \frac{l\cdot n}{N}\right)\\
&=\sum_{l=0}^{N-1}w_{kl}\exp\left(2\pi i \frac{l\cdot n}{N}\right)
\end{align}
ここで,$w_{kl}$の定義は,
\begin{equation}
w_{kl}=
\begin{cases}
a_{kl}& 0\leq l \leq N/2-1\\
a_{k(l-N)}& N/2 \leq l \leq N-1
\end{cases}
\end{equation}
となる.よって$\sum_l$は長さ$N$の複素IDFTに変形できた.
次に,$\sum_k$ : $x$方向の級数展開について
$$
\begin{align}
\psi_{mn}
=\sum_{k=-N/2}^{N/2-1}
\psi'_{kn}
\exp\left(2\pi i \frac{k\cdot m}{N}\right)\
\end{align}
$$
である.$\psi_{mn}$は実数なので,$\psi_{(-k)n}'=\psi^{\prime *}_{kn}$の制約が生じる.
\begin{align}
\psi_{mn}
&=\sum_{k=-N/2}^{N/2-1}
\psi'_{kn}
\exp\left(2\pi i \frac{k\cdot m}{N}\right)\\
&=\sum_{k=0}^{N/2-1}
\psi'_{kn}
\exp\left(2\pi i \frac{k\cdot m}{N}\right)
+\sum_{k=1}^{N/2}
\psi'_{(-k)n}
\exp\left(2\pi i \frac{(-k)\cdot m}{N}\right)\\
&=\sum_{k=0}^{N/2-1}
\psi'_{kn}
\exp\left(2\pi i \frac{k\cdot m}{N}\right)
+\sum_{k=1}^{N/2}
\psi'^*_{kn}
\exp\left(-2\pi i \frac{k\cdot m}{N}\right)\\
&=\psi_{0n}
+2\sum_{k=1}^{N/2-1}\left[{\rm Re}(\psi'_{kn})\cos\left(2\pi\frac{k\cdot m}{N}\right)-{\rm Im}(\psi'_{kn})\sin\left(2\pi\frac{k\cdot m}{N}\right) \right]
\end{align}
自由度$N$の(1次元)実IDFTは通常,上式に$k=N/2$の項を加えた次の式を計算するように(FFTで)実装されている1.
$$
\begin{align}
f(m)
&=F_{0}
+\left(-1\right)^mF_{N/2}
+2\sum_{k=1}^{N/2-1}\left[{\rm Re}(F_{m})\cos\left(2\pi\frac{k\cdot m}{N}\right)
-{\rm Im}(F_{m})\sin\left(2\pi\frac{k\cdot m}{N}\right) \right]
\end{align}
$$
よって上式で$F=\psi',;;F_{N/2}=0$とすれば,$\sum_k$の項は長さ$N$の実IDFTを用いて計算できる.この逆は,$x$方向の実DFT$\rightarrow$$y$方向の複素DFTを行えば良い.
結局,式(8),(13)のような"フーリエ空間で定義される展開係数"$\leftrightarrow$"物理空間で定義される値"の変換は,2次元DFTを利用して計算できる.
非線形項の評価
次式で定義される展開係数$F_{kl}$の効率的な評価方法を述べる.
(再掲)
$$
\begin{align}
F_{kl}
&=\frac{1}{(2\pi)^2}\int_0^{2\pi}\int_0^{2\pi}
\left(
u\frac{\partial \zeta}{\partial x}+v\frac{\partial \zeta}{\partial y}
\right)e^{-ikx}e^{-ily}dxdy\
\tag{9}
\end{align}
$$
変換法
フーリエ空間で計算するよりも,IFFTを使って物理空間上の格子点で値を評価しFFTでフーリエ空間に戻す方が計算量が削減できる.ただし,この時aliasing誤差が生じるため,これを取り除く工夫をしなくてはならない.
計算量の削減
非線形項を次のように変形することで,FFTの回数を減らすことができる1.
\begin{align}
u\frac{\partial \zeta}{\partial x}+v\frac{\partial \zeta}{\partial y}
&=\frac{\partial u\zeta}{\partial x}+\frac{\partial v\zeta}{\partial y}\\
&=\frac{\partial }{\partial x}
\left\{
u\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)
\right\}
+\frac{\partial }{\partial y}
\left\{
v\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)
\right\}\\
&=\frac{\partial }{\partial x}
\left(
\frac{\partial uv}{\partial x}
+\frac{1}{2}\frac{\partial v^2}{\partial y}
-\frac{1}{2}\frac{\partial u^2}{\partial y}\right)
+\frac{\partial }{\partial y}
\left(
\frac{1}{2}\frac{\partial v^2}{\partial x}
-\frac{1}{2}\frac{\partial u^2}{\partial x}
-\frac{\partial uv}{\partial y}\right)\\
&=\left(\frac{\partial^2}{\partial x^2}-\frac{\partial^2}{\partial y^2}\right)(uv)
+\frac{\partial^2}{\partial x \partial y}(v^2-u^2)
\end{align}
$uv,v^2-u^2$それぞれに対する展開係数を$A_{kl},B_{kl}$と置けば,
\begin{align}
F_{kl}=-(k^2-l^2)A_{kl}-klB_{kl}\tag{19}
\end{align}
と計算できる.
変形前
$a_{kl}\rightarrow u,v,\zeta_x,\zeta_y$でIFFT4回と$(u\zeta_x+v\zeta_y)\rightarrow F_{kl}$のFFT1回:合計5回
変形後
$a_{kl}\rightarrow u,v$でIFFT2回と$uv\rightarrow A_{kl},v^2-u^2\rightarrow B_{kl}$のFFT2回:合計4回
aliasing誤差の除去
2次元における2次の非線形項の変換法による評価でどのようなaliasing項が表れるかを確認し,そのaliasing項の除去方法について述べる7689.
aliasing項
例として$w(x,y)=u(x,y)v(x,y)$の評価を考える.$u,v$を有限なフーリエ級数で展開し,展開係数$\hat{u}$, $\hat{v}$を定義する.
$\hat{u},\hat{v}\rightarrow \boxed{{\rm IFFT}} \rightarrow u,v\rightarrow w\rightarrow \boxed{{\rm FFT}} \rightarrow \hat{w}$を考える.
\begin{align}
&u(x_m,y_n)=u_{mn}=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
\hat{u}_{kl}e^{ikx_m}e^{ily_n},\\
&v(x_m,y_n)=v_{mn}=\sum_{k=-N/2}^{N/2-1}\sum_{l=-N/2}^{N/2-1}
\hat{v}_{kl}e^{ikx_m}e^{ily_n}.\tag{20}
\end{align}
$w_{mn}=u_{mn}v_{mn}$より,物理空間の格子点上で得られた$w_{mn}$を2次元DFTすれば,$\hat{w}_{kl}$が得られるように思える.しかし実際にはaliasing項による誤差が生じる.
\begin{align}
\widetilde{w}_{kl}
=\frac{1}{N^2}\sum_{m=0}^{N-1}\sum_{n=0}^{N-1}
w_{mn}e^{-ikx_m}e^{-ily_n}\\
\end{align}
上式に式(20)を代入する.
\begin{align}
\widetilde{w}_{kl}
=\frac{1}{N^2}\sum_{m}\sum_{n}\sum_{k_1}\sum_{k_2}\sum_{l_1}\sum_{l_2}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}e^{i(k_1+k_2-k)x_m}e^{i(l_1+l_2-l)y_n}\tag{21}
\end{align}
各$\sum$の範囲は省略した.ここで,$e^{ikx_m}$について考える.式(11)の格子点座標の定義を代入すると,
\begin{align}
\sum_{m=0}^{N-1}e^{ikx_m}
=
\sum_{m=0}^{N-1}\exp\left(ik\frac{2\pi }{N}m\right)
=
\sum_{m=0}^{N-1}\exp\left(\frac{2\pi ik}{N}m\right)\tag{22}
\end{align}
初項$1$,公比$r=\exp(2\pi ik/N)$,項数$N$の等比級数とみなせる.よって,等比級数の和の公式より
\begin{align}
\sum_{m=0}^{N-1}e^{ikx_m}
=
\begin{cases}
0& r\ne1\\
N& r=1
\end{cases}
\end{align}
当然$e^{ily_n}$も同様に考えられる.$C_1,C_2$を整数として
\begin{align}
k_1+k_2-k=C_1N,\;\;l_1+l_2-l=C_2N
\end{align}
の場合のみ値を持つ.$k,l$の取りうる範囲から,$C$の取りうる範囲は,$C_1,C_2=0,\pm1$.結局,次のように書き下すことができる.
\begin{align}
\widetilde{w}_{kl}
&=\sum_{C_1,C_2=0,\pm1}\sum_{k_1,k_2,l_1,l_2}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}
\delta_{k_1+k_2-k,C_1N}\delta_{l_1+l_2-l,C_2N}\\
&=\sum_{C_1,C_2=0,\pm1}\sum_{k_1+k_2-k=C_1N,\\l_1+l_2-l=C_2N}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}\\
&=\sum_{k_1+k_2-k=0,\\l_1+l_2-l=0}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}
+\sum_{k_1+k_2-k=\pm N,\\l_1+l_2-l=0}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}
+\sum_{k_1+k_2-k=0,\\l_1+l_2-l=\pm N}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}
+\sum_{k_1+k_2-k=\pm N,\\l_1+l_2-l=\pm N}
\hat{u}_{k_1l_1}\hat{v}_{k_2l_2}\tag{23}
\end{align}
第一項が求めたい$\hat{w}$で,第二,三項がsingle-aliased contributions,第四項がdouble-aliased contributionと呼ばれる7.
de-aliasing method
式(23)で表れた3つのaliasing項を取り除く.
$$
-N<k_1+k_2-k<N,;;;-N<l_1+l_2-l<N\tag{24}
$$
この二つの条件を満たせば良い.つまり,$k_1,k_2,k,l_1,l_2,l$の絶対値が大きくなると条件式(24)が満たされない可能性がある.よって,絶対値の大きな波数に対する成分を強制的に$0$とすれば良い.ここで,非線形項を評価する物理空間の自由度を変えずフーリエ空間の有効自由度を減らす方法を2/3-rule.物理空間の自由度を増やしてフーリエ空間の自由度を保つ方法を3/2-ruleと呼ぶ.
2/3-rule
$-N'\leq k_1,k_2,k,l_1,l_2,l\leq N'$ならば式(24)を満たすと考えると,
$${\rm Max}(k_1+k_2-k,l_1+l_2-l)=3N'$$
$${\rm Min}(k_1+k_2-k,l_1+l_2-l)=-3N'$$
より,$3N'<N$ならば良い.最大有効波数は$N'=N/3=(2/3)\cdot(N/2)$となるので,2/3-ruleと呼ばれる.展開係数$\hat{u}$,$\hat{v}$を次の条件で打ち切ってからIFFTで$u,v$を求め,その積$w$をFFTにかければ,$\hat{w}=\widetilde{w}$となる.
\begin{align}
\hat{u}_{kl},\hat{v}_{kl}=0\;\;\; {\rm if}\;\;|k|>N/3 \cup|l|>N/3
\end{align}
この手法では,計算量はほぼ増加しないことがメリットであり,フーリエ空間の自由度が減少することがデメリットである.
3/2-rule
フーリエ空間での自由度が$N^2$であることを担保したい.2/3-ruleとは逆に,$-N/2\leq k,l\leq N/2-1$の区間の外側に$0$成分を付加して,aliasing項が表れないようにする.必要な$0$成分の範囲は,
$$
\begin{align}
-N'<\frac{3}{2}N-1\leq k_1+k_2-k,l_1+l_2-l\leq\frac{3}{2}N-2<N'
\end{align}
$$
よって$N'=3N/2$として,次図で示す中心部の緑色の領域を除く赤線で囲まれた領域を0詰めした展開係数$\hat{u}_{k'l'}$
\begin{align}
\hat{u}_{k'l'}=
\begin{cases}
\hat{u}_{kl}\:\:&-N/2\leq k,l\leq N/2-1\\
0&{\rm otherwise}
\end{cases}\\
(-3N/2\leq k',l'\leq 3N/2-1)
\end{align}
を用いて,自由度$(3N/2)^2$の2次元IFFTで$u_{m'n'},v_{m'n'}$を求める.ここで$m',n'$の範囲は$0\leq m',n'\leq 3N/2-1$である.$w_{m'n'}=u_{m'n'}v_{m'n'}$を求めて,自由度$(3N/2)^2$の2次元FFTで$\widetilde{w}$を求める.求めたい$\hat{w}$はその一部で,
\begin{align}
\hat{w}_{kl}=\widetilde{w}_{kl}\;\;(-N/2\leq k,l\leq N/2-1)
\end{align}
となり,$\widetilde{w}_{kl}$の残りの成分の情報は捨てることになる.この手法のメリットは,フーリエ空間の自由度が減少しないことであり,デメリットは計算コストの増加である.今回はこちらの手法を用いる.
時間積分法
古典的な4段4次精度Runge-Kutta法と積分因子法(integrating factor technique)7を組み合わせて時間積分を行う110.
4段4次Runge-Kutta法
次の微分方程式を4段4次のRunge-Kutta法を用いて積分することを考える.
$$
\begin{align}
\frac{df}{dt}=\mathcal{F}(t,f)
\end{align}
$$
時刻$t^n$で$f=f^n$とすると,$t^{n+1}=t^n+\Delta t$における$f^{n+1}$は次の様に評価できる.
\begin{align}f^{n+1}&=f^n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)+\mathcal{O}(\Delta t^5)\tag{25}\\k_1&=\Delta t\mathcal{F}(t^n,f^n)\\k_2&=\Delta t\mathcal{F}(t^n+\frac{1}{2}\Delta t,f^n+\frac{1}{2}k_1)\\k_3&=\Delta t\mathcal{F}(t^n+\frac{1}{2}\Delta t,f^n+\frac{1}{2}k_2)\\k_4&=\Delta t\mathcal{F}(t^n+\Delta t,f^n+k_3)\\\end{align}
積分因子法
次の微分方程式を考える.
$$
\begin{align}\frac{df}{dt}=\mathcal{F}(f)-\alpha f\end{align}
$$
ここで$\mathcal{F}(f)$は$f$に関する非線形関数,$\alpha$は定数係数である.これは以下の様に式変形できる.
\begin{align}\frac{df}{dt}+\alpha f&=\mathcal{F}(f)\\\left(\frac{df}{dt}+\alpha f\right)e^{\alpha t}&=\mathcal{F}(f)e^{\alpha t}\tag{26}\\\end{align}
この形のもとで,時間微分を差分化するのが積分因子法 (Integrating factor technique) である.例えば,1次精度のEuler陽解法を用いると次式のようになる.
$$
f^{n+1}=e^{-\alpha \Delta t}\left(f^n+\Delta t\mathcal{F}(f^n)\right)
$$
RK法と積分因子法の組み合わせ
式(26)に式(25)を用いると,次のように時間微分を差分化できる110.
\begin{align}f^{n+1}&=e^{-\frac{\alpha}{2}\Delta t}\left(e^{-\frac{\alpha}{2}\Delta t}(f^n+\frac{1}{6}k_1)+\frac{1}{3}(k_2+k_3)\right)+\frac{1}{6}k_4+\mathcal{O}(\Delta t^5)\\k_1&=\Delta t\mathcal{F}\left(t^n,f^n\right)\\k_2&=\Delta t\mathcal{F}\left(t^n+\frac{1}{2}\Delta t,e^{-\frac{\alpha}{2}\Delta t}\left(f^n+\frac{1}{2}k_1\right)\right)\\k_3&=\Delta t\mathcal{F}\left(t^n+\frac{1}{2}\Delta t,e^{-\frac{\alpha}{2}\Delta t}f^n+\frac{1}{2}k_2\right)\\k_4&=\Delta t\mathcal{F}\left(t^n+\Delta t,e^{-\frac{\alpha}{2}\Delta t}\left(e^{-\frac{\alpha}{2}\Delta t}f^n+k_3\right)\right)\\\end{align}
なお,線形項 ($\alpha f$) の評価を常に$\frac{\Delta t}{2}$の時間刻み幅で行うために,文献10の式を若干変形している.
計算例
3つの計算例を示す.共通して自由度$N=2^{10}$,4並列で計算した.Paraview4を用いて渦度を可視化した.
###Kelvin-Helmholtz不安定性
demo01 | Kelvin-Helmholtz不安定性 |
---|---|
$\nu$ | $0$(=非粘性) |
$\Delta t$ | $5.0e-3$ |
$t_{end}$ | 20 |
total step | 4000 |
計算時間 | 約12分 |
初期条件
\begin{align}
\zeta(x,y)=&4.0\left\{
\exp\left[-(20x-10\pi)^2\right]\cdot(1+0.01\cos(2.0(y+0.5\pi))\\
-\exp\left[-(20x-30\pi)^2\right]\cdot(1+0.01\cos(2.0(y+0.5\pi))
\right\}
\end{align}
可視化
###同符号の渦の合体
文献1にある計算例を用いた.
demo02 | 同符号の渦の合体 |
---|---|
$\nu$ | $2.0e-6$ |
$\Delta t$ | $2.0e-2$ |
$t_{end}$ | 100 |
total step | 5000 |
計算時間 | 約15分 |
初期条件
\begin{align}
\zeta(x,y)=\;\;&
\exp\left[\left(\cos(x-x_1)+\cos(y-y_1)-2\right)/\sigma^2\right]\\
-&\exp\left[\left(\cos(x-x_2)+\cos(y-y_2)-2\right)/\sigma^2\right]-C\\
\end{align}
なお定数は以下の通り与える
$$
\sigma=\pi/10,;;x_1=y_1=0.8\pi,;;x_2=y_2=1.2\pi,;;C=0.0322447
$$
可視化
###2次元等方性乱流の減衰
San,O.,Staples,A.E.11の数値実験(6. Two-dimensional decaying turbulence)を行う.
demo02 | 2次元乱流の減衰 |
---|---|
$\nu$ | $1.0e-3$ |
$\Delta t$ | $2.0e-4$ |
$t_{end}$ | 10 |
total step | 50000 |
計算時間 | 約2.5時間 |
初期条件
San,O.,Staples,A.E.11と同様
可視化
エネルギースペクトル
ソースコード
上で示した3つの計算例を含むソースコードはGithubで公開している.
https://github.com/toya42/psm4cfd
本プログラムの特徴は,次のようにまとめられる.
- 安全で移植が容易な記法で実装 (Fortran2008.macOS&Ubuntuで動作確認)
- 高速なFFT (Intel MKLの利用)
- 並列計算が可能(スレッド並列:対応,プロセス並列:未対応)
- 大きな自由度の計算も可能(コンパイルオプションの変更だけで$N\geq 2^{14}$に対応)
- MITライセンスで公開(自由に使ってください)
Intel MKL以外のFFTライブラリを使用したい場合
fft2d_mkl.F90
を置き換えれば良い.入出力は一般的なCCE形式12なので,FFTW13などへの変更は容易だと思われる.
終わりに
スペクトル法で非圧縮性流体の2次元周期流を計算するプログラムを作成した.V&Vが不十分,プロセス並列に未対応等と課題は残っているが,せっかくなのでGithubで公開した.
今後の予定としては
- 3次元周期流
- 2次元ポアズイユ流
- 3次元チャネル乱流
に挑戦したいと考えている.
-
C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Methods Fundamentals in Single Domain, Springer-Verlag (2006). ↩
-
Intel Math Kernel Library (https://software.intel.com/en-us/mkl) ↩
-
Paraview (https://www.paraview.org/) ↩ ↩2
-
大宮司久明 (2009) 『数値流体力学大全』(http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/) 『第12章 非圧縮性流れの解法―渦度輸送方程式と流れ関数方程式を解く方法』(http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/Chapter12.pdf) ↩ ↩2 ↩3
-
C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Methods Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer-Verlag (2007). ↩ ↩2 ↩3
-
中央大学理工学研究科物理学専攻中野研究室 講義資料 (https://www.phys.chuo-u.ac.jp/labs/nakano/tokuron2/sec61(08).pdf) 2019.09.04閲覧 ↩
-
岐阜大学准教授 田中雅宏(2006) 講義配布資料のまとめ (https://www1.gifu-u.ac.jp/~tanaka/numerical_analysis.pdf) 2019.09.04閲覧 ↩
-
北川真帆,村上洋一.数理解析研究所講究録(2012), 1800:226-235,京都大学数理解析研究所,https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/173013 ↩ ↩2 ↩3
-
San,O.,Staples,A.E.:High-order methods for decaying two-dimensional homogeneous isotropic turbulence. Comput. Fluids 63, 105–127 (2012) ↩ ↩2
-
Developer Reference for Intel® Math Kernel Library - Fortran DFTI_PACKED_FORMAT (https://software.intel.com/en-us/onemkl-developer-reference-fortran-dfti-packed-format), fortranに関するメモ FFT(http://www.ton.scphys.kyoto-u.ac.jp/~michikaz/fortran.html#level5) 等 ↩
-
FFTW (http://www.fftw.org/) ↩