数理解析研究所講究録
第 1822 巻 2013 年 24-31

24

低マッハ数の圧縮性乱流と非圧縮性乱流について

名古屋工業大学大学院創成シミュレーション工学専攻
中嶋大樹 (Taiki NAKAJIMA)
畑中祥吾 (Shogo HATANAKA)

渡邊威 (Takeshi WATANABE)
後藤俊幸 (Toshiyuki GOTOH)

Nagoya Institute of Technology

1

はじめに
非圧縮性乱流の数値的研究において，空間，時間の解像度が大変重要であることは良く知られている．コルモゴロフ長さと時間はそれぞれ $(Re^{-z}, Re^{-2})31$ で小さくなるので，高レイノルズ数乱流の直接数値計算では，大きな計算資源を必要とする．このため，現在の大規模計算においては，並列計算機を用いることが一般的になっている．非圧縮性流体の場合，ボアソン方程式を解く必要があり，この解法として代表的な手法には差分法，スペクトル法などが挙げられる．差分法では反復解法を用いてボアソン方程式を解くことが多いが，収束するまでのデータが無駄になってしまい，精度もスペクトル法に比べると悪くなってしまう．それに対し，スペクトル法では精度は良いがノード間のグローバルな通信が必要となり大規模計算には不向きである．圧縮性乱流についてもスペクトル法による計算が木田，三浦により行われており [1], 定常乱流における圧縮成分の運動エネルギーと内部エネルギー，圧力勾配それぞれのスペクトルの相互関係について述べられているが，空間解像度は大きくなく通信の問題はなかった．圧縮性流体では，全ての変数が陽に時間発展の方程式で与えられるため，ボアソン方程式を解く必要が無くなる．そのため，計算手法としては差分法が主に用いられる．差分法では空間的にローカルで解くことができるため通信量が抑えられ大規模計算に向いている．だが，やはり精度の面ではスペクトル法に比べると劣ってしまう．そこで差分法として，計算負荷は低いが精度の高い結合コンパクトスキームを用いることが一つの方法として考えられる．流体の圧縮性の効果はマツハ

数の 2 乗に比例して入ってくるので，マッハ数をを小さく設定すれば，実質的に非圧縮性流体を扱うことになるのではないかと考えられる．もしそうなら，通信負荷の低い計算で非圧縮性流体の数値計算が可能となるので，本研究ではどの程度のマツハ数までなら非圧縮性流体の性質を持つのか又，圧縮性乱流そのものの統計性質はどの様なものなのかを探る．

252問題の定式化
密度 , 速度 , 圧力 , ストレステンソル$u$$\rho$

ギー ,$\epsilon$$\mathcal{L}$$p$

温度丁，熱伝導率, 代表的な速度を$\mathcal{U}$$\kappa$$\sigma_{ij}$,トータルエネルギー $E$ , 内部エネル, ガス定数 $R$ と表し，流体運動を特徴づける代表的な長さをとして，方程式に用いる変数を以下のように無次元化する．
(1)$x’= \frac{x}{\mathcal{L}}, t’=\frac{t}{\mathcal{L}/\mathcal{U}}, u’=\frac{u}{\mathcal{U}}, \rho’=\frac{\rho}{\rho_{0}}, p’=\frac{p}{\rho_{0}\mathcal{U}^{2}}$上式でダッシュ付きの変数は無次元量を表す．これを用いて，連続方程式，運動方程式，エネルギー方程式，状態方程式，トータルエネルギーはそれぞれ以下のようになる [2]. ここで繰り返された添字については和を取るものとする．ただしダッシュは省略している．
$\frac{\partial\rho}{\partial t}+\frac{\partial\rho u_{i}}{\partial x_{i}}=0$,(2)$\rho(\frac{\partial u_{i}}{\partial t}+\frac{\partial u_{j}u_{i}}{\partial x_{j}})=-\frac{1}{\gamma M^{2}}\frac{\partial p}{\partial x_{i}}+\frac{1\partial\sigma_{ij}}{Re\partial x_{j}}$,(3)$\frac{\partial E}{\partial t}+\frac{\partial(E+p/\gamma M^{2})u_{i}}{\partial x_{i}}=\frac{1}{\alpha}\frac{\partial}{\partial x_{i}}(k\frac{\partial T}{\partial x_{i}})+\frac{1}{Re}\frac{\partial\sigma_{ij}u_{j}}{\partial x_{i}}$$p=\rho T$

,

(4)

,

(5)

(6)

$E= \frac{p}{(\gamma-1)\gamma M^{2}}+\frac{1}{2}\rho(u_{j}u_{j})$となる．ただし，Re(レイノルズ数), Pr(プラントル数), $M$ (マッハ数),$Re= \frac{\rho_{0}\mathcal{L}\mathcal{U}}{\mu},$$M= \frac{\mathcal{U}^{2}}{c^{2}},$$Pr= \frac{\mu}{\rho_{0}\kappa},$$\alpha$はそれぞれ，$\alpha=PrRe(\gamma-1)M^{2}$(7)と表される．境界条件は，周期境界条件を用いる．初期条件については，非圧縮でランダムな波数空間での速度場を与える．初期エネルギースペクトルは，$E(k)=16 \sqrt{\frac{2}{\pi}}u_{0}^{2}k_{0}^{-5}k^{4}exp(-2(\frac{k}{k_{0}})^{2}), k_{0}=6, u_{0}=1$(S)$\frac{3}{2}u_{0}^{2}=\int_{0}^{\infty}E(k)dk$(9)である．これで得た $u(k, 0)$ を FFT を用いて実数空間へ移す． は初期条件のみ非圧縮条件であるのですべて 1 としているが， を決める必要がある．圧力のボアソン方程$\rho$$p,$$\epsilon$

式は，$\frac{\partial^{2}p}{\partial x_{i^{2}}}=-\gamma M^{2}(\rho\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{j}}{\partial x_{i}}-\frac{1}{Re}\frac{\partial^{2}\sigma_{ij}}{\partial x_{j}\partial x_{i}})$(10)26となる．これを $SOR$ 法で解いたのち，式 (6) にその値を代入し， を計算している．こ$p_{0}=1$ を加えて $p=1+p$ として式 (6) にのは，ゆらぎ成分であり，静力学分の圧力，$\epsilon$$p$代入する．

3 数値計算法
一般に差分法において，精度を高くするためにはそれだけ多くの参照点が必要になる．
そのため，精度が高くなればなるほど計算時間は長くなってしまう．それに対し，コンパクトスキームでは，一階微分値を陰的に解くことにより，少ない参照点で高精度な解が得られる．今回はそれを更に発展させた結合コンパクトスキームを用いる [3,4].

3.1

コンパクトスキーム

まず，コンパクトスキームでは以下の式を用いる．
$\beta f_{1-2}’+\alpha f_{1-1}’+f_{i}’+\alpha f_{i+1}’+\beta f_{1+2}’$$=c \frac{f_{1+3}-f_{i-3}}{6h}+b\frac{f_{i+2}-f_{i-2}}{4h}+a\frac{f_{i+1}-f_{i-1}}{2h}$

,

(11)

$\beta f_{i-2}"+\alpha f_{i-1}"+f_{i}"+\alpha f_{+1}"+\beta f_{i+2}"$$=c \frac{f_{i+3}-2f_{i}+f_{1-3}}{9h^{2}}+b\frac{f_{i+2}-2f_{\dot{\iota}}+f_{i-2}}{4h^{2}}+a\frac{f_{i+1}-2f_{1}+f_{i-1}}{h^{2}}$(12)

式 (11), (12) がそれぞれが一階微分，二階微分を求める際に用いる式である．上式の両を中心としたテイラー展開を行い， の乗数が同じ項についての係数比較を行うことで，各係数を求める．ここでどこまでの項を考慮するかによって，精度が決まって辺に点$h$$i$くる．

3.2 結合コンパクトスキーム
コンパクトスキームを発展させた手法が結合コンパクトスキームである．このコンパク
トスキームとの大きな違いは，一階微分値に加え二階微分値も用いることである．このことで更に精度が上がり，また，一階微分，二階微分が一度で得られる．本研究では，8 次精度の結合コンパクトスキームを用いる．実際に使う式は以下のようになる．

$180f_{i}’+51(f_{i-1}’+f_{i+1}’)+9h(f_{i-1}"+f_{i+1}")$$= \frac{107}{h}(f_{i+1}-f_{i-1})-\frac{1}{h}(f_{i+2}-f_{i-2})$

,

(13)

$138(f_{\dot{|}+1}’-f_{i-1}’)-h(18f_{\dot{\iota}-1}"-108f_{i}"+18f_{i+1}")$$=- \frac{702}{h}f_{i}+\frac{352}{h}(f_{1+1}+f_{1-1})-\frac{1}{h}(f_{i+2}+f_{i-2})$

(14)

27
式 (13) を一階微分，式 (14) を二階微分にそれぞれ用いる．本研究ではこれらの式を以下のブロック行列を並列計算で解く．各プロセッサがそれぞれの行を担当し計算する．

$L_{2}C_{1}$$A_{j}C_{2}$$L_{N-1 ,A_{N}}$$L_{j}$$\ovalbox{\tt\small REJECT} C_{N}A_{2}L_{1}$$C_{N-1}A_{1}]\{\begin{array}{l}X_{1}X_{2}|X_{j}|X_{N-1}X_{N}\end{array}\}=[D_{N-1 ,D_{N}}D_{j}D_{2}D_{1}\ovalbox{\tt\small REJECT}$

,

(15)

$A_{N-1}^{C_{j}}$$L_{j}=\{\begin{array}{llll}a_{2}^{j}\dot{\nu}_{1} c\dot{i}b_{2} c_{2} \ddots a_{n}^{j} b_{n}\end{array}\},$$A_{j}=\{\begin{array}{llll}0 \cdots 0 a_{1}^{j}0 \cdots 0 0| |0 0\end{array}\},$ $C_{j}=\{\begin{array}{llll}0 \cdots 0 0| \cdots 0 00 .

|c_{\dot{n}}^{i} 0 \cdots 0\end{array}\}$4 計算結果3 つの Mach 数について減衰乱流の計算を行った．計算パラメータの詳細を表 1 に示

す．小さいスケールについての非圧縮成分の目安としてはエンストロフィー , そして圧縮成分については $divu$ が目安となるので，$divu$ の分散や規格化された 3 次モーメン$\Omega$$\Theta$

トをそれぞれ調べた．
$\Omega=<(rotu)^{2}>, \Theta=<(divu)^{2}>, S_{divu}=\frac{<(divu)^{3}>}{<(divu)^{2}>^{\frac{3}{2}}}$

ここで $<>$ は体積平均である．運動エネルギーの時間発展を非圧縮成分，圧縮成分に分けて図 1 に示す．非圧縮成分については，どのマッハ数においてもほぼ同じ値を示しており，時間と共に減衰している．それに対し，圧縮成分ではマッハ数が高くなるほど，運動エネルギーも大きくなっていることが分かる．減衰の仕方はマッハ数が低いほど振動の振幅が強くなっている様子が見られる．の時間発展を図 2 に示す．これから分かるように，マッハ数の違いによる差は見られなかった．初期条件で非圧縮性としているため，回転運動が支配的にエンストロフィー$\Omega$起きている．そのため，初期では渦度の値は大きくなっており，時間の経過と共に減衰している様子が分かる．

乱流レイノルズ数の時間発展を図 3 に示す．これにも渦度と似た様子が見られ，マッハ数の違いによる差はほとんどなかった．初期で大きな値を示しており，時間の経過と共に減衰している．マッハ数 0.05, .1 と 0.3 の間でわずかに違いが見られ，マッハ数が高い$O$ものの方が$R_{\lambda}$は低くなっていた．

エネルギースペクトルの比較を図 4,5 に示す．非圧縮成分については，マッハ数が高28

表1パラメータ

い方が減衰が速くなっており，その傾向は高波数側でより強く出ている．圧縮成分では，マッハ数が高い条件の方が全体的に大きな値を示している．
次に $divu$ の二乗平の時間発展を図 6 に示す．初期の値の上昇については $Ma=0.05,0.1$$Ma=0.3$ はかなり大きな値を示している．その後の値の減衰の仕方についてと比べると，

も $Ma=0.05,0.1$ においては$0$$Ma=0.3$ では，有限の値にとどにかなり収束しているが，まっている．
$divu$ の歪度
続いて，$S_{d1vu}$$Ma=0.05,0.1$ についてだの時間発展を図 7 に示す．まず，が，初期で正の方向，つまり膨張の傾向が見られ，その後は の周りで細かく揺れ動く様子を示している．それに対し $Ma=0.3$ では初期に負の値，つまり圧縮の傾向が見られて$O$$k_{0}u_{0}t=1$ 近傍では初期の様子よりも強く圧縮の運動が起こっており，その後はいる．

$0$の周りで振幅している．また，変動の特性時間についても違いが出ており，マツハ数が高くなるほど，時間変化がゆつくりとなる傾向が見られた．
マッハ数 0.3 での $divu$ の吻 $k_{0}t=1.2$ においての様子を実際に可視化したものを図 8,

図 9 に示す．図 8 が膨張の様子，図 9 が圧縮の様子である．図 8 では，divu $<0$ の部分を赤色で示し，$divu>0$ の部分を絶対値が大きくなるに連れて，色が緑から青に変化す$divu<0$ の部分を絶対値がるように描いた．図 9 では，$divu>0$ の部分を青色で示し，大きくなるに連れて，色が緑から赤に変化するように描いた．これを見ると，図 8 では膨張を強く示す部分が広く分布しており，図 9 では圧縮を強く示す部分は狭い範囲に密集していることが分かる．図 7 では uokot $=1.2$ の時に圧縮の傾向が強く見られているため，圧縮が起きている範囲は狭いが，その強さは大きいことが分かる．
$divu$ の確率密度分布を規格化せずに示したものが図 10,ては，どの時間どの空間でも $divu$ は$0$

$Ma=0.05$ におい
11 である．に近い値を示しており，非圧縮の性質を示してい$Ma=0.3$ では初期に裾野が少し開き，左に広がった非対称な形をしている．それに対し，る．これは図 7 において初期に$S_{divu}$が負の値を示していることや，図 8, 9 において圧縮成分が強く出ていることに対応している．そして，その後徐々に減衰している．
$\rho$$divu$ の結の揺らぎの確率密度分布を規格化せずに図 12,13 に示す．これについては，$Ma=0.05$ では全ての時間において果と似た傾向が出ている．

については初期に裾野が開き，その後減衰していく．

$0$$Ma=0.3$
を示しており，

29

ル $a=0.3E$ in$Ma=0.1E$ in$Ma=0.05E$ in–––$Ma=0.3Ecom-$$Ma=0.1Ecom-$$Ma=0.05Ecom-$$Ma=0.3$ Total Energy –$Ma=0.l$ Total Energy –$Ma=0.05$ Total Energy –$0 1 2 \mathcal{U}d\zeta 0^{t}3 4 5 6$
図1

図2

図4

運動エネルギーの時間発展

渦度の時間発展

エネルギースペクトル $(Ma=0.05)$

図3

図5

乱流レイノルズ数の時間発展

エネルギースペクトル $(Ma=0.3)$

30

図

図8$6$$divu$$<divu^{2}>$ の時間発展の膨張の様子 (uokot $=1.2$ )

図9

図7$divu$ の歪度の時間発展$divu$の圧縮の様子 (uokot $=1.2$ )$divu<0$ の部分を$(divu>0$ の部分を青，絶対値が大きくなるに連れて，緑から赤)$divu>0$ の部分を
( $divu<0$ の部分を赤，絶対値が大きくなるに連れて，緑から青)$J0^{0}$$J0^{0}$$J0^{J}$$l0^{-1}$$\tilde{Q_{\backslash }}J0^{-3}\wedge:_{o}^{d}>\approx J0^{2}$$\tilde{Q_{\backslash }}10^{\vee}\sim_{l}:_{o}^{\mathcal{T}\dashv}\approx 0_{3}^{-2}>$$10^{-4}$$J0^{-4}$$10^{-5}$$10^{-5}$$-3$$-2$$-,$$0$

123$-3$

図 10$divu$ の確率密度分布 $(Ma=0.05)$

$-2$

$-1$

$0$

1

2$divu$$divu$

図 11$divu$ の確率密度分布 $(Ma=0.3)$331$p$

図 125$\rho$$p$のゆらぎの確率密度分布 $(Ma=0.05)$

図 13$\rho$のゆらぎの確率密度分布 $(Ma=0.3)$

まとめ

どの値についてもマッハ数を小さくするほど圧縮の傾向の現出が小さくなっており，予想された振る舞いを示している．運動エネルギーについては，マッハ数が低い条件ほど圧縮成分におけるエネルギーの揺らぎが大きく，圧力項とのエネルギー推移が起こっていることが分かった．波数空間での圧縮成分，非圧縮成分それぞれのエネルギースペクトルでは，マッハ数が低い場合，非圧縮成分と圧縮成分の間に 4 桁ほどの違いがあり，圧縮$divu$ については，限りなく成分は無視できる値であると言える．に近い値をとっている．どの統計量についてもマッハ数が .1 以下では，圧縮成分は無視できるほど小さい値を示すことが分かった．今後の展望としては，本研究の計算と非圧縮性流体の計算を実際に比較し，どれだけ非圧縮性流体の振る舞いを再現できているかの検証や，スペクトル法での圧縮性流体計算との比較による，結合コンパクトスキームの制度の検証などが挙げら
$0$

$O$

れる．

6 文献
[1] H.Miura and .kida,”Acoustic energy exchange in compressible turbulence” Phys.
$S$

Fluid 71732(1995)[2] R.Samtaney, D.I.Pullin,and B.Kosovic,”Direct numerical simulation of decaying
compressible turbulence and shocklet statistics” Phys. Fluid 131415(2001)[3] K.Mahesh,”A Family of High Order Finite Difference Schemes with Good SpectralResolution” Journal of Computational Phys. 145332(1998)[4] T.Nihei and K.Ishii,”Parallelization of a Highly Accurate Finite Difference Scheme
for Fluid Flow Calculations”Theoretical and Applied Mechanics Japan 5271(2003)