数理解析研究所講究録
第 1800 巻 2012 年 184-194

184
急拡大部をもつ管路内の流れの安定性と遷移
同志社大学エネルギー変換研究センター
同志社大学理工学部

村上 真也 (Shin-ya Murakami)
水島 二郎(Jiro Mizushima)

Energy Conversion Research Center, Doshisha University
Faculty of Science and Engineering, Doshisha University

1

はじめに

これまで、流れの安定性理論は流れが定常流から振動流へと遷移する条件を明らかにし、流体運動を支配する方程式の解には安定な解と不安定な解があることを示し、層流から乱流への遷移の機構を
徐々に明らかにしてきた。 しかし、物体後流や境界層・円管ポワズイユ流などについて、線形安定性を支配するオア・ゾンマーフェルト方程式から導かれる臨界レイノルズ数は、実験から得られる定量的な値と著しく異なる値を予測していた。 これは、我々が基本的な流れの安定特性さえも理解できていないことを意味する。 たとえば、 円柱を過ぎる流れが振動流となる臨界レイノルズ数については、

オア・ゾンマーフェルト方程式は実験値の 1/10 程度の小さな値を予測する [1]。
このような状況の下、 プラズマ物理において、主流が存在するときに伝播する波の安定性について、対流不安定撹乱と絶対不安定撹乱という概念が生まれた [2, 3, 4, 5, 6]。ある主流を波束型の撹乱が伝播するとき、波束型撹乱の最大振幅位置が伝播する群速度が正で、 かつ波束型撹乱の中心波数のモードの線形増幅率が正であればその撹乱を対流不安定撹乱と呼び、 主流は対流不安定 (Convective
instability, Convectively unstable) であるという。 これは、 ある有限領域に不安定な撹乱が発生しても、 時間が経てば主流によって撹乱がその領域から流れ去ってしまい、やがてその有限領域は元の状態に戻ることを意味している。 一方、 波束型撹乱の群速度がゼロまたは負で、 かっ中心波数のモードの線形増幅率が正のとき、 その撹乱を絶対不安定撹乱であると呼び、 その主流は絶対不安定(Absolute instability, Absolutely unstable) であるという。 もし流れが絶対不安定であれば、 不安

定撹乱は有限領域に留まり、流れ場に不安定性を引き起こし、流れのレジームの遷移を引き起こす。
これに対して、 Takemoto and Mizushima(2010)[7] は、 これら二つの不安定性の概念を一般化し、パッシブ不安定性 (passive mode in-stability) とアクティブ不安定性 (active mode instability) を定義した。 その結果、 円柱を過ぎる流れにおいては、 全体不安定は波束の振幅の最大位置の伝播速度がゼロになるのではなく、 波束の後端 (tail-) の伝播速度がゼロになるため起きるという描像を得た。 これを詳しく説明する。波束の最大振幅位置は有限の速度で伝播し、 波束の波長が時間的に拡大している状況において、 ある時刻において波束の最大振幅位置よりも後ろの位置う。 点$x_{0}$

における振動振幅$a$$x_{0}$

で波の伝播を観測しよの時間成長率の符号は、波束が正の有限速度で伝播することによる負の寄与と、 波束の波長の拡大による正の寄与によって決まる。 波束の最大振幅位置は正の速度で伝播する場合でも、 波束の波長の単位時間あたりの拡大率が十分に大きい場合は、 波束の後端が負の速度185で移動することになり、 固定点$x_{0}$

における時間振動の成長率は正をとることを意味する。 円柱を過ぎる流れはこのような状況に相当している。
ここでは、 工学や工業においてよく現れる周期的な空間構造を有する管路を簡単化したモデルとして、 急拡大部をもつ管路を流れる 2 次元非圧縮流体の安定性について考える。 急拡大部と急縮小部を周期的にもっ管路では、 急拡大部と急縮小部で構成される管路を 1 ユニットとして、 ユニット数が増えるにしたがって、 より低い臨界レイノルズ数で全体不安定が起こることが知られている [8]。主流に乗って伝播する波束の振幅の線形増幅率 (パッシブ不安定性) とレイノルズ数の関係と、 全体不安定の関係を神田 (2010)[8] は検討し、 急拡大部をもつ管路を流れる流れについても円柱を過ぎる流れの場合と同様の結果を得た。 しかしながら、 数値計算法が主に 2 次精度に基づいていること、 時間や空間刻みを変えたときの数値解の収束性が $+$ 分に検討されていなかったことが問題であった。 この報告では、 周期的に急拡大部をもつ管路流れのパッシブ不安定性を調べた神田 (2010) の数値実験の数値的収束性を検討するために、 新たに開発した空間について 4 次精度の差分に基づいたコードによって、 先行研究 [8] の数値実験の一部を追試した結果を述べる。 追試したのは、 周期的な急拡大部
を $N=1\sim 5$ 個もつ管路の、 周期境界条件および、 系の左端で流入条件、 系の右端で流出条件を課した開放系の数値実験で、 安定性の遷移に関わる臨界レイノルズ数を調べた。

2 基礎方程式と計算方法
2.1

系の形状

急拡大部をもつ管路として、 図 2 に示した管路を用い、 この管路を流体が左から右に流れる様子を調べる。 座標系として、 管路の中心を $y=0$、

管路の左端を $x=0$ とする直交直線座標系を用いる。

急拡大部をもつ管路は、 図 1 に示すように、 狭窄部と急拡大部の組合せを 1 ユニットとして $N$ ユニット(図 2 の場合は $N=5$ ) がつながっている。

狭窄部の

$x$ 、

$y$

方向の幅をそれぞれ路の形状は (1) 拡大比 $E=D/d$$l$、$2d$ 、

急拡大部の

$x$ 、

$y$

方向の幅を

$L$
、

$2D$

とすると、 管(2) 急拡大部のアスペクト比 $A=L/(2D)$ (3) 狭窄部の無次
元長さ $s=l/d$ の 3 つの無次元パラメターで表せる。 このとき、 1 ユニットの方向の長さは、、$x$$\ovalbox{\tt\small REJECT}=2(s+EA)$である。 本研究では、 これまでによく研究されている $E=3$、$A=7/3$、$s=0.5$の場合 (このとき、 $\sim=15$ ) を調べた [8, 9, 10, 11, 12]。

代表長さを狭窄部の管の$y$方向の幅の 1/2 、つまり$d$として、代表速度を流入部の最大速度 $U$ とし、 レイノルズ数 $Re$ を $Re:=Ud/\nu$ と定義した。 ここで$\nu$は粘性係数である。

また、 流れを周期的に配置された観測点で点 $P_{i}(i=0,1, \ldots, 4)$ で観察するとき、 これらの点の$x$座標を周期境界条件の場合は、 (13, 28, 43, 58, 73)、流入・流出条件の (18.5, 33.5, 48.5, 63.5, 78.5) と
した。

186
22 基礎方程式
基礎方程式として、 2 次元非圧縮性流体の渦度方程式と渦度 $\omega(x, y, t)$ と流れ関数 $\psi(x, y, t)$ に関するボアソン方程式$\frac{\partial\omega}{\partial t}+J(\omega,\psi)=\frac{1}{Re}\nabla^{2}\omega$

(1)

,

(2)

$\nabla^{2}\psi=-\omega$を用いた。 ここで、 $J(f,g):=\partial(f, g)/\partial(x, y)=\partial_{x}f\partial_{y}g-\partial_{y}f\partial_{x}g$ はヤコビアン、 $Re$ はレイノルズ数である。 この場合、流れ関数と速度ベクトル $u=ui+vj$ との間には$u= \frac{\partial\psi}{\partial y}$の関係がある。 ここで、

$i,$

$j$

渦度と流れ関数の定常解はそれぞれ $x,$$\overline{\omega}(x, y),$$y$$\overline{\psi}(x, y)$,(3)$v=- \frac{\partial\psi}{\partial x}$方向の単位ベクトルである。

を求めるときは、 (1) 式から時間微分項をゼロと置いた

方程式系$J( \overline{\omega},\overline{\psi})=\frac{1}{Re}\nabla^{2}\overline{\omega}$

(4)

,

(5)

$\nabla^{2}\overline{\psi}=-\overline{\omega}$

を用いた。

図1

$=$

急拡大部と急縮小部をもつ管路の 1 ユニット分を表した図。

識．想

$=$

・韓 r3”-, 紫 v i y$=$$\hat$$=\grave$d 煽繍・．m‘ず 醗$\hat$$\backslash \nwarrow$.
$\cdot\ovalbox{\tt\small REJECT}_{\backslash tY}*$
でち

サお$P_{1}$$P_{2}$$P_{4}$$P_{3}$$\Psi,\ldots$.-. な$-:^{\backslash }\dot{A}$$\backslash \backslash$$P_{5}$

方向の計算領域は急拡大部を有する管路部分 (流入口から最後のユニットの急縮小部の左端) の 2 倍程度の長さをとっているが、ここでは省略して描いて

図2

いる。

$N=5$ のときの管路全体の形状。

$x$

187
23 境界条件
境界条件としては、 壁ですべりなし条件を適用する。 管路左端 $x=x_{0}=0$ の流入口では十分に発達したボアズイユ流 $u(x_{0}, y)=1-y^{2}$ を仮定し、 管路右端の流出口 $x=x_{1}$ ではゾンマーフェルトの放射条件$\frac{\partial\psi}{\partial t}+c\frac{\partial\psi}{\partial x}=0$を用いる。 ここで、$c$として流出口における

$x$

,

(6)

$\frac{\partial\omega}{\partial t}+c\frac{\partial\omega}{\partial x}=0$

方向速度 $u(x_{1}, y)$ を用いた。

24 数値計算法
数値計算は、 空間については等間隔格子点上の有限差分法を用いた。 空間の離散化については粘性

項には 4 次精度の中央差分を用い、 移流項は 4 次精度の荒川ヤコビアンを用いた [13, 14]。 4 次精度の荒川ヤコビアンは、 エネルギーとエンストロフィーの両方を 4 次精度で保存する移流項の差分近似である。 4 次精度の荒川ヤコビァンは壁の内部の渦度を必要とするため、 壁付近で特別な取り扱いが必要である。 ここでは、 壁ですべりなし条件が満たされるように、 壁に隣接した格子点上の渦度によって、 壁の上で壁に沿った方向に誘導される流れがゼロとなるように、 壁のすぐ内側の点での渦度

は壁のすぐ外の格子点の渦度と同じ渦度をもつとした。
壁の上で渦度を求めるためには、 壁の上で流れ関数がゼロであることを用いて構成した 4 次精度の2 階微分の差分近似を用いた。 壁の角での渦度をどのように取り扱うかは重要であるが (4 節参照)、
ここでは簡単に、 壁の角の周囲の格子点の 4 点の平均値を用いた。

定常問題を解く際には、 定常渦度方程式とボアソン方程式に対して SOR 法を用いた。 定常解には$y$方向に対称性をもつ場合ともたない場合がありうる。

関数値が

$0$

$y$

方向に対称性をもつ場合は、 $y=0$ で流れなので、 反復計算中に $y=0$ での流れ関数値を$0$とし、 対称条件が満たされるように解いた。解の収束条件には、 ボアソン方程式と定常渦度方程式に反復解子点で最も大きい絶対誤差が$10^{-10}$$\psi$

を代入したときに全ての格以下となることを採用した。

時間発展問題を解く際には、 時間積分に 4 次ルンゲークッタ法を、 ボアソン方程式は前処理付き

BiCGStab 法を用いて解いた ([15, 16])。前処理は、 対角要素によるスケーリングと、 並列化して計算を行うために局所 ILU(0) 前処理を用いた ([15, 16])。並列化実装には OpenMP を用い、領域を単純に分割して各 CPU コアに割り当てた。 収束条件は、 残差ベクトルの$L_{2}$

ノルムが $10^{-30}$ 以下となることとした。 このとき、解を離散化したボアソン方程式に代入したときの残差の絶対値の最大値
は、 おおよその計算時刻において $10^{-12}$ 程度であった。

空間刻みとしては $\triangle x=\triangle y=0.1$ を用い、 時間刻みとしては主に $\Delta t=0.01$ を、 一部の高レイ

ノルズ数の計算については $\Delta t=0.0025$ を用いた。

188
3 結果
ユニット数が $N=1\sim 5$ の場合について、 周期境界条件の場合と流入流出条件の場合の流れの安定性をそれぞれ調べた。特に、 流入・流出条件の場合の流れの安定性について述べる。

3.1 流入流出境界条件の場合の流れの安定性
図 3 は、 ユニット数が $N=1\sim 5$ の場合に、点 $P_{i}(i=1,2, \ldots, 5)$ でのズ数依存性をプロットしたものである。

$v$

が

$0$

$y$

方向速度

$v$

のレイノル

でないということは、すなわち、安定な定常蛇行流が存在するレイノルズ数を示す。
ユニット数がどの場合でも、 レイノルズ数がおよそ 40 未満のときには対称定常流が安定となるが、 およそ $Re\approx 45\sim 50$ でピッチフォーク分岐が起こり、安定な蛇行流が発生する。 また、 およそ$Re\approx 63\sim 68$

で逆ピッチフォーク分岐が起こり、安定な蛇行流が存在しないレイノルズ数領域が再び現れる。 また、 各々の$v$のパターンは、総ユニット数によらずに、何番目のユニットであるかによって決まっているように見えることに注意したい。 これは、各々のユニットの形状は周期的に繰り返されて同じであっても、 流れ場は周期的でないこと、 そして、 あるユニット目の流れ場は、 より手前のユニットの流れの影響によって規定されていることを示唆する。 $N=2,3,4,5$ の 2 ユニット目
以降は、

$v$

のレイノルズ数依存性が不連続に変わっていることから、 この逆ピッチフォーク分岐は亜臨界であることを示唆している。
図 4 は、 総ユニット数 $N$ とレイノルズ数のパラメタ空間内の時間的安定性を示す。 計算の結果、安定なら $O$、

不安定なら$\cross$、

数値誤差の範囲内でほぼ同じ振幅で振動しているときや判別に必要な計算時間が不十分など判別できなかった場合を
$N=1$

、

$N=2$ と、 $N=4$、$\triangle$で表した。 N-Re パラメタ空間内での中立曲線は、$N=5$ で不連続な印象を受ける。 逆に、 $N=2,3,4$ はあまり変化していないような印象を受ける。

4 議論$N=1$ で流入流出境界条件を用いた計算のピッチフォーク分岐の臨界レイノルズ数
ピッチフォーク分岐の臨界レイノルズ数$Re_{p1\text{、}}$

ホップ分岐の臨界レイノルズ数$Re_{h}$$Re_{pO\text{、}}$

逆の先行研究と本報告を表に 5 表す。 なお、 本報告の計算ではホップ分岐の臨界レイノルズ数を特定するには至らなかったが、 $Re=1000$ で安定、 $Re=1100$ で不安定であった。

これらのことから、少なくとも 3 つの独立したコードの計算結果を見る限り、低レイノルズ数の計算では結果は定量的にほとんど変わらず、 数値計算法の細部に依存しないことが予想される。

189$\mathscr{K}$$*_{\aleph_{\infty}}$$\theta^{\aleph^{\aleph^{\phi}}}\infty\propto$$1^{\theta^{\rho\infty\propto}}$$*_{\infty}$$\aleph*_{\infty}$$\varphi^{\phi}\propto$$\ovalbox{\tt\small REJECT}^{*\infty}$$i^{t^{\theta^{\mu\infty*}}}$$\infty$,$*_{*_{\infty}}$$\aleph**\infty$

455055606570$Re$ , 3rd unlt

455055606570$\aleph_{\aleph\infty\propto}$

45 50 55 60 65 70$Re$ , lSt unit

図3

$”$

5055606570$Re$ , 2nd unit$Re$

, 4th unit

455055606570$Re$ , 5th unit

流入流出境界条件の下での対称定常流の安定性。 上から $N=1,2,3,4,5$ の場合。 左端から 1 ユニット目、 2 ユニット目と順に並べた。

また、 高レイノルズ数では今回の報告で行った計算自体が壁付近の境界層を十分解像しているかどうかは疑問であるため、 結果の正しさは主張できないが、 少なくとも、 数値計算法境界条件の取り

扱い方解像度に流れの安定性が大きく依存する結果であるといえる。
$N=5$ で流入流出境界条件を用いた計算では、神田 2010[8] では、 全体不安定の臨界レイノルズ数はおよそ 2175 と求められていたが、 本報告の計算では、 およそ 200 から 230 の間にあることまでは分かり、 互いに矛盾しない結果となった。

壁付近で薄い境界層が形成されるため、 壁付近での格子間隔 (解像度) が境界層の構造を自然に表現できるかどうかを決める。 この境界層は壁から離れた領域の運動にある程度の影響を与えるだろう。 現在の格子幅では、 この境界層の構造を十分に表現できているとはいい難い。 今回は格子幅を変

190
1600
$\cross$

1400$\cross$

1200$0\S\cross$

1000

2

800

$0$

$\Vert$$\circ 8$600$\hat{O}\cross$

$0$

400

$0$

$o$

200

$\ovalbox{\tt\small REJECT}$

$0$

0123

4

5

6

$N$

図4

総ユニット数 $N$ と $Re$ 数のパラメタ空間内の流入流出境界条件の下での流れの時間的な安定性。 安定な場合を $O$ で、 判別できなかった場合をMizushima,$et$神田 (2010)[8]
本報告
図5

al. $(1996)[9]$$\triangle$

で、 不安定な場合を$\cross$で表した。

$Re_{p0}$$Re_{p1}$$Re_{h}$

47.70
493
489

65.24
657
681

843
8454
$\leq 1100$$N=1$ で流入・流出境界条件を課した場合の流れの安定性に関わる臨界レイノルズ数。

えた計算を行っていないので、 結果の数値的収束性を検討していない。
また、 壁の角で渦度が不連続になる、 すなわち、 角が特異点になっているという問題もある。壁の角は二枚の別々の壁が合わさる端であるが、 各々の壁に沿って渦度を外挿して求めると、 壁での渦度は一般に異なる結果となる。 そこで、 角の渦度を単純に取り扱う方法は、次のように色々な方法が考えられる:(1) 角を構成する二つの壁面のそれぞれから片側差分で渦度を推定してそれらの平均をとる、 (2) 1 と同様だが、 壁に沿った線上にある流体の内部領域の格子点を用いて片側差分をとる、 (3)1 と同様だが、 片側差分のかわりに中央差分を用いる。 今回の計算では、 (3) で 2 次中央差分を用いた。 すなわち、 角の渦度を、隣接した周囲の 4 つの格子点の平均値とした。 もっと工夫された方法として、 この特異点に起因する誤差を小さくとどめるには、 特異点周辺の格子幅を小さくすることである [17]。この二つの問題は、壁付近で密な格子幅をもつ不等間隔格子や、有限要素法を利用し、壁付

近で密な面積要素をもつ要素を生成することで実現できるが、 等間隔の差分法を用いたまま壁付近の解像度を高める方法もある [18,19]。この方法は、再帰的に格子を分割して 1 レベルずつ内側の渦度/流れ関数の値を境界値問題として解く数値的方法と、その再帰的な方法に Moffat 渦の解析解 [20]を援用する方法を提案した [19]。彼らは再帰の回数を解析解を援用することで減らすことができることを示した。 ここで注意したいのは、 全体の格子幅を一定 $(\Delta x=\Delta y=0.05)$ にした条件で行った

191
数値解についてのみ彼らは議論しており、全体の格子幅を小さくしたときの数値解や、 あるいは壁付
近の格子幅だけを密にするような方法による数値解との比較が行われていないことである。 彼らはこの方法を急拡大する管路の流れに適用している [18]。
今回は比較的荒い解像度 格子点幅 $\triangle x=\triangle y=0.1)$ を用いて計算を行ったが、 数値的収束性を吟$($味する必要がある。 今回の計算は、 必要に応じて OpenMP にて最大 4 並列で行ったが、 特に時間発展問題について数値的収束性を吟味するには、 効率的な数値計算法を用いるか、 高速な計算機を用いる必要がある。

5

まとめ
急拡大部と急縮小部をもつ流れについて、 アスペクト比 A$=$

7/3、拡大比 E$=$

3、急縮小部の無次

元長さ $s=0.5$ の場合について、 ユニット数を変化させて主に神田 (2010) の一部の数値計算の追試を、 神田 (2010) とは独立に作成したコードを用いて行った。

ピッチフォーク分岐の臨界レイノルズ数は定量的にはほとんど変わらず、 数値計算法や境界条件の取り扱いの細部にはあまり依存しないことが示唆された。 それに対して、 ホップ分岐の臨界レイノル
ズ数は比較的大きく変化した。 したがって、 高レイノルズ数 (およそ 700 以上) の場合には、 異なる

数値計算法や境界条件の取り扱い方によって数値収束性を十分に検討する必要があるだろう。
前節で述べたように、 数値的収束性を検討する必要性があるが、 現在の方法で数値的な収束性を十分に検討するには、 現状の環境ではやや力不足であり、 単純に使用する計算機の性能を上げるか、 使

用するアルゴリズムをより並列度の高い、 あるいは計算量の少ないアルゴリズムに変更するか、 コードの実装を工夫するかのいずれかが必要である。
効率的な数値計算法としては、 単純に格子点幅を一様に細かくする方法に加えて、 上に述べたような角付近で密な格子を用いる方法、 また、 壁付近でのみ密な格子を用いるような非一様格子を採用する方法、 ボアソン方程式の求解アルゴリズムとして別の方法 (例えば代数的マルチグリッド [16, 17])を採用することなどが考えられる。
高速な計算機を用いる具体的な方法としては、 MPI を使って計算機間並列を併用するか、 GPGPUなどのメニーコアアーキテクチャを用いることなどが考えられる。
今後、 GPGPU を用いて、 さらに [18, 19] のように角付近で密な格子を用いたコードを用いて、数値的収束性を検討する予定である。

急拡大部と急縮小部をもつ管路の流れは、他のさまざまなな研究と比較されるべきである。 たとえば、 急拡大部を一つだけもつ流れはよく研究されている [18, 21, 22, 23, 24, 25, 26]。急縮小部を一つだけもつ流れもよく研究されている [27]。また、 急拡大部をも持つ流れは、 二つのバックステップ(backward-facing step) をもつ流れ [23, 28, 29] と解釈することもできる。 これらの間を関連付けて

192
理解することは、 理論的に興味ある話題である。

謝辞
本研究は、 科学技術振興機構 (JST). 戦略的創造研究推進事業 (CREST) の研究領域 : 「数学と諸分野の協働によるブレークスルーの探索」の「現代数学解析による流体工学の未解決問題への挑戦」の援助を受けている。
図の作成には Asymptote http: //asymptote sourceforge net $/)$ gnuplot
$($

、

$($

httP: //gnuplot info$/)$ 、

地球流体電脳倶楽部の GPhys ライブラリ (http: $//ruby$.gfd-dennou.org/products/gphys/) を用いた。数値実験には、 IO ライブラリに地球流体電脳倶楽部の gtool (http://www.gfd-

dennou org/library/gtool/) 及び Unidata の NetCDF (http: $//\mathfrak{m}vw.unidata$ ucar edu/software/netcdf/)を用いた。 これらの開発者及びコミュニティの方々に感謝する。

参考文献[1] Taneda, S., 1963: The stability of two-dimensional laminar wakes at low Reynolds numbers,
J. Phys. Soc. Japan., 18, 288–296.
[2] Briggs, J. R., 1964: Electron-stream interaction with plasmas MIT Press, Cambridge, Chap.
2.
[3] Betchov, R. and Criminale, W. O., 1966: Spatial instability of the inviscid jet and wake,
Phys. Fluids, 9, 359–362.

[4] Gaster, M., 1968: Growth of disturbances in both space and time, Phys. Fluids, 11, 723727.
[5] Huerre, P., Monkewitz, P. A., 1990: Local and global instabilities in spatially developing
flows, Annu. Rev. Fluid Mech., 22, 473–537.
[6] 水島二郎，藤村薫，2003: 流れの安定性，朝倉書店，東京．
[7] Takemoto, Y., Mizushima, J., 2010: Mechanism of sustained oscillations in a fluid flowing
past a circular cylinder obstacle, Phys. Rev. $E,$ $82$ , 056316.
[8] 神田徹平，2010: 周期的急拡大管路流れの不安定性と遷移，修士論文，同志社大学大学院．
[9] Mizushima, J., Okamoto, H., Yamaguchi, H., 1996: Stability of flow in a chamel with asuddenly expanded part, Phys. Fluids, 8, 2933–2942.

[10] Mizushima, J., Shiotani, Y., 2001: Ttansitions and instabilities of flow in a symmetricchannel with a suddenly expanded and contracted part, J. Fluid Mech., 434, 355–369.
[11] Mizushima, J., Yoshida, S., Takaoka, M., 2006: Deflection-Flipping Waves in a Symmetric

Chamel with Spatially Periodic Expansions and Contractions, J. Phys. Soc. Japan., 75,
113401.

193[12] Takaoka, M., Sano, T., Yamamoto, H., Mizushima, J., 2009: Convective instability of flow
in a symmetric channel with spatially Periodic structures, Phys. Fluids, 21, 024105.
[13] Arakawa, A., 1966: ComPutational Design for Long-Term Numerical Integration of theEquations of Fluid Motion:

Two-Dimensional IncomPressible Flow. Part 1, J. Com-

put. Phys., 1, 119–143.

[14] Krishnamurti, T. N., Bedi, H. S., Hardiker, V. M., 1997: An Introduction to Global Spectral
Modeling, Oxford Univ. Press, New York.
[15] Van der Vorst, H. A., 2003: Iterative Krylov Methods for Large Linear Systems, Cambridge Univ. Press, Cambridge.
[16] Saad, Y., 2003: Iterative Methods for SParse Linear Systems, 2nd edition, SIAM, Philadelpha, PA.

[17]J. H. ファーツィガー，M. ペリッチ著，小林敏雄，大島伸行，坪倉誠訳，2003:

コンピュータによる流体力学，シュプリンガー．ジャパン，pp. 176 –178.
[18] Hawa, T., Rusak, Z., 2001: The dynamics of a laminar flow in a symmetric channel with asudden expansion, J. Fluid Mech., 436, 283–320.
[19] Hawa, T., Rusak, Z., 2002: Numerical-Asymptotic Expansion Matching for Computing a

Viscous Flow Around a SharP Expansion Corner Theoret. Comput. Fluid Dynamics, 15,
265–281.
[20] Moffat, H. K., 1969: Viscous and resistive eddies near a sharp corner, J. Fluid Mech., 36,
117–129.
[21] Shapira, M., Degani, D., Weihs, D., 1990: Stability and existence of multiple solutions forviscous flow in suddenly enlarged channels, ComPut. Fluids, 18, 239–259.
[22] Feam, R. M., T. Mullin, T., Cliffe, K. A., 1990: Nonlinear flow Phenomena in a symmetricsudden expansion, J. Fluid Mech., 211, 595–608.
[23] Gagnon, Y., Giovannini, A., H\’ebrard, P., 1993: Numerical simulation and physical analysis
of high Reynolds number recirculating flows behind sudden expansions, Phys. Fluids $A,$ ,2377–2389.
[24] Drikakis, D., 1997: Bifurcation Phenomena in incompressible sudden expansion flows,$5$

Phys. Fluids, 9, 76-87.

[25] Alleborn, N., Nandakumar, K., Raszillier, H., Durst, F., 1997: Further contributions on thetwo-dimensional flow in a sudden expansion, J. Fluid Mech., 330, 169–188.
[26] Rusak, Z., Hawa, T., 1999: A weakly nonlinear analysis of the dynamics of a viscous flow
in a symmetric channel with a sudden exPansion, Phys. Fluids, 11, 3629–3636.
[27] Chiang, T. P., Sheu, T. W. H., 2002: Bifurcationms of Flow Through Plane Symmetric

Chamel Contraction, J. Fluid Eng., 124, 444–451.
[28] Kaiktsis, L., Karniadakis, G. E., Orszag, S. A., 1996: Unsteadiness and convective insta-

bilities in two-dimensional flow over a backward-facing step, J. Fluid Mech., 321, 157187.

194[29] Sch\"afer, F., Breuer, M., Durst, F., 2009: The dynamics of the transitional flow over a
backward-facing step, J. Fluid Mech., 623, 85–119.