数理解析研究所講究録
第 661 巻 1988 年 53-67

53

3 次精度上流差分による乱流遷移の数値シミュレーション\S 1.

鳥取大学工学部

河村哲也(Tetuya Kawamura)

計算流体力学研究所

岩津玲磨(Reima Iwatsu)

はじめに層流から乱流への遷移には大別して二つの型がある.
ばれるもので,

”

一つは”

突発型と呼

微小擾乱に対して流れは安定であるが十分高いレイノルズ数では乱流となり中間のレイノルズ数では乱流と層流の入り混った状態が存在するというもので,$C$

Po seu lle 流,

円筒内

$i$

$i$

Coue $tte$ 流 (Taylor-

回転同軸円筒内のou $ette$ 流) で外円筒の回転速度が卓越した場合などの乱流遷移がその代表的な例になっている.

”

もう一つはトル遷移型

スペクと呼ばれるもので,”の出現後種々の流れを経て最終的に乱流に移行するもので,$Tay1$微小擾乱

自由せん断流やor-Co $uette$ 流で内円筒の回転速度が卓越した場合の遷移などが代表例である.

一方,二次元$Po$

I $seui11e$ 流や平板境界層の遷移はその中間的なタイプと見ることができる.

本研究の目的は,レー乱流への遷移を三次精度の上流差分を用いてどの程度シュミト可能なのかを調べることであるが,1)遷移に対してはある程度の成果をうることができた

トル遷移型の代表例である$Tay1$

or-Couet $te$ 流は,$T$

Po Ise I11中間型である二次元ので,$u$

本研究ではスペク

ay lor-Cou $ette$ 流の外円筒静止の場合を調べる.

Tay lo$r$の先駆的な研究

2)の後,種々の理論的・実験的な研究が行われてきた流れであるが,最近の電子計算機の進歩により,な研究も増加しつつある.

ose

その中で

ル法を用いた三次元計算を行い,る波をともなう”Wav”
$y$

$T$

$M$

ay lo

$r$$r$

et al3)やMar cu)はスペク トな流れを数値的に得ることに成功している.

計算では一対の$T$

軸方向に基本波数$\alpha$

一方,そに関する周期性をay 1or 渦の挙動を調べていることになっている.

本研究でも軸方向に周期境界条件を課しているが一周期内に$T$

$4$

$s$

数値的渦が不安定となって方位角方向に進行す

れらの論文では無限に長い円筒を仮定し,課しているため,流の$e$ay1or 渦がはいることのできる長さであるため,-1-10 対程度の周期内では渦の大きさが異なっ54
周期性に少し乱れが生じた流れ (図 6 (a) ) や,た流れも取り扱うことができ,周期内での渦の個数の変化も取り扱えるため,より実験条件に近い計算になっている.

\S 2.

基礎方程式と計算方法

基礎方程式は連続の式 (1)

ここで$C$

法5)

V

$N$

av $ie$ r-S $to$ kes 方程式 (2)である:$divV=0$(1)$\partial V/\partial t+(V\cdot\nabla)V=-gradp+\frac{1}{Re}\Delta V$(2)は流速,に従って,と$P$は圧力,$Re$ はレイノルズ数を表す.

(1)(2) を連立させて解くかわりに,式を考慮して得られる圧力に関する

$Po$

$i$

この計算では MA(1)(2) の発散をとってsson 方程式 (3), (2) を連立させて解を求める:(3)$\Delta P=-div(V\cdot\nabla)V+\frac{D}{\Delta t}$

ここで右辺第 2 項のいる時間における項は,$div$$\Delta t$

は(2) を時間積分する時の時間きざみで $D$ は着目してV の差分近似値(したがってゼロではない) である.

(2), (3) 式を差分近似して解いて得られるV が,常にこの(1) を近似していることを保証する補正項となっている.
同軸二重円筒問の流れを問題とするため用いるのが便利である (図 1).

このとき(3), (2) に円柱座標系 ( $r,$$\theta$, z) を(3), (2) 式は次のようになる:$\frac{\partial^{2}P}{\partial Z^{2}}+\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial P}{\partial r})+\frac{1}{r^{2}}\frac{\partial^{2}P}{\partial\theta^{2}}$$=- \frac{\partial}{\partial Z}(Vz\frac{\partial Vz}{\partial Z}+Vr\frac{\partial Vz}{\partial r}+\frac{V\theta}{r}\frac{\partial Vz}{\partial\theta})+\frac{1}{r}\frac{\partial}{\partial r}\{r(Vz\frac{\partial Vr}{\partial Z}+Vr\frac{\partial Vr}{\partial r}+\frac{V\theta}{r}\frac{\partial Vr}{\partial\theta}-\frac{V\theta^{2}}{r})\}$$+ \frac{1}{r}\frac{\partial}{\partial\theta}(Vz\frac{\partial V\theta}{\partial Z}+Vr\frac{\partial V\theta}{\partial r}+\frac{V\theta}{r}\frac{\partial V\theta}{\partial\theta}+\frac{VrV\theta}{r})+\frac{1}{\Delta t}(\frac{\partial Vz}{\partial Z}+\frac{1}{r}\frac{\partial}{\partial r}(rVr)+\frac{1}{r}\frac{\partial V\theta}{\partial\theta})$-2-(3)

55$(2a)’$$\frac{\partial Vz}{\partial t}+Vz\frac{\partial Vz}{\partial Z}+Vr\frac{\partial Vz}{\partial r}+\frac{V\theta}{r}\frac{\partial Vz}{\partial\theta}=-\frac{\partial P}{\partial Z}+\frac{1}{Re}\Delta Vz$$\frac{\partial Vr}{\partial t}+Vz\frac{\partial Vr}{\partial Z}+Vr\frac{\partial Vr}{\partial r}+\frac{V\theta}{r}\frac{\partial Vr}{\partial\theta}-\frac{V\theta^{2}}{r}=-\frac{\partial P}{\partial r}+\frac{1}{Re}(\Delta Vr-\frac{Vr}{r^{2}}-\frac{2}{r^{2}}\frac{\partial V\theta}{\partial\theta})$$(2b)’$$\frac{\partial V\theta}{\partial t}+Vz\frac{\partial V\theta}{\partial Z}+Vr\frac{\partial V\theta}{\partial r}+\frac{V\theta}{r}\frac{\partial V\theta}{\partial\theta}+\frac{VrV\theta}{r}=-\frac{1}{r}\frac{\partial P}{\partial\theta}+\frac{1}{Re}(\Delta V\varphi-\frac{V\theta}{r^{2}}+\frac{2}{r^{2}}\frac{\partial Vr}{\partial\theta})(2c)’$

A を Vz

ただしまたはVr または Vとして$y$$\Delta A=\partial^{2}A/\partial Z^{2}+(1/r)\cross\partial(r\partial A/\partial r)/\partial r+(1/r^{2})\cross\partial^{2}A/\partial\theta^{2}$である.

壁$(r=r_{\Xi} ’ r=r_{b})$の近くでは差分格子を細かくとる必要があるため半径方向に次の変換をおこなう:(4)$r=r_{a}+(r_{b}-r_{a})\cross(1+\frac{e^{x}-1e^{c}+1}{e^{x}+1e^{c}-1})\cross\frac{1}{2}$

ただし$-c$$\leqq$$x$

1 $og$ $(1 +0.95)$
なお$r_{a}$

(3)

$-$

$\leqq$$/$

$c$

で$C$

は格子の集中度をきめるパラメータでたとえば$(1 -0.95)$ととる.

( $0.95$ を1 に近づけばより集中する.

)は内円筒の半径である.

式はすべで 2 次精度申心差分で近似した。 一方 (2)線形項を除き,$-$式は,空間微分は非2 次精度中心差分で近似し非線形項は 3 次精度上流差分$f \frac{\partial U}{\partial X}=f\frac{-U_{j+2}+8(U_{i+1}-U;_{-1})+U_{i-2}}{12\Delta X}+|f|\frac{U_{i+2}-4U_{i+1}+6U;-4U_{i-1}+U_{i-2}}{12\Delta X}(5)$で近似した.
なわちa $u/$

また時間積分は 2 次精度は Adams-Bashforth 法を用いて行った.
$\alpha t=f(u)$に対して$U^{n+1}=U^{n}+\Delta t(3f(u^{n})-f(u^{n-1})/2$-3-(6)す56を用いる.

まず初期条件または前の時間ステップでの速度より (3)解を求める手順は,近似式の右辺を計算し,S.

$0$. R. 法により圧力の Po ss

して得られる (2)$-$の方程式を解き,圧力を求得られた圧力と現在の速度より(6) を参照$i$

次の時間ステップでの速度は,める.

$-$の近似式より求める.

$on$

この手続きを必要な時間ステップ繰り返

時間発展させるときに 2 次精度の Adams-Bas $hf$ or $th$ 法を用いたためす.

なお,

(8)

よりわかるように一つ前の時問ステップでのまた最初のステップではこの方法は使えないため,分を用いた.

格子システムは,$f$を記憶しておく必要がある.
このときに限り Eu 1er 前進差速度と圧力を同一点で評価する通常格子 (図 1)を用いた.
壁面での境界条件は,ては粘着条件を$N$速度に対しては粘着条件 (V$=$回転速度),圧力に対しav er-S tokes 方程式た代入して得られる
$i$

(7)

$\partial P/\partial r=V\theta^{2}/r-1/Re\cross\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial Vr}{\partial r})$を差分近似した式より壁面での圧力を求める.

ただし(7) で 2 階微分を求めるとき仮想点での値を決める必要がおこるので図 1 のように決めた.
精度の差分を用いたため,非線形項に 3 次やはり仮想点を必要とするのでこの場合も図 1 のように決めた.
他の面での境界条件は周期条件を課し,また初期条件としては$C$

ou $ette$ 流の厳密解 (8) を用いた.

$Vr=Vz=0$$V \theta=\frac{(V_{\theta a}-V_{\theta b})r_{a}^{2}r_{b}^{2}1}{r_{b}^{2}-r_{a}^{2}r}+\frac{V_{\theta b}r_{b}^{2}-V_{\theta u}r_{a}^{2}}{r_{b}^{2}-r_{a}^{2}}r$

ただし

Vea,

Ve

$b$

内外円筒の回転速度である.

\S 8. 計算結果

-4-(8)

57
$Tay$

$1$

or-Co $uette$ 流の流れのパターンを決めるパラメータは多数ありすべてを変化させて調べることは現実的でないため,場合のみ取り扱う.

これは,比/$r_{b}$$r_{a}$

ノルズ数
$F$

,

流れの支配パラメータは内円筒と外円筒の半径円筒間のギャ ッ プと円筒の高さ

Voa$Re=\Gamma a$ens $ter$ macher $et$ a1$(r_{b} -r_{a})$
6)

$/$

$v$

:

けた計算を行うことは,$(r_{b} -r_{5})$$/$$h$

動粘性係数)の実験に従い,$/\Gamma e$$\Gamma b$である.

$1.14$ ,$=$$h/$

ただし,

本計算では軸方向($i$ay 1 $or$ 渦がはいるので,(Iv) 経験的事実として,間のかかる圧力の $Pois$ so 方程式の収束が,$n$同じ計算時間が許されているとき,値に比べ$T$23 点である

23 である.

(ただしレイノルズ数が理論的に予測されるは$Re$$\Delta t$は203,$P$

な8 メ ッ シュ内側)

ar $tic$ le Pa $th$476,

238,それぞれ64,

半径方43,45,

1904 の場合を行った.

$T$$Re=110$ , 238 のとき$y$

110 の場合のみ,周方向に

なお

Cr $iti$ ca 1 なレイ ノルズ数 ( ay lo 渦が発生するレイノルズ数)

119 であるので上記の値はTay 1or 渦や Wav面より110,) 計算はレイ ノルズ数が7 》により調べられ,13% ほど高くなることがわかっている.のとき計算に用いた格子数は仮想点まで含めると軸方向に向に$Co1e$無限長とした理論計算の予測

$r$

$=20$

本計算法で最も計算時より多くの格子点をとった計算ができること

ay lo 渦が発生する臨界レイノルズ数が,$(r_{b} -r_{a})$そし周期条件を課した法がはやいため,なお有限長の円筒を用いた影響は実験的に$h=$

殴周期条件を課し

1 対の Tay 1 $or$ 渦のみの挙動を調べる計算よりは実験に近て現実的だが重要な理由として,な周期条件を選$(iii)$ 低面の影響を調べる場合の比較データになること.

くなっていること

avy

底面をつ10 対程度のたため一応無限の長さの計算になっているがひとつの周期内に$W$$=20$$(I^{\cdot}b -r_{\text{\‘{e}}})$

境界条件を少し修正すれば可能であるが,にもよるので完全に対応した計算が困難であることなどである.

本研究では(i) 実験的に知られているように流れのパターンは履歴 (ヒステ

んだ理由として$T$およびレイ実験とは完全な対応はついていない.にも周期条件を課しているため,リシス)($v$の比$h$

レイノルズ数を変化させて計算を行った.に固定し,外円筒静止の理論的にも実験的にも数値的にも最もよく研究されこの場合,ているケースである.

本研究では内円筒回転,Tayloのの$Re$

1/1000$r$$r=$$0.92$ ,2,それ以外は$r$4,16 倍となっている2.5/1000 である.

渦をわかりやすく表示するため,

consを可視化に用いた.

$t$. の画内で

(なお)

外円筒近く

(表

Vr を無視した場合のただし流速はある特定の時間のものを固定し-5-

58
このようにすると図 2 より,て使った.

きこまれる線に沿って粒子が集まり,筋の間に一対の$y1$

or 渦があれば,流れが内側にひ黒い筋ができると考えられる.

Wav

流れがそして筋とになれば筋も波う$y$)(これが実験で撮られる可視化写真に最も近いと思われる.

Vr まで含めて$P$

ar $tic1e$ Pa $th$

存在をよりはっきり示すため,$T1me$

をつ くると定常な

ay 1 $or$ 渦の場合でも,$T$

Wav

2 次元的に表示すると波うって見える.

線軌道を描くため,方向の速度のa$T$

ay lor 渦が存在することになる.

$T$

つと考えられる
なお,もし流れの中に特定の点を指定して,$history$ も合わせ表示した.

Vの値が異なる 8 点でのを表示した.

$r$

本研究では

もし定常な

$T$

$R$

V

$r$

$T$

ay 1 $or$ 渦のその点での半径,は同じで$\theta$

ay 1 $or$ 渦が存在すれば,はそれぞれの 3 点で時間的に変化しないある一定値をとるが,ていたとするとな$y$らもしはある周期的に変化する値をとることになる.

$W$

$z$

Vr

avy になっ

またその図より進行波の速度も計算することができる.
最初にテス ト計算としてより低いレイノルズ数で,$Rc$

この場合に限り初期条件は静止とした.

おこなった.の速度ベク トルを示すが,予想どおり図には示されていないが,Vz格子数も少ない計算を

図 3 に軸に垂直なある面でCo $uette$ 流が実現されている.

Vr もゼロであることや,もこのことは,どの断面でも速度分布が同じであることからもわかる.

図 4 にレイノルズ数 $Re=238$ の場合の計算結果を示す.
$P$

at

$P$

ar $tic$ le Pa $th$ である.

各$h$

$T$

で,図の左側は$\theta$$=c$ons . の面で切った切口で V
$t$

この図より

11 対 (22 個)

ay lor 渦はほぼ同一の大きさであるが,ることがわかる.

図 4(6) はごとに表示したもので,あたりまではV$r$$=0$

りまではわかった

Vr$=0$であるが,$Re=$$T$$=0$$131$であるが,$z$$T=15$

それ以降,ay 1or 渦ができていること,の異なる 8 点での( $T=2\pi$$0$

Vr の値を各時間で内円筒が一回転する)でない値をもち$(= 1.1 {\rm Re})$$T=80$ 程度で一$(T=0)$に内円筒の回転速度の2.5$\chi$の大き) 次にこの計算で同じレイノルズ数であるが,-6-

図にTay lor 渦になることがCoue $tte$ 流と重ね合わせている.

Tay lor 渦流を初期条件にとり,なお,も計算をおこない $T=40$ あたそれ以後値をもつようになり(ただし両ケースとも初期

ar $ti$ cle

中央付近のものが若干大き くなっていが一定での場合$P$としたときの定常な Tay 1or 渦ができていることがわかる.

さの撹乱を一様乱数で与え,得られた$\theta$

この図より定値をもつようになり,は示さなかったが$r,$の$o$

図 4(a) はこれに上と59
同じ

2.5

から

$T=240$ 経た時間での Par $tic1e$ Patる

(なおの撹乱を加えた場合の計算結果を図 5 に示す.

$\chi$

$Re=$

$131$

以降周期的に速度が変化して,

Wav

$y$

$T$

$=$

進行波の進行速度は図より計算でき

$T=$

た

$120$

$T$と$T=$ay lor 渦が$160$

$T$

$=$

Tay 1 $or$ 渦の存在がわかの
$160$

$P$

ar $tic$ le Pa$\theta$

一定で z のになって、、ることを裏付けている.

この$r$から$170$

180 の間でおおよその値を見積
図 5(6) , (7) で$T=$より以$120$

図には示さなかったが,を比較すると$th$$r,$$T=$

160 以降では波形が異なっている.

$T=$

avy な

V の時間変化を示したものででの1/2 程度になっている.

ると内円筒の回転速度の前と$N$の場合は Wavy にならない) 図 5(b), (C) は異なる 3 点 (図 4(6) と同じ点ではない)
$50$

で,

$h$

図 5(a) は撹乱を加えて$T=$$120$のとき11 対あっ9 対 (図 5(a) と同じ) になっているのでこの間でのとき渦の融合と再配列がおこっていることになる.
Tay 1or 渦はいったんできるとかなり安定で,この計算でおこなったように撹乱を加えない限り早い時間内で Wavy にはならない.
あるがCou $ette$ 流を初期条件に用いると,$T$

図 6 (a) は $Re=476$ の結果で.

ay 1or 渦ができ撹乱を加えずに計算を続けると図 6(a) に見られるように周期性にところどころ乱れが出てくるが,以降,渦が融合してまた図 6 (b) に見られるようにTay lo 渦流にもどる傾向がある.
$r$

しかし撹乱を加えるとWavy な流れとなる.

さらに高いレイ ノルズ数の計算をおこなうためには,格子数を増加させる必要があると思われるが定性的な議論だけなら可能と思われるので,$Re=$$1904$の場合の計算をおこなった.

構造は残っているものの,いる.

またVr

それ結果を図 7 に示す.

同じ格子数で

図 7 (a) よりWav$y$な

図 5(a) や図 6(b) に比べてかなり乱れたものになっても図 7(b) に見られるように周期構造よりかなりへだたっている.

このことは計算を長時間続けても同じで,ことを示している.

-7-

別の流れに移る過渡的な状態ではない60

\S 4.

まとめ

8 次精度の上流差分を用いて Tay 1or-Co $uette$ 流で外円筒静止,合の数値シ ミュレーションをおこなった.
Tay 1 $or$ 渦流,

Coue $tte$ 流,

Wav

$y$

$Tay1$

な

その結果,内円筒回転の場レイノルズ数の値に応じて,or 渦流を数値的に得ることができた. 計

算結果を用いたより詳しい解析 (たとえば V$r$の波形のスペク トル分析) は別のSX $-2$

計算流体力学研究所の NEC機会に行う予定である.

なお計算には,を用いた.

10 万ステップ ( $T=$ 100 または 250) で約 3 時間であった.

計算時間は引用文献

1)

Kawamur , T. & Kuwahar , K. 1985:

I nne$r$$a$

rec $tSi$ mu 1 at on $of$ aI. 1923:

G.

wo $ro$ ta $tingcy1in$ der$s$,$St$

ab $i1i$ ty

Ph $i1$.

$Tr$

$of$

a $ns$ .

a $vi$ sc us 1 I $quid$ contaI ne
$0$

Roy.

$S$

urb $u1$ en

Mo er, R. D. , Mo $in$ , P. & Le on ar , A. 1983: A
$d$

$s$

$t$

$s$

pectr alnu mer cal$i$method for the $Navier-Stokes$ Eq. with appl ication to $Taylor-Couette$J.

flow,4)$N$

compt.

phys.

Marcu , P. S.
$s$

52 524-544$1984b$ ;$Si$

mu la $ti$ on $ofT$ ay 1 or-Co $uettef1$ ow, Par $t2$ .

umer cal results for wavy vor ex$t$$i$ne tr av ell ng wave, J.

$f$

low$N$

ume $ri$ ca 1 calc

$w$

ith

$0$

$i$

Fluid Mech. 146, 655)

Har 1 $ow$ , F. H. & We 1 ch, J. E. :$-dIpendent$
Phys.

$t$

be ween$d$oc. A 223 289-343
$-$8)

$T$

$i$

Fl ow by Fi ni te-D I fference Method, A 1AA pap er 85-0376

2) Tay 1 $or$ ,$t$

$Di$

$a$

$vIs$

cous in compr $ssI$ bl$e$$u1$

at on $ofti$ me$i$fl ow of fluid wi th free sur ace,$f$

Fluid 8 (1965) 2182-2189-8-61

Fe ns $ter$ mac $her$ , P. R. , Sw nn $ey$ , H. L. & Go llu , J. P. 1979: Dy na

6)
$i$

$b$

$i$

ns ab $i1iti$ es an
$t$

$d$

$the$

$tr$

an $sition$$to$chao $ticT$ ay 1 $or$$v$

or $tex$

$f1$

ow,

J.

Mech. 94 103-128
7)

Co le, J. Aeffects,

1976:

$T$

ay 1or-vo $rt$ ex$inst$ ab $i1$ I $ty$

J. Fl uid Mech. 75, 1-15.

図1

格子システムと境界条件-9-an$d$

ann$u1$us-l ngt
$e$

$h$

$mi$

ca 1$F1uid$

62$\searrow$$\searrow$$\searrow$$\nearrow$$\nearrow$$\nearrow$$\searrow$$\searrow$$\nearrow$$\nearrow$$\searrow$$\nearrow$

図2

Par $tic$ le

$P$

at

$h$の概念図

VELOC I TY

図3

軸に垂直な画内での速度分布 $(Re= 110)$

-10-

63

図 4 (a)${\rm Re}=238$.

$T=200$でのPar $tic$ le Path

速度の時間変化 $(Re=238)$-11-641

図 5 (a)$T$

ay lo 渦流に撹乱を与えた後 $T=240$$r$$(Re= 238)$-12-での$P$

ar $tic$ le Path

6\geq プ$\triangleright 3$$\vee\infty$$-f3-$

66

図 6(a)

図 6(b)$Re=476$ ,$T=$$100$

でのPar I le Pa $th$
$t$

$c$

$Re=476$ , 撹乱を加えた場合の Part \dager $c1e$ Pat $h(T=200)$-14-67