LoginSignup
31

More than 3 years have passed since last update.

趣味でスペクトル法_#003_2次元周期流

Last updated at Posted at 2019-11-30

概要

 文献12を参考に,2次元非圧縮性流体の数値解析プログラムをFortranで作成した.2次元周期流を仮定し,渦度輸送方程式をスペクトル法で解いた.離散フーリエ変換にはIntel MKL3を用いた.

image-0100.jpeg

$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}$の残りの成分の情報は捨てることになる.この手法のメリットは,フーリエ空間の自由度が減少しないことであり,デメリットは計算コストの増加である.今回はこちらの手法を用いる.

padding.jpg

時間積分法

 古典的な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}

可視化

demo01.gif

同符号の渦の合体

文献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
$$

可視化

demo02.gif

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と同様

可視化

demo03.gif

エネルギースペクトル

es.jpg

ソースコード

 上で示した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次元チャネル乱流

に挑戦したいと考えている.


  1. 石岡圭一(2004)『スペクトル法による数値計算入門』東京大学出版会  

  2. C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Methods Fundamentals in Single Domain, Springer-Verlag (2006). 

  3. Intel Math Kernel Library (https://software.intel.com/en-us/mkl

  4. Paraview (https://www.paraview.org/

  5. 大宮司久明 (2009) 『数値流体力学大全』(http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/) 『第12章 非圧縮性流れの解法―渦度輸送方程式と流れ関数方程式を解く方法』(http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/Chapter12.pdf

  6. 木田重雄,柳瀬眞一郎 (1999) 『乱流力学』朝倉出版 

  7. C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Methods Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer-Verlag (2007). 

  8. 中央大学理工学研究科物理学専攻中野研究室 講義資料 (https://www.phys.chuo-u.ac.jp/labs/nakano/tokuron2/sec61(08).pdf) 2019.09.04閲覧 

  9. 岐阜大学准教授 田中雅宏(2006) 講義配布資料のまとめ (https://www1.gifu-u.ac.jp/~tanaka/numerical_analysis.pdf) 2019.09.04閲覧 

  10. 北川真帆,村上洋一.数理解析研究所講究録(2012), 1800:226-235,京都大学数理解析研究所,https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/173013 

  11. San,O.,Staples,A.E.:High-order methods for decaying two-dimensional homogeneous isotropic turbulence. Comput. Fluids 63, 105–127 (2012) 

  12. 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) 等 

  13. FFTW (http://www.fftw.org/)  

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
31