数理解析研究所講究録
921 巻 1995 年 152-162

152

円形噴流の 3 次元離散渦法数値解析
北大工

井門康司

(Yasushi Ido)

北大工

木谷勝(Masaru Kiya)

1. はじめに離散渦法は、渦度のある領域を渦要素の分布によって表し、各渦要素の位置と渦度を時間の関数として計算するものである。 3 次元離散渦法

について現在様々な試みがなされているが、粘性効果をどのように取り入れるか、計算の進行に伴う渦要素の再配分による渦要素数の急増をどう抑えるかなど、 いくつかの問題点が残っている。 粘性効果の導入方法としてこれまで Random Walk$\text{つ渦度に上限を設ける方法^{}(3})\text{、}$$\text{法^{}(1)_{\backslash }}$Core Spreading$\text{法^{}(2)}\text{、}$

渦要素の持渦度方程式に粘性項の積分表示式を加える方法 (4) などが提案されている。 これらの手法のうち、 Random Walk 法

および粘性効果の積分表示式による方法では、多数の渦要素が計算領域に存在しないとその効果はあまり期待できない。 そこで渦要素の渦度に上限を設ける方法と Core Spreading 法とについて、 それぞれが計算結果にどのような影響を与えるかについて検討を行った。 ここでは円形噴流を取り上げ、 3 次元離散渦法とパネル法を用いて解析した。

2. 3 次元離散渦法と計算モデル
渦要素は円筒形をしており、$\alpha$番目の渦要素の中心位置を$x^{\alpha_{\text{、}}}$

渦度

153

を \mbox{\boldmath $\omega$}\alpha 、体積を$d^{3}x^{\alpha}$とする。渦要素の強さをの渦要素のまわりの渦度分布は\mbox{\boldmath $\omega$}\alpha \Leftarrow )$\gamma^{\alpha}=\omega^{\alpha}d^{3}x^{\alpha}$で表すと、$=p(|x-x^{\alpha}|/\sigma_{\alpha})\gamma^{\alpha}/\sigma_{\alpha}^{2}$

こで表される。 ただし $p(\xi)$ は Winckelmans-Leonard によって提案された球対称関数で $p(\xi)=8/15\pi(\xi^{2}+1)^{7/2}$ である (4)。各渦要素の位置および渦度の時間的変化率はそれぞれ Biot-Savart の法則(1)$\frac{dx^{\alpha}}{dt}=-\frac{1}{4\pi}\sum_{\beta}\frac{r_{\alpha\rho^{2}}+(5/2)\sigma_{\beta}}{(r_{\alpha\rho^{2}\beta}+\sigma^{2})^{5/2}}r\mathrm{x}\gamma^{\beta}\alpha\rho$

および渦度方程式$\frac{d\gamma^{\alpha}}{dt}=\frac{1}{4\pi}\sum_{\beta}1\frac{\gamma^{\beta}\cross\gamma[\alpha r+\alpha\rho 2(5/2)\sigma\rho]2}{(r_{\alpha\beta^{2}\beta^{2}}+\sigma)^{5}/2}$(2)$- \frac{3\gamma^{\rho\alpha\beta}\mathrm{X}r(\gamma\cdot r)\alpha\alpha\beta 1r\rho+\alpha(27/2)\sigma\rho]2}{(r_{\alpha\beta}+22)^{7}\sigma_{\beta}/2}\}$から計算される。 ただし、$r^{\alpha\beta\alpha}=x-X^{\beta}\text{、}r_{\alpha}\rho=|r^{\alpha\beta}|$である。 また時間

積分は 2 次精度の Adams-Bashforth 法を用いて計算した。

図 1 に計算モデルおよび座標系を示した。 ノズル半径を$R_{\text{、}}$

ノズル長さを $L(=0.262R)$ とする。 円形ノズルの壁面境界は周方向に 48 分割

軸方向に 2 分割の計 96 個の四角形パネルを用いて構成した。 さらにノズル内部に 3791 個の吹き出しを配置して、 ノズル出口外縁から渦輪を時々刻々放出して混合層を表現した。放出渦要素はノズルの先端から距
離 $e=0.061R$ だけ離れた位置に導入され、パネル渦要素および放出渦
.

.

. .

要素の循環は、 パネル上の参照点での法線方向速度成分がゼロの条件と Kelvin の循環保存定理とから決定した。 計算は $x=-L$ の面に対して$x<-L$ の領域に鏡像を仮定し、 $x\geq-L$ の領域のみを計算対象とした。

154

また、計算の空間分解能を維持するため、渦要素の長さが初期長さの 2倍以上となった渦要素についてはこれを 2 つの渦要素に分割した。

3. 計算結果と検討
3-1 渦要素の渦度の最大値制限
この計算では図 1 中の噴流出口流速$U_{1}$に対して外部流れ $U_{2}=U_{1}/5$を与え、速度比 $\lambda=(U_{1}-U_{2})/U_{m}$ なる複合噴流としている。 ただし平均流速 $U_{m}=(U_{1}+U_{2})/2$ である。 また渦要素を移動させる時間間隔は\Delta t*=Um ムヴ R $=0.06$ で、放出渦要素を導入する時間間隔は$\Delta t_{v}^{*}=2\triangle t^{*}=0.12$である。 アスタリスクは無次元時間であることを表す。 ここでは渦要素の面面の最大値$\omega_{\max}$

を初期値$\omega_{0}$の 4 倍、 8 倍に制限した場合および制限を設けない場合についてそれぞれ計算を行い、 結果の比較検討を行う。

まず渦要素の上瓦に制限を設けることの物理的意味を考えてみる。 いま粘性流体 (動粘性係数歪速度$\alpha$$\nu$) 中で–様な渦度$\omega_{0}$を持つ半径$\sigma_{0}$の渦管がの軸対称伸長流れによって引き伸ばされたとき、渦度は中心で有限な最大値$\omega_{\max}\alpha\omega_{0}\alpha\sigma_{0}2/\nu \text{をもつ}(5)$。このことは、非粘性流体の場

合 $(\nu=0)$ には心血が発散し、粘性流体では渦伸長による渦度の増加と粘性による半径方向への渦度の拡散とが釣り合って渦度が有限になることを意味している。 また、 Chen ら (6) は差分法によって$Re=80\mathrm{o}\mathrm{o}$での円形噴流の計算を行い、時間進行に伴う渦度の最大値の変化を示しているが、 その値は初期値の 4 倍を越えていない。 これらのことから渦要素の155渦度に最大値制限を与えることは粘性効果を模擬する –つの方法であると考えられる。

図 2 は渦度の最大値制限を初期値のそれぞれ 4 倍、 8 倍、無制限とした場合の時間進行に伴う渦要素数の変化を示したものである。 図より明らかなように、最大値制限を設けることによって渦要素の分割による渦要素数の急激な増加が抑えられている。 これにより、長時間にわたる渦構造の発展を計算することが可能となる。

図 3 および図 4 に上記のような最大感度制限を与えた場合の渦要素分布および感度分布をそれぞれ示した。図 3(a) および図 4(a) に見られる
ように $t^{*}=8.4$ では渦度の最大値制限によらず流れの軸対称性がほぼ保たれており、 渦要素分布、渦度分布ともに大きな差異はない。 –方、 時
刻 $t^{*}=16.2$ の場合には図 3(b) および図 4(b) に見られるように流れの軸対称性はくずれ、 3 次元性が顕著となっている。 しかしながら、 この場合にも渦度の最大値制限によらず比較的大きな渦構造に差異はない。
すなわち、 円形ノズルから出た渦輪はしばらく軸対称性を保っているが、$x/2R=$. $1.5$ の近傍で渦輪の合体が起こる。 さらに下流では図 5 の y-z 断面における渦虫の$x$

方向成分分布図に見られるように縦渦が発生している。 また、 $x/2R>6.0$ の領域では 3 次元性が急速に現れて乱流化が進み、渦要素の分割が頻繁に生じている。 このため渦度の最大値制限を設けた場合の方が無制限の場合よりも分割回数が少なく、尺度の最大値が無制限の場合に見られるような渦要素数の急激な増加が抑えられている。

156

これらの比較的大きな渦構造の定性的な発展は、渦度の最大値制限を設けても無制限の場合と比べてほとんど相違は見られない。 しかし粘性効果の導入方法としては、歪速度は局所的に決まるものであるため、 渦要素の最大値を局所的に決める必要があるなど、 問題点が残されている。

3.2 Core Spreading 法による計算
つぎに粘性効果の導入法として Core Spreading 法採用した場合の計算結果を示す。 本節での計算結果は、 図 1 で外部流速 $U_{2}=0$ としたものであり、 タイムスチップは \triangle t*=0.03、渦放出時間間隔はである。導入された渦要素の渦核半径時間的に値を変化させる。 ただし、$\nu$$\sigma$$\triangle t_{v^{=}}^{*}2\triangle t^{*}=0.06$

は d\mbox{\boldmath $\sigma$}/dt=2\nu \epsilon /\mbox{\boldmath $\sigma$}にしたがっては動粘性係数で、定数$\in$は本計算

では 1 とした。

図 6 および図 7 はそれぞれレイノルズ数が $Re=U_{1}R/\nu=1000_{\text{、}}5000_{\text{、}}$

10000 での渦要素分布および渦度分布を示している。 また、 図 8 は y-z断面における渦度のる。

$Re=10\mathrm{o}\mathrm{o}$$Re=50\mathrm{o}\mathrm{o}$$x$方向成分、すなわち縦渦成分の分布を示していでは流れの軸対称性が比較的保たれているのに対して、

あるいは 10000 では $x/2R=1\sim 1.5$ 付近で渦輪の合体が生じている。 図 6 および図 7 より、 レイノルズ数が 5000 や 10000 の場合に比べ$Re=1\mathrm{o}\mathrm{o}\mathrm{o}$では、 ノズル直径の 7 倍下流近傍での渦度の大きさが明らかに小さくなっている。 これは Core Spreading によって渦度が拡散し、渦要素の引き伸ばしによる渦度の増加が抑えられていることを示している。 また、 この近傍ではレイノルズ数が 5000 あるいは 10000 の場合に

157は軸対称性がくずれ、 図 8 に見られるように縦渦が発生している。

数値誤差を様々な周波数成分を含む擾乱に対応するものと考えると、計算の進行とともに不安定性を最も助長する周波数成分の擾乱が成長するものと考えられる。 そのため、計算開始愛しばらくは流れの軸対称性が保たれているが、 この成長した擾乱によって軸対称性がくずされる。

この擾乱の影響が上流側に伝播し、渦輪の合体が生じるものと考えられる。 しかし、 レイノルズ数が小さい場合には Core Spreading によって密度が拡散されるため、 この擾乱の影響が敏感には効かなくなるものと考えられる。

4. まとめ本稿では渦要素の渦度の最大値制限と Core Spreading という 2 つの粘性効果導入法について検討した。皇恩の最大値制限は比較的大きな渦構造を定性的に見る場合には有効であり、渦要素数の急増を抑制して長時間計算を可能とする。 しかし、粘性効果の導入方法としてはまだ検討を要する。

–方、 Core Spreading 法は粘性拡散を表現するのに適しているが、長時間計算を行い、平均流速などの統計量を得るためには、渦要素数の急増を抑える方法を別途導入する必要がある。今後、 3 次元離散渦法の検証を行うため、長時間にわたる渦構造の発展過程を計算し、平均流速、乱れ強さ、 レイノルズ応力などを実験結果と比較することが課題である。

158

参考文献

(1) Chorin, AJ., J. Fluid Mech., Vo1.57 (1973), 785.
(2) Shirayama, S., Kuwahara, K. and Mendez, R., AIAA Paper,85-1488 (1985).
(3) 木谷・井門・太田, 日本機械学会論文集, 61-583(4) Winckelmans,$\mathrm{G}.\mathrm{S}$$\mathrm{B}(1995)$.

. and Leonard, A., J. Comput. Phys., Vol.109(1993), 247.
(5) Batchelor, G.K., An Introduction to Fluid Dynamics, (1967), 272,Cambridge University Press.
(6) Chen, J.H., Lienau,$\mathrm{J}.\mathrm{J}$. and Kollmann, W., Proc. 9th Symp.

Turbulent Shear Flows, VolJ (1993), 9-4-1.

159

図1

計算モデルと座標系$t^{*}$

図2

時間進行に伴う渦要素数の変化$/\mathrm{h}_{-}1)$$[|||(\mathrm{a}- 2)||||||||||||||||_{||}!_{\mathrm{i}\mathrm{i}1},\cdot|.,i|||i_{\mathrm{t}}(.’$””’”,$|‘$$||||:\mathrm{I}_{;,:}|;_{\mathrm{i}\mathrm{i}}|;||_{\mathrm{i}}^{!}|.|_{:}^{:}\mathrm{i}:.|::\mathrm{I}|‘’ \mathrm{i}|\mathrm{I}||_{\mathrm{i}}|^{1}|_{\mathrm{i}_{\mathfrak{l}}}|\mathrm{t}l|\iota \mathrm{i}||\mathrm{i}$$;\ovalbox{\tt\small REJECT}_{\rangle}|‘,||\#^{\mathfrak{g}_{1}}|;1||^{\mathrm{I}}|_{||}^{\mathrm{t}}\mathrm{I}\mathfrak{i}\mathrm{i}’|$$.\mathrm{i}\mathrm{I}|1|_{\backslash ,;1_{1}}|;||’|_{\mathrm{i}_{;}^{\mathrm{I}}}1\downarrow$$;||^{;}|!\dagger \mathrm{i}^{;}\mathrm{i};|:_{\mathrm{i}}‘$

“”

”

.

$[||||||(\mathrm{a}-3)||||_{1}||||\mathrm{i}||‘|_{\mathrm{i}’}^{|||}|$”’’”$|_{1}^{1}qi’|;\mathfrak{i}\mathrm{i}|\mathfrak{l}|;|||||\mathrm{i}b|,4?|;||^{!}|\mathrm{I}|||\mathrm{T}||_{1},||$$\ovalbox{\tt\small REJECT}_{-1.\overline{0}^{\backslash }}^{0}1.5\wedge\sim\approx \mathrm{a}$$:||;|;;;\iota_{11_{1}^{:_{\mathrm{i}}1}}^{1_{1}^{1}}|i^{\mathfrak{i}}\mathrm{i}^{;}!||^{1}\mathrm{i}\mathrm{I}_{;}!|;;1;‘||_{\mathrm{i}}^{1};_{i:_{\mathrm{t}_{;}},\mathrm{t}1}‘,’.!_{}i|,:||1|_{1}^{1};|.\backslash |:,:!|;:;^{l}|\mathrm{j}_{1};:\mathrm{i}^{\mathrm{I}}||_{\mathrm{i}}^{\mathrm{I}}||^{1}|_{b}||,|||,$$\mathrm{i}|_{||}$

$0$

..

24

2

$0$

$x/(2R)$

図3648$x/(2R)$

最大渦度制限の違いによる渦要素分布の変化(a) $t^{*}=8.4,$$(\mathrm{b})16.2$(1)$\omega_{\max}/\omega_{0^{=}}4.\mathrm{o},$(2) $8.0,$ (3)$\infty$.

160$\mathrm{x}_{\circ\dot{\mathfrak{D}}}^{\S}\sigma_{\mathrm{o}_{\mathrm{O}}}0$$pPo_{\mathfrak{G}\Im}\circ$$\circ\Im^{\Phi_{\mathrm{O}}}\mathrm{b}^{\circ}\mathrm{O}\circ$.

$)$$a_{\sim}$$\circ$$\circ$$.\ovalbox{\tt\small REJECT}$

15$-15$.

th-ti$\ovalbox{\tt\small REJECT}arrow \mathrm{k}_{\iota^{\backslash }}^{\mathrm{Q}}$。

‘$\overline{\approx}$$\circ.\mathrm{o}\mathrm{o}\mathrm{o}$$\wedge\approx$

$\{$

$0$

$\{$$0$$\backslash \vee 0_{\vee}\neg$珂$\ovalbox{\tt\small REJECT}_{\tilde{\infty}_{\mathrm{O}}\ovalbox{\tt\small REJECT}}\ovalbox{\tt\small REJECT}^{a}$15$-1.5$$0$

2$x/(2R)$

4$0$

286

4$x/(2R)$

図4
感度の等高線図( -y 断面) 太線は正の$x$

渦度、 細線は負の渦度を表す。

(a) $t^{*}=8.4,$(1)(2) 8.0, (3)$(\mathrm{b})16.2$$\omega_{\max}/\omega_{0}=4.0$$\infty$,.

図5Pz 断面における渦度の 方向成分の$x$分布$(t^{*}=16.2)$(a),(b) 8.0, (c) ,(1) $x/2R=3.1$ ,(2) 5.1, (3) 7.0.
$\omega_{\max}/\omega_{0}=4.\mathrm{o}$$\infty$161(1)$(’\mathrm{C}^{-}\ddot{\grave{p}}$$.-=\urcorner^{\circ}0\Im C\circ$.

$\cap$$\{$$\backslash$$\ovalbox{\tt\small REJECT}^{0_{\mathrm{O}}}$\copyright\supset$.\mathrm{O}a\wedge$$0$$\mathrm{o}_{O^{o}}\mathrm{O}$$0\mathit{0}$(2)$(_{\tilde{\bigwedge_{\mathrm{q}}^{(-}}}..\overline{\circ}.\mathrm{d}|\overline{\supset}\acute{\backslash }$$-\underline{\wedge}^{---}---arrow \mathrm{g}\copyright^{\sim}\Theta \mathit{0}\ovalbox{\tt\small REJECT}$$[egg0]_{[egg0]}\mathrm{q}\gamma$$\{$<)$\mathrm{i}_{\text{ノ^{}-}}\ovalbox{\tt\small REJECT}\circ- c\approx\circ’ oO\circ$\copyright.

$\mathrm{O}\circ$$\zeta[egg0]\supset$$\mathrm{o}$(3)$— arrow\theta\Phi 0\frac{()}{}\cup \mathrm{c}^{\backslash }’\epsilon 8;3$$\mathit{7}_{\vee}^{J\mathfrak{r}},$.

$.0;^{\llcorner\prime}.,\cdot-\backslash .:’\circ 3\prime^{-}-_{J}\backslash \mathrm{Q}$$\{$$\ovalbox{\tt\small REJECT}_{0}^{\circ}\frac{\mathrm{o}}{-}\dot{\mathrm{o}}0_{\mathrm{r}}^{\alpha}arrow\dot{\sim}.\cdot.\succ\sim$$\grave{\text{ノ}}\ovalbox{\tt\small REJECT}$

D@\approx- も (\copyright \supset$.\ell_{\mathit{0}}7$,$0_{\dot{\mathrm{G}}^{0-}}^{\circ \mathrm{Q}_{\mathrm{Q}}}.s_{\text{。}}\circ\dot{\grave{\check{\dot{o}}}}$$0^{\frac{\circ 1}{.01.02.03.04.0\chi_{2}’\mathrm{R}5.06.07.08.09}}.’.0$

図6

渦要素分布 $(t^{*}=16.2)$(1) $Re=1000,$ (2) 5000, (3) 10000.

図7

渦度の等高線図 $(t^{*}=16.2)$(1) $Re=1000,$ (2) 5000, (3) 10000.

162

(a-3)

(a-1)

(b-2)

(b-1)

(b-3)

驚.

(c-1)

$\mathfrak{s}_{\dot{\oplus}}^{\mathrm{v}^{\mathrm{Q}^{\otimes.\epsilon}}}0\mathrm{t}\dot{\ovalbox{\tt\small REJECT}}.\ovalbox{\tt\small REJECT}\prime 0\mathrm{a}_{\mathrm{b}}\ovalbox{\tt\small REJECT}_{o^{9}}^{\mathrm{e}^{l}}\mathrm{t}^{\emptyset}\delta.\dot{\S}^{\mathit{0}}r_{\ovalbox{\tt\small REJECT}}^{\Phi},$.

(c-2)$\mathrm{q}_{\mathrm{i}}\circ 2.0\exists_{1}\vee$.
’.

$\sim:$$\dot{\circ}.p^{\mathrm{o}_{\tilde{\theta}}^{\circ}}\mathrm{a}_{\approx}^{\partial_{\searrow}\cdot\Im}...\mathit{0}^{\mathrm{e}_{l^{\circ}}}\dot{\mathrm{Q}^{\dot{6}}.}.\cdot..\cdot.\cdot$$- 2.0\aleph-2i.\mathrm{L}_{y/2R}0^{\cdot}\overline{2.}0\approx$

図8Pz 断面における渦度の 方向成分の分布 $(t^{*}=16.2)$(a) $Re=1000,$
5000, (c) 10000,(1) $x/2R=3.0,$ (2) $5.0,$ (3) $7.0$ .
$x$$(\mathrm{b})$