数理解析研究所講究録
第 1567 巻 2007 年 111-122

111

矩形曲がり管内流における非定常解の 3 次元構造
岡山大学大学院自然科学研究科
渡辺毅 (Takeshi WATANABE), 柳瀬眞一郎 (Shinichiro YANASE)
Graduate School of Natural Science and Technology, Okayama University

1

はじめに
矩形断面曲がり管内流は, 平面クェット流, 平面ボアズイユ流, 円管ボアズイユ流などと並んで,流体運動の基本的な性質を理解するための基礎研究として古くから研究されてきた流れの一つである. 一般に曲がり管の場合, 管路の持つ曲率によって流体に遠心力が働くために, 曲がり管に特有の断面二次流れが形成されることが知られている. これは,流速が管路中央部で大きく管壁付近では小さいことから, 管中心付近で遠心力が最大となり, これによって管中心付近の流体が曲がり管の外壁側に押し出されることによってつくられるものであり, Fig. 2 に示すような複雑な分岐構造を持っている.

こうした分岐解析は Winters $(1987)^{1)}$ によって初めて詳細に研究され, いくつかの典型的な定常解と分岐構造が数値的に求められた. Winters が行ったのは曲率 $\delta=0.02(\delta$ の定義については後述) の場合についてであったが, 後に多くの数値解析が行われ, Winters によって示された分岐構造が曲がり管内流において一般的なものであることが明らかになった. こうした研究のうち最近行われたものをいくつか挙げると, 例えば Wang and Yang$(2005)^{2)}$ が,$\delta=0.01$ に対する実験的,理論的研究を行い, 2 次元数値計算結果と実験結
果が良く一致することを示している. また, Mondal ら $(2006)^{3)}$ は, $\delta=0.1$ の場合に対して徹底的な解析を行い, $\delta=0.25$ 付近で著しい分岐構造の変化が発生することを見出している.

これまでに行われてきた矩形曲がり管内流れの研究はほぼすべて, 2 次元性の仮定のもとに行われたものである. これは数値的には, 流れが管軸方向に一様であるとの仮定をおいて計算することに対応し, また実験的には, 管路のどこか一断面に注目し,ところで,その断面における断面二次流れを可視化することに対応している. この仮定は定常解の場合は妥当的であると考えられるが, 非定常解の場合, 軸方向の一様性を保ったまま振動
する流れは極めて特殊であり, 軸方向にも何らかの構造を持つことが予想される. この点について, Mondal ら $(2007)^{4)}$ は, 2 次元計算で得られた非定常解は 3 次元的には進行波になっていることを予想している. 本研究の直接の動機は, この予想の妥当性を数値的に検証することにあり, 計算を実行した結果, 確かに進行波が生成されることが明らかに
なった.

2

基礎方程式

Fig. 1 に示す座標系 $(x, y, z)$ において, 方向に一様な圧力勾配 によって駆動される流れを考える. 曲がり管の中心から管中心までの距離を R. 管断面の幅および高さの半分
$z$

$G$

112

Fig. 1. 座標系の距離を として, 曲率 を $\delta=d/R$ で定義する. 基礎方程式は, 非圧縮粘性流体を仮定
した連続の方程式と Navier-Stokes 方程式および渦度方程式の無次元形
$d$

である.

$\delta$

ここに $u,$

$t,$

$P,$

$\nabla\cdot u=0$(1)$\frac{\partial u}{\partial t}+(u\cdot\nabla)u=-\nabla p+\nabla^{2}u+Gk$(2)$\frac{\partial w}{\partial t}+\nabla x(w\cross u)=\nabla^{2}w$

(3)

$w,$

次元渦度. 無次元圧力勾配,, ゐはそれぞれ無次元流速, 無次元時間, 無次元圧力, 無はおよび方向の単位ベクトルを表し,

$G$

$z$

$\nabla$$\nabla=\frac{\partial}{\partial x}+\dot{\mathcal{J}}^{\frac{\partial}{\partial y}+k\frac{1\partial}{1+\delta x\partial z’}}$で与えられる. ただし $i,$
上ですべりなし条件$i$$\nabla^{2}$$\nabla^{2}=\frac{\partial^{2}}{\partial x^{2}}+\frac{\delta\partial}{1+\delta x\partial x}+\frac{\partial^{2}}{\partial y^{2}}+\frac{1\partial^{2}}{(1+\delta x)^{2}\partial z^{2}}$はそれぞれ $x,$$y$

方向の単位ベクトルである. 境界条件は壁面$u(\pm 1,y,z,t)=u(x, \pm 1, z,t)=0$を課す. なお各量の無次元化はそれぞれ$u’= \frac{\nu}{d}u$,$t’= \frac{d^{2}}{\nu}t$,$p’= \frac{\rho\nu^{2}}{d^{2}}p$,$(x’,y’,z’)=(dx,dy,dz)$(4)

によって行い, 唯一の駆動力である無次元圧力勾配 $G$ は$G= \frac{d^{3}}{\rho\nu^{2}}\frac{\partial p’}{\partial z’}$によって与える. ただしプライム (’) を付した変数は次元を持つ量を表す. またれぞれ動粘性係数, 密度である. また Dean 数は $Dn=\sqrt{2\delta}G$ で定義される.

$\nu,$$\rho$

はそ

113

3

数値解法

3.1

2 次元計算2 次元性, すなわち$\partial/\partial z=0$を仮定すれば; 曲がり管座標において連続の方程式 (1) は$\frac{\partial u}{\partial x}+\frac{\delta u}{1+\delta x}+\frac{\partial v}{\partial y}=0$と表示されるから,これを用いて流れ関数$u= \frac{1\partial\psi}{1+\delta x\partial y}$

で導入し, これを渦度方程式 (3) のことにより, 解くべき方程式$z$$\psi$,を$v=- \frac{1\partial\psi}{1+\delta x\partial x}$

成分と Navier-Stokes 方程式 (2) の$z$成分に代入する$( \nabla_{2}^{2}-\frac{\delta\partial}{1+\delta x\partial x})\frac{\partial\psi}{\partial t}$$= \nabla_{2}^{4}\psi-\frac{2\delta}{1+\delta x}\nabla_{2}^{2}\frac{\partial\psi}{\partial x}-\frac{3\delta^{3}\partial\psi}{(1+\delta x)^{3}\partial x}+\frac{3\delta^{2}\partial^{2}\psi}{(1+\delta x)^{2}\partial x^{2}}-\frac{1\partial(\nabla_{2}^{2}\psi,\psi)}{1+\delta x\partial(x,y)}$$- \frac{\delta\partial\psi\partial^{2}\psi}{(1+\delta x)^{2}\partial x\partial x\partial y}+2\delta w\frac{\partial w}{\partial y}$$+ \frac{\partial\psi}{\partial y}\{\frac{2\delta}{(1+\delta x)^{2}}\nabla_{2}^{2}\psi+\frac{\delta\partial^{2}\psi}{(1+\overline{\delta}x)^{2}\partial x^{2}}-\frac{3\delta^{2}\partial\psi}{(1+\delta x)^{3}\partial x}\}$$\frac{\partial w}{\partial t}=G+\nabla_{2}^{2}w+\frac{\delta\partial w}{1+\delta x\partial x}-\frac{\delta^{2}w}{(1+\delta x)^{2}}-\frac{1\partial(w,\psi)}{1+\delta x\partial(x,y)}-\frac{\delta w\partial\psi}{(1+\delta x)^{2}\partial y}$が得られる. 解くべき変数はを用いて$\psi,$$u$の 2 個である. これら 2 変数を, 展開関数$\Psi_{t}$および$\Phi_{l}$,$\psi(x, y, t)\approx\sum_{l=0}^{L}\sum_{m=0}^{M}\Psi_{l}(x)\Psi_{m}(y)\psi_{lm}(t)$$w(x, y, t) \approx\sum_{l=0}^{L}\sum_{m=0}^{M}\Phi_{l}(x)\Phi_{m}(y)w_{lm}(t)$のように展開し, これに選点法を適用して解く. ただし $L,$
たおよび
は, Chebyshev 多項式
を用いて$\Psi_{l}$$\Phi_{l}$$\Psi_{l}(x)=(1-x^{2})^{2}T_{l+1}(x),$と定義され,$\Lambda’I$

は打ちきり項数である. ま$\tau_{\iota}$$\Phi_{l}(x)=T_{l+1}(x)-\frac{1}{2}\{1-(-1)^{l}\}x-\frac{1}{2}\{1+(-1)^{l}\}$

自動的に境界条件を満たしている. なお
満たす他の定義として$\Phi_{t}$(6)の定義については同様の条件を$\Phi_{l}(x)=(1-x^{2})T_{l+1}(x)$も考えられるが, 計算を実行した結果, 時間発展計算に必要となる逆行列の精度が, 式(5) で示す定義を用いた場合の方がはるかに良くなることが明らかになったため式 (5) に示す定義を選んだ. ただしこの精度の問題は 2 次元計算ではほとんど差がなく, 3 次元計算の場合に顕著となるものである.
2 次元計算では, Newton-Raphson 法を用いて定常解を, また求めた定常解について線形安定性解析を行い, さらに時間発展計算を実行する. 時間発展計算においては, 非線形
項に Adams-Bashforth 法, 高次微分項に対して Crank-Nicolson 法を適用する.

114
3.2

3 次元計算3 次元計算では, まず流れ場 $u$ と圧力場

$P$$u(x,y,z,t) \approx\sum_{n=-\#}^{2}e^{iknz}u_{n}(x,y,t)u_{-1}$

を

$z$

方向に

$p(x,y,z,t) \approx\alpha_{-1}\sum_{n=-\#}^{2}e^{iknz}p_{n}(x,y,t)$は波数, $N$ は打ち切り項数である.
のようにフーリエ展開する.. ただし は虚数単位,を用いては, さらに 2 次元計算で用いたのと同じ展開多項式
$k$

$i$

$\Phi_{l}$$u_{\mathfrak{n}}$(6)$u_{n}(x,y,t) \approx\sum_{l=1m}^{L}\sum_{=1}^{M}\Phi_{i}(x)\Phi_{m}(y)u_{lmn}(t)$のように展開される. ただしこのように展開した場合, 流れ場の境界条件としてすべりな
し条件は課されるものの, 連続の方程式から導かれる境界条件(7)$\frac{\partial u_{n}}{\partial x}|_{x=\pm 1}=\frac{\partial v_{n}}{\partial y}|_{y=\pm 1}=0$は自動的には課されない. この点については Appendix A.1 で述べる.
次に変数の消去法を考える. 解くべき変数は流れ場の 3 成分 $u_{n}=(u_{n}, v_{n}, w_{n})$ およびとおよびを消去し, 実際にはの 4 個であるが, 計算の効率化のために
圧力場の 2 変数についての方程式を解く.
の消去は連続の方程式 (1) を用いる. 式 (1) は第 モードに対して
変数$v_{n}$$p_{n}$$u_{n}$$p_{n}$$w_{n}$$n$$u_{n}$$\frac{\partial u_{n}}{\partial x}+\frac{\delta u_{n}}{1+\delta x}+\frac{\partial v_{n}}{\partial y}+\frac{iknw_{n}}{1+\delta x}=0$を含む項を左辺, それ以外を右辺に移項してと書かれるから,展開し, さらにこれに選点法を適用することで$u_{n}$$u_{n}$

を展開多項式$\Phi_{l}$で$\sum_{l,m}[\frac{\partial\Phi_{l}}{\partial x}\Phi_{m}+\frac{\delta}{1+\delta x}\Phi_{l}\Phi_{m}]_{1j}u_{lmn}=[-\frac{\partial v_{n}}{\partial y}-\frac{iknw_{n}}{1+\delta x}]_{\dot{t}\dot{j}}$とする. ただし添字毎は選点$(x_{t}, y_{j})$行列を両辺に左から乗じることで$u_{n}$

における値を表す. この式で左辺の展開多項式の逆の展開係数 $u_{lmn}$ が(8)$u_{lmn}= \sum_{1j}\{[\frac{\partial\Phi_{l}}{\partial x}\Phi_{m}+\frac{\delta}{1+\delta x}\Phi_{l}\Phi_{m}]_{ij}\}^{-}1[-\frac{\partial v_{n}}{\partial y}-\frac{iknw_{n}}{1+\delta x}]_{\dot{J}}$によって表され, これによってが消去される.
次に Navier-Stokes 方程式の および 成分を直接積分することによってる. Navier-Stokes 方程式 (2) の および 成分の第 モードはそれぞれのように$v_{n}$と$u_{n}$$w_{n}$

$x$

$x$

$y$

$y$

$p_{n}$$n$$\frac{\partial u_{n}}{\partial t}+[(u\cdot\nabla)u-\frac{\delta w^{2}}{1+\delta x}]_{n}=-\frac{\partial p_{n}}{\partial x}+\nabla^{2}u_{n}-\frac{\delta^{2}u_{n}}{(1+\delta x)^{2}}-\frac{2\delta V^{a}knw_{n}}{(1+\delta x)^{2}}$

を消去す

115

$\frac{\partial v_{n}}{\partial t}+[(u\cdot\nabla)v]_{n}=-\frac{\partial p_{n}}{\partial y}+\nabla^{2}v_{n}$と書かれるから, まず成分について, 圧力項を左辺, それ以外を右辺に移項して

$x$

$x$

で積

分することにより圧力場が$p_{n}= \int_{-1}^{x}\{-\frac{\partial u_{n}}{\partial t}-[(u\cdot\nabla)u-\frac{\delta w^{2}}{1+\delta x’}]_{n}+\nabla^{2}u_{n}-\frac{\delta^{2}u_{n}}{(1+\delta x’)^{2}}-\frac{2\delta iknw_{n}}{(1+\delta x)^{2}}\}dx’+f_{n}(y)$に関する不定性 $f_{n}(y)$ を残して決定される. そこで寄与のある項のみを同様に積分することで
のように

$y$

$y$

成分のうちこの $f_{n}(y)$ に$f_{n}(y)= \int_{-1}^{y}[\frac{\partial^{2}v_{n}}{\partial x^{2}}+\frac{1}{1+\delta x}\frac{\partial v_{n}}{\partial x}]_{x=-1}dy’$が得られ,これによって圧力項が消去される. なお具体的な数値積分の方法についてはについてAppendix で詳述する. 以上のようにして$v.,$$Pn$ が消去されたため,Navier-Stokes 方程式 (2) と渦度方程式 (3) の 6 式のうち, まだ使っていない 4 式から適当な 2 式を選び連立させて解く. 本研究においては Navier-Stokes 方程式の 成分$B$$u_{n}$.

$w_{n}$$z$$\frac{\partial w_{n}}{\partial t}+[(u\cdot\nabla)w+\frac{\delta uw}{1+\delta x}]_{n}=-\frac{iknp_{n}}{1+\delta x}+\nabla^{2}w_{n}-\frac{\delta^{2}w_{n}}{(1+\delta x)^{2}}+\frac{2\delta iknu_{n}}{(1+\delta x)^{2}}$と渦度方程式の

$z$

(9)

成分

$\frac{\partial\omega_{zn}}{\partial t}+[(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+u\cdot\nabla I^{\omega_{z}}]_{n}$$+[ \frac{2\delta w}{1+5x}\frac{\partial w}{\partial y}-\frac{\delta w\partial_{t^{1}}}{(1+\delta x)^{2}\partial z}+\frac{1\partial w\partial v}{1+\delta x\partial x\partial z}-\frac{1\partial w\partial u}{1+\delta x\partial y\partial z}]_{n}$(10)$= \{\nabla^{2}-\frac{\delta^{2}}{(1+\delta x)^{2}}\}\omega_{zn}+\frac{2\delta ikn\omega_{xn}}{(1+\delta x)^{2}}$の 2 式を用いる. ただし$\omega_{x},$$\omega_{z}$は渦度の $x,$$z$$\omega_{x}=\frac{\partial u1}{\partial y}-\frac{1\partial v}{1+\delta x\partial z}$,

成分であり, 曲線座標系において$\omega_{z}=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$と表示される.

この 2 式を用いて圧力勾配 $G=3400$ の場合に対して時間発展計算をおこなう. ただし
非線形項に Adams-Bashforth 法, 高次微分項に対して Crank-Nicolson 法を適用し, アラ
イアス誤差は 3/2 法によって取り除く.
なお本節で述べたような数値計算手法を選択するに至るまでに, 他の数通りの数値計
算方法を試したものの, それらはどれもうまくいかなかった. このことに関する詳細はAppendix A2 で述べる.

116

4

計算条件

本研究においては, $x,$ $y,$ 方向の打ちきり項数 $L,$ $M,$ $N$ をそれぞれ $L=24,$ $M=24$ ,$N=64$ とする. なお $L,$ $M$ は, 2 次元計算によって得られる定常解の収束性によって決定した. 打ちきり項数を上記のように選んだ場合, $G=3400$ における定常解に対する高$z$

次モードからの寄与はのように評価されることから, 十分なモード数であると判断される.
方向の打つ切り項
数 $N$ については, 図 1 によりゐ $n>1$ からの寄与が極めて小さいことが分かるため, 十分$z$

なモード数であると判断した.
また曲率 および波数んは, 同時に行っている可視化実験の条件に合わせるために
0036, ゐ $=0.048$ を選ぶ. ゐの値は曲がり管 3/4 周分に対応するものである.
$\delta$

5

結果

5.1

2 次元計算$\delta=$

本研究の直接の動機は, 2 次元計算によって得られた非定常解が 3 次元的に進行波を形成しているかどうかを確認することであった. このため, まずは同じ条件の曲がり管について 2 次元計算を行い, どの圧力勾配領域で非定常解が得られるかを検討し, この結果を用いて 3 次元計算の計算条件を決定した.
Fig. 2 は, 2 次元計算によって得られた定常解の分岐図である. 図中で, 解が存在する
領域を I,II,III の 3 種類に分割しているが, I は, ただ一つの安定な定常解が存在する領域,
II は, 安定な定常解が存在せず, かつ, 1 個の対称な不安定解と反対称な 1 対の不安定解の, 合計 3 個の解が存在する領域, また III は, 安定な定常解が存在せず, かつ, 3 個以上の不安定な定常解が存在する領域を示している.
2 次元時間発展計算を実行した場合, 領域 I ではただ一つの安定な定常解に収束し, II,
III では非定常解が得られるが, 特に領域 II において, 存在する 3 個の定常解のまわりに,Fig. 3 に示すような, 非常に大きく変動する安定な周期解が構成される. 実際, Fig. 3 における $t=0.158$ と 0.554 の解は, Fig. 2 に示す 4 渦対称解に非常に近く, また $t=0$ と0.396 の解は 2 渦反対称解に非常に近いことがわかる. これは単に定性的に良く似ているというだけでなく, 渡辺ら $(2006)^{5)}$ は解の状態ベクトルを定量的に検討し, 確かにこれら
3 個の不安定な解の周りに構成される周期解になっていることを示している. また曲率は異なるが, Wang and Yang $(2005)^{2)}$ が 2 次元数値計算によって求め, 実験との一致を確認したのも, この解と同じ領域 II に存在し, 定性的に同じく安定な周期解に対してである.
以上の検討から, 3 次元計算を実行する圧力勾配として $G=3400$ を選ぶことにする.

117

圧力勾配

$G$

Fig. 2. 2 次元計算による分岐解析結果. 1st branch および 3rd branch の太線部は安定な定常解.

Fig. 3. $G=3400$ における 2 次元周期解 (周期 $T\approx 0.791$ )5.2

3 次元計算3 次元計算を実行することにより確かに進行波が存在することが確認された. Fig. 4 は,得られた進行波を流脈線によって可視化したものである. 流れは右から左に向かっており, 右手前側が外壁側となっている. 図の中央付近にある, 流脈線が大きく盛り上がっている部分が流れ方向に伝播していく. この盛り上がった部分を除く部分, つまり図中の右端側と左端側ではほぼ対称な流れの状態となっていることから, 進行波の発生によって,ほぼ対称な流れの状態における空間対称性が不安定となってくずれ, 対称性のくずれが進

行波として伝わっていくことが分かる.
進行波の性質を定量的に把握するためにスペクトル解析を実行した結果が Fig. 5 である. スペクトル解析を行う上では, 例えば各波数とこのエネルギーなど, 流れ場全体を代表する量をとるのが望ましいが, 実際に計算してみたところ, 波数ごとの全領域におけるエネルギーはほとんど変化しておらず, 状態ベクトルが複素平面内で単振動しているような状態であることが明らかになった. このため進行波の性質を把握するための量とし$u_{n}$

118

Fig. 4. 進行波の流脈線による可視化

Fig. 5. 進行波のスペクトル解析結果

2 個の代表速度 $v_{n}(0,0.5)$ と $w_{n}(0,0)$ を選び, これらの実部に対してスペクトル解析を行った.
Fig. 5 においてまず注目すべき点は, 図 (A) と (B) がそれぞれ異なる代表速度を用いて得られたものであるにもかかわらず, 両方とも, スペクトルがほぼ同じ傾きを持つ直線上に分布していることである. このことから, 計算を行った条件における曲がり管内流においては, ただ一つの進行波成分が存在していることが分かる. この傾きを最小二乗法によって求めることにより位相速度 $u_{\phi}=501.31$ が求められた. ただし Fig. 5-(A) の $n=0$
モードはこの直線上に乗らないため, 位相速度を求める上では考慮していない.
一方, (A) では $n=0$ モードの振動成分が無視できない大きさを持っているのに対し,(B) ではこれがほとんど見られない. これは, (A) は上下方向流速であるため, 断面 2 次
流れの変化に軸方向一様成分が存在することを示すのに対し, (B) は軸方向流速であるため, 全領域にわたる流量の変化は存在しないことに対応していると考えられる.
て,

119
5.3

実験結果との比較

まず 2 次元計算結果について考える. Fig. 3 に示すように, 圧力勾配 $G=3400$ に対して 2 次元計算で得られる周期解の周期は $T\approx 0.791$ である. 本研究では式 (4) のように無次元化しているため, これを次元のある量 $T’$ に戻すことを考える. 実際に用いる量として, 筆者らが行っている曲がり管内流の可視化実験に用いている実験装置の値として$d=0.0125[m]$ , \mbox{\boldmath $\nu$}=llO-6[m2/s](水の物性値) を用いると$\cross$$T’=T \frac{d^{2}}{\nu}=124[s]$と求められる. ところが現時点で, このような長い周期の振動解は実験によっては観測されていない. 一方で Fig. 5-(B) から進行波の一周期に対応する角周波数は $\omega=73.747$ で与えられることから, 3 次元計算によって得られる周期は $T_{3D}=13.312[s]$ となり, 2 次元計算の場合と約 10 倍も異なっているが, 実験値に近い値と考えられる. 現時点で筆者らが行っている曲がり管内流の可視化実験装置は断面のアスペクト比が 2(縦:横 $=2:1$ ) であり本研究で計算を行った断面形状とは異なるため, 今後実験装置のアスペクト比を 1(正方形) にする力\searrow あるいは数値計算の条件を長方形の場合に変更するかのどちらかを行い,実験によって得られる周期を詳しく調べて比較検討を行う予定である.

A

数値計算法に関する補足

A. 1$u_{n\prime}v_{n}$の展開関数

流れ場砺において,に対してはすべりなし条件のみを課せばよいから, 式 (5) に示
す展開関数
を用いて展開すれば十分であるが,に対してはすべりなし条件だけでなく式 (7) に示す断熱条件をも課さなければならないから, 本来は式 (5) に示す展開関を用いて数$w_{n}$$\Phi_{t}$$u_{n\prime}v_{n}$$\Psi_{l}$$u_{n}= \sum_{l=1}^{L}\sum_{m=1}^{Af}\Psi_{l}\Phi_{m}u_{lmn}$,$v_{n}= \sum_{l=1}^{L}\sum_{m=1}^{M}\Phi_{l}\Psi_{m}v_{lmn}$のように展開する必要があるように思われるが, 本研究では式 (8) を用いてを消去しているため, この消去の手続き中で断熱境界条件式 (7) に対する十分条件が課されているため, 断熱境界条件を展開関数中に入れる必要はない.
$u_{n}$

実際に計算の各ステップごとに数値的収束性を確認するために連続の方程式 (1) に対す
る残差 $R$ を$R=\mapsto^{n---\frac N\sum_{2}^{\not\in-1}(\frac+\frac+\frac+\frac)^{2}}$で定義して確認してみると,$R=O(1\cross 10^{-12})$ で十分に収束していることからように展開するのは妥当であると考えられる.

, 式 (6) の120(a) Navier-Stokes 方程式の

$z$

成分

(b) 渦度方程式の

$z$

成分

Fig. 6. 解の収束性

A.2

解くべき変数と方程式の選択

3.2 節で述べた数値計算法を選択するに至るまでに, 他に数通りの方法を試したものの,いずれもうまく計算を実行することができなかった. そこで本節では, 他に試した方法について述べ, それがなぜうまく計算出来ない力\searrow という点について検討した結果を示す.
曲がり管内流の 3 次元計算を行うにあたって最初に立てたプランは$\bullet$$\bullet$$\bullet$

方向に Fourier 展開すると管軸方向流$n=0$ モードに対しては, 2 次元計算と同じ 2 式を用い, 流れ関数
速 $w0$ を解く$n\neq 0$ モードに対しては,を消去し, 渦度方程式 (3) の連続の方程式を用いてを解く
成分に 成分又は 成分のどちらかを連立させて$z$$\psi_{0}$$z$$w_{n}$$x$というものであった.

$y$$u_{n},$$v_{n}$

しかしこのプランに沿って時間発展計算を実行すると, わずか数ス

テップの計算で発散してしまう. そこで 2 次元計算によって得た定常解を渦度方程式 (3)の各成分に代入し, 解になっているかどうかを確かめたところ, 特に管断面の角の部分で解の誤差が異常に大きくなっていることが明らかになった. この異常なふるまいの様子
を Fig. 6 に示す. 図中 (a), (b) は 2 次元計算で定常解を求めるために用いている 2 式であり, (a) は $10^{-10}$ 程度, (b) は程度のオーダーで非常に良く収束している. ここで,これらが十分に収束しているなら, 渦度方程式 (3) の $x,$成分も収束しているはずだと考えられるが, 図の (c), (d) から明らかなように, 収束していない. 全体としては収束しているが, 管壁部, 特に角の部分で極端に大きくなっている. 時間発展計算が発散するのは, このように管壁部で解が異常値を取ることが原因であると考えられる. そこでこの現$10^{-9}$

$y$

121

象を検証するために, 簡単な例として関数$f(x,y)= \frac{1}{25+x^{2}+y^{2}}$を考え, これを正方形領域で Chebyshev 多項式補間してその高次微分を調べたところ, 同様に, 壁および角の部分で極端に大きな値を取る現象が再現された. またこの現象は, モー

ド数を増やしても改善されないことも同時に明らかになった. このことは第 1 に, 微分階数は可能な限り低い方が良いことを, また第 2 に, 解くべき方程式を選択するにあたって,$n=0$ モードとそれ以外で異なる方程式を用いるのは困難であることを意味している.
ところが, 解くべき変数を 2 個として解くべき方程式を 2 個選ぶ場合, Navier-Stokes 方
程式の 成分に駆動力となる圧力勾配項が含まれていることから, これはどうしても選ばなければならず, そうすると, $n\neq 0$ のとき, 変動圧力の項を解かざるを得ない.
以上の検討から, 32 に示したような数値解法を用いて解く方法を選んだ. ただ, この問題の解決方法として, 例えば選点法ではな \langle Galerkin 法を用いる, あるいは展開関数
として Chebyshev 多項式ではない関数を用いるなど, 他の方法を試行し, より良い方法を考える必要がある.
$z$$\partial p_{n}/\partial z$$B$

圧力項の数値積分方法

精度良く数値積分を行おうとした場合, 通常は Gauss-Legendre 積分が用いられる. ところが展開関数に Chebyshev 多項式を用いている場合, Chebyshev 多項式の補間点と Legendre 多項式の補間点は異なるため, 数値操作が煩雑になる. もし, 区間 [-1, 1] での定積分を求めるなら, Chebyshev 多項式の展開係数から, ただ一度だけ Legendre 多項式の補間点における値を求めれば十分であるが, 不定積分を行おうとした場合,

1. Chebyshev 多項式の展開係数から Legendre 多項式の補間点における値を求める
2. Legendre 多項式の重みを $w_{i}(i=1,2, \cdots, I)$ として, 任意の関数 $f(x)$ の不定積分$F(x)$ の補間点における値を$F(x_{j}) \equiv\int_{-1}^{x_{j}}f(x’)dx’=\sum_{i=1}^{j}f(x_{i})w_{i}$,$x_{j}= \sum_{i=1}^{j}w_{i}-1$,$j=1,2,$ $\cdots,I$のように求める

3. この新しい補間点$x_{j}$

は Chebyshev 多項式の補間点でも Legendre 多項式の補間点でもないため, こんどこの補間点における値 $F(x_{j})$ を用いて Chebyshev 多項式の展開係数を求める4. 得られた展開係数を用いて Chebyshev 多項式の補間点における値を求めるという, かなり煩雑な操作を行う必要がある. 特に上述の手続き中 3 において, Chebyshev多項式の補間点でない点の値を用いて展開係数を求めていることから, この操作において精度が悪くなっている可能性が考えられる.
研究会で報告した時点では上述のような積分方法を用いていたが, その後改良を加えた. その結果, 今まで $n=0$ モードに残っていた振動成分がはるかに小さくなり, 無視し122て良い程度の大きさになったため, 今後, 実験データを得る必要性も含め,さらに詳細に検討する必要がある.
新しい積分方法は次のようなものである.
の
まずは区間 $[-1, x]$ におけるによる積分を考える.
の定積分であるから $p_{n}(x=-1)=0$ を満たす. さらにの境界条件は $n\neq 0$ のとき, 式 (9) から, 壁
面上で$\partial p_{n}/\partial x$$x$$\partial p_{n}/\partial x$$p_{n}$$p_{n}$$p_{n}= \frac{1+\delta x}{i\text{ゐ}n}(\frac{\partial^{2}w_{n}}{\partial x^{2}}+\frac{\delta\partial w_{n}}{1+\delta x\partial x}+\frac{\partial w_{n}^{2}}{\partial y^{2}})$であるから,$p_{n}(\pm 1, \pm 1, t)$=0(複号任意) を満たす. したがって Chebyshev 多項式男を用いて$p_{n}(x, y, t)= \sum_{l=1}^{L}(1+x)(\frac{2-x-y^{2}}{3})T_{l+1}(x)p_{ln}(y, t)$のように展開できるはずである. そこで上式を項別に微分して選点法を適用すれば$\frac{\partial p_{n}}{\partial x}|_{x=x_{j}}=\sum_{l=1}^{L}[\frac{1-2x_{j}-y^{2}}{3}T_{l+1}(x_{j})+(1+x_{j})(\frac{2-x_{j}-y^{2}}{3})T_{l+1}’(x_{j})]p_{ln}$と書ける. ただしプライム (’) は による微分を表す. この式で, 展開関数の逆行列を左が求められ, 最終的に求めたいが,を用いて表される.
から乗じることでのによる積分については, $x=-1$ 上で積分を行うため, $y,$だけの関同様に
$y=\pm1$
数であり,ですべりなし条件と同様の条件が課せるから$x$$p_{ln}$$\partial p_{n}/\partial y$$p_{n}$$\partial p_{n}/\partial x$

$t$

$y$

$p_{n}(y,t)= \sum_{m=1}^{M}\Phi_{m}(y)p_{mn}(t)$を求める. これによりと展開できるはずである. これを用いて の場合と同様にして補間点の変更を行う必要がなく, かつ, 角の部分での境界条件が保証されるため高精度で積分を行うことができる.
$x$$p_{mn}$参考文献
1) K. H. Winters: A bifurcation study of laminar flow in a curved tube of rectangularcross-section, J. Fluid Mech. (1987) vol. 180, pp. 343-369
2) Wang, L. and Yang, T., 2005, Periodic oscillation in curved duct flows, Physica $D$

200, 296-302.
3) R.N.Mondal: Isothermal and Non-isothermal Flows thouugh Curved Ducts with
Square and Rectangular Cross Sections, Ph. D. Thesis, Okayama university (2006)
4) Mondal, R. N., Kaga, Y., Hyakutake, T. and Yanase, S., 2007, Bifurcation diagram
for two-dimensional steady flow and unsteady solutions in a curved square duct, toaPPear in Fluid $Dyn$ . Res..
5) 渡辺毅 t R. N. Mondal, 加賀義人, 柳瀬眞一郎, 2006, 矩形曲がり管内流における振動解の構造, 『日本流体力学会年会 2006 講演要旨集\Delta , AM06-06-010.