数理解析研究所講究録 1231 巻 2001 年 44-55

44立方体箱内自然対流の遷移に関する数値的研究(Shinichiro YANASE)
岡山大学大学院自然科学研究科 岡田隆徳 (Takanori OKADA)
岡山大学大学院自然科学研究科 太田貴憲 (Takanori OHTA)

岡山大学工学部

1.

柳瀬眞一郎

序論

熱対流現象は流体力学, 流体工学, 熱工学にまたがる重要な現象で, 古くから研究されているが, 現在においても理論的に興味深いばかりでなく, 例えば多くの物質の製造工程に対して大きな影響を与える点で, 応用面でも重要な役割を果たしている. その中でも静止流体が温度勾配の上昇によって開始する自然対流の臨界点を求めることは, 特に基本的な問題である. この問題は最初無限平板間の熱対流に対して解かれ, 有名な Chandresekhar の本に述べられている (2). 側壁がある場合も多くの研究があるが, 完全に周囲を境界で取り囲まれた領域での熱対流の臨界点を求めるのは決して容易ではなく, 最近, 全ての境界面で完全熱伝導(温度固定) 条件の下で, 線形安定性によって臨界点が求められた (1).
今回の研究の目的は, エり現実的な側壁が断熱壁の自然対流の臨界点を求め, 超臨界状態での対流の構造を数値計算によって調べようとするものである. スペクトル法は非常に精度が高く , 流れを精密に計算するのには適しているが, 境界の形がたとえ単純でも閉領域のように境界の数が多い場合に適用することは実用上非常に困難である. 一方, ナヴイエ.

$\text{ス}$}$\backslash -$クス方程式を数値的に解くために古くから開発されてきた差分法は, 任意の幾何形状に対して柔軟にコーデイングができるため, 複雑な形状を持つ境界内の流れを計算するのに適している. しかし乱流の数値解析等, 特に高い精度が要求される場合では, 差分法を直接数値計算に利用することに対しては, しばしば精度不足が指摘されることがある. 実は, 層流であっても非常に高い精度の計算が要求される場合があり, それは自然対流などの流れの遷移点付近を計算する場合で, 流れ場がパラメータの微妙な変化によってその振る舞いを大きく変えるため, 時には乱流を計算するのに匹敵する精度が必要となる.
本研究では, スタガード格子上で, 空間 4 次精度の差分法(3) を用いて計算を行った.

最初に, 標準的な検定問題である立方体キャビティー内流れの数値計算を行い, スペクトル法による計算結果 (4) と比較した. 同時に,$\omega-\psi$法により, 高次精度差分法 (空間 6 次精度) で解かれた解(5)とも比較し, その精度を確認した.

本研究で注目した点は, 適当な座標変換を用いて格子点分布を非一様にすることが解の精度向上にどの様に貢献するかという問題である. これに基づいて, 立方体箱内の熱対流を初期値問題として解き, その結果から臨界 Rayleigh 数を推定した. この結果は, 既に得られている側壁で完全熱伝導条件の下で線形安定性によって求められた臨界 Rayleigh 数の値 (1) と矛盾しないものである. 熱伝導方程式を含むナヴイエ.

$7\backslash$}$\backslash -$クス方程式の時間発展の計算には Kim-Moin 法を用いた (6). この方法は比較的簡便で精度が高く, 今回の研究には適当であると考えられる.

45
2.
2$\cdot 1$

基礎方程式

図 1 に示す立方体箱内流れの計算を行う.

基礎方程式

Fig. 1 立方体箱内流れの配置図

重力の方向は一$y$方向とする. 無次元化した連続の式, ナヴイエ. ストークス方程式, およひ熱伝導方程式を以下に示す. なお, 長さは立方体箱の一辺の長さ $L$ で無次元化した. また, 速度はキャビティー流に対しては上壁面の一定ずれ速度$U_{0}$, 熱対流に対しては$\kappa$を温度伝導率として$L/\kappa$で行った (ここでは熱対流に対する基礎方程式のみを示す).
$\nabla\cdot u=0$,$\cdot$.. ..... .... .......... . ... . . ...... . .... .. ..... .. ........... .. . . .. ..........

$(1)$$\frac{\partial u}{\partial t}+(u\cdot\nabla)u$$=-\nabla\tilde{p}+Pr\nabla^{2}u+RaPrTj,$ $\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdot\cdot(2)$$\frac{\partial T}{\partial t}+(u\cdot\nabla)T=\nabla^{2}T$.

ここで, $i,$ $j,$ $k$ はそれぞれ $x,$ $y,$

圧力で, 圧力$p$$\cdot$$z$. .. ... .. ... .... . .. ... ....... ... .... .. . .. ..

方向の単位ベクトル, 流速 $u_{-}=ui$. $+vj+wk,$

$T$

は温度,

$\tilde{p}$(3)は一般化と$\tilde{p}=p+RaPr(y-\frac{1}{2}y^{2})$. ...... .. .. . ....... ......... . .. ........ .... .. . .. .. .. .... ....の関係がある. $Pr=\nu/\kappa$ は Prandtl 数,$\nu$.——–.

$\nabla=\partial/\partial xi+\partial/\partial yj+\partial/\partial zk$

は動粘性係数で, Rayleigh 数 & の定義は後ほど述べる.

$(4)$で, Boussinesq 近似を用いてぃる.

46
方程式 (2), (3) に Kim-Moin 法を適用すると次のようになる.
$\frac{\hat{u}-u^{n}}{\Delta t}=\frac{1}{2}(3\mathrm{H}^{n}-\mathrm{H}^{*-1}’)+\frac{1}{2}Pr\nabla^{2}(\hat{u}+\mathrm{U}^{n})$,........................................... ................. .. .. . . (5)$\frac{u^{n+1}-\hat{u}}{\Delta t}=-\nabla\phi^{n+1}$,$\cdot$........................ .. ....... ........ ............. .... .. .. ..(6)$\frac{T^{1*+1}-T^{||}}{\Delta t}=\frac{1}{2}(3\mathrm{H}_{\mathrm{r}}^{n}-\mathrm{H}_{\tau}^{n-1})+\frac{1}{2}\nabla^{2}(T^{\mathrm{B}}+T^{n-1})$.

................................................................ .. (7)
ここで, $u^{n},$$\mathrm{H}^{n}$$\mathrm{H}_{\tau}^{n}$$\Delta t$$\tilde{p}$’$T^{*}$$=$$=$

はそれぞれ第$n$

時間ステップにおける $u,$ $T$ の値,$-(u^{l}’\cdot\nabla)u^{n}+RaPrT^{\mathfrak{n}}j,$$-(u^{n}\cdot \mathrm{v})T^{\mathfrak{n}},$$\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdot(8)$$\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots(9)$は時間刻み幅を表す. \^u は中間速度 (中間ステツプの非物理的速度) を表す.

$\phi$

は擬圧力と呼ばれ,$\phi$との関係は式 (5), (6) から \^u を消去し式 (2) と比較することで......................... ............ .. .................... ........ (10)$\tilde{p}=\phi-\frac{\Delta t}{2Re}\nabla^{2}\phi$となる. また連続の式は.................................................................. .. .. . (11)$\nabla\cdot u^{n+1}=0$である.

式 (6) の両辺の発散を取り, 連続の式 (11) を考慮するとポアソン方程式$\nabla^{2}\phi^{n+1}=\frac{1}{\Delta t}\nabla\cdot$.

..................................... ..................... .. .. .. . (12)が得られる.

2$\cdot 2$

スタガード格子と座標変換関数$\Delta x,$ $\Delta y,$ $\Delta z$

流体運動の数値計算においては, 対象とする領域を各辺がの格子に分割し, 速度と圧力を異なる点で定義するスタガード格子がよく用いられる. 今回の解析では, これを適用した. スタガード格子を用いると, 圧力の異常な振動を防止することができる (7). さらにボアソン方程式を解く際に, 境界条件としてノイマン条件を用いる場合は, スタガード格子を用いると自然に条件を入れることができることも有利な点である. 流体方程式の数値解法では境界層の分解能を上げるため, 境界近辺の格子間隔を小さくする様な座標変換を行うことがしばしば行われる. 座標変換のための関数として,$h(x)= \frac{1}{2}\{\frac{1}{\alpha}\tanh$

ただし,C。$(2x-1)+1\}\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots(13)$$C_{o}= \frac{1}{2}\log\frac{1+\alpha}{1-\alpha}$

47が広く用いられている.

$\alpha$は $0<\alpha<1$ で, 0 に近づくほど一定の格子間隔に, 1 に近づくほど境界近辺において格子間隔は小さくなる. Fig 2 に $\alpha=0.82,0.92$ ,0.98346 に対する $h(x)$ のグラフを示す.

Fig. 2 幾つかの$\alpha$の値に対する $h(x)$ の関数形

しかし, ここでは示さないが, ポアソン方程式に対して具体的な例について調べた結果によれば, たとえ解の振る舞いが境界付近でなめらかであっても, 解の精度を上げるためには不等間隔格子を取った方が精度が向上することが明らかになった. さらに,$\alpha$の最適値はスキームが 2 次精度かか 4 次精度かにょっても異なることがわかった. 本論文ではスキーム全体の制度を最も高くするため, 4 次精度による流れの計算は $\alpha=0.92$で行い, 2 次精度の計算は $\alpha=0.98346$ でおこなった. $\alpha=0.98346$ の値は平面ポアズイユ流の数値シミュレーションに対する最適値として用いられている値である (8).

3.

キャビティー流への適用

前節までに説明した方法を, 立方体キャビティー内流れの計算に適用した. この流れは, 数値解法の標準的な検定問題であるばかりでなく, ラグランジアンカオスの典型例としても有名である (9).
初期条件として上壁面上 $(y=1)$ で$x$軸方向に無次元速度 $1(u=1)$ をあたえ, それ以外の場所では速度成分 $u$ を 0 とする. また境界条件は, 上壁面で $u=1$ , それ以外の壁面上で $u=0$ である.

本研究では, レイノルズ数$Re=1\mathrm{O}\mathrm{O},$$400$,1000 の 3 つの場合について計算を行った. これらのレイノルズ数に対しては定常解が存在することが知られている (9) ので時間発展方程式の解の漸近状態として定常解を求めた. 収束判定基準は, 各速度成分に対して$u_{res}+v_{res}+w_{re*}<1.0\mathrm{x}10^{-6}$.............. .....................

.....——.

(14)

48..... .. ... .. .... .. ...... ........ . .. ....... .. .. . (15)である. なお $v_{ree},$向,$z$$w_{re\epsilon}$も同様に定義される. ここで $N=N_{l}N_{y}N_{z}$ で,$N_{x},$ $N_{y},$ $N_{z}$はそれぞれ$x$方向,$y$

方方向の格子点数である. 時間刻みはすべての場合において $\Delta t=5\mathrm{x}10^{-3}$ とした.

キャビティー流では, 過去に Ku ら (4) がスペクトル法を用いた高精度な計算を行っていて, 信頼性が高いとされている. また, 西田 (5) は $\omega-\psi$ 法を用いて高い空間精度 (6 次) で計算を行っている. 本研究では, 格, 400, 1000 の場合において子点数を $(30 \mathrm{x}30\mathrm{x}30)$ として $\alpha=0.92$ の計算結果をレイノルズ数 &=1

文献 (4) と, $Re=400$ ,1000 の場合において, 文献 (5) と比較した. 比較は, 立方体の上下, 左右壁面の座標中心を結ぶ直線上での速度分布を用いた.

$\mathrm{v}$$\mathrm{y}$$\mathrm{u}$$\mathrm{x}$(a) 垂直面中心線(b) 水平面中心線

Fig. 3 $Re=1W\mathrm{O}$ に対する速度分布

まず,$Re=1\mathrm{O}\mathrm{O}$のときの速度分布は, 2 次精度においても, 4 次精度においてもスペクトル法と比較してほとんど差が見られず, 正確な解が得られた. $Re=400$ の場合, 水平面中心線 (直線 $x=z=0.5$ ) 上, およひ垂直面中心線 (直線 $y=z=0.5$) 上の速度分布を見ると, スペクトル法と$\omega-\psi$

法め速度分布はほぼ同じであり, 本計算での 4 次精度の速度分布もそれにほとんど重なるような形となった. Re=10 のときの速度分布を Fig.3(a), (b) に示す. これを見ると, レイノルズ数の増加に伴い, 2 次精度解は他の解法と比較して誤差力] 大しているのがわかるが, 4 次精度解,$\omega-\psi$法, スペクトル法の結果は, ほぼ重なり合っており,2 次精度と比べてレイノルズ数の増加に伴う精度の劣化が少ないと思われる.
次に$\alpha$の最適値を求めるため, $\alpha=0.82,0.92,0.98346$ に対する解の垂直面中心線上の速度分布を部分的に拡大して比較した. Fig .4(a) によれば, $\alpha=0.92$ のとき, 格子点数を $(30 \mathrm{x}30\mathrm{x}30)$ から $(40 \mathrm{x}40\mathrm{x}40)$ ,49
さらに $(50 \mathrm{x}50\mathrm{x}50)$ へと増加させると, 速度分布は非常に良い収束を示していることがわかる. 収束解は$\omega-\psi$法による解に近くスペクトル法による解とは少し異なっている. これから, 収束解が正確であるかどうかを断定することは出来ないが, 少なくとも収束性の上で良い性質を持っている. 一方, Fig.4(b) によれ
ば, $\alpha=0.82$ または $\alpha=0.98346$ の場合には, 格子点数を増加させてもあまり良い収束が得られていないことがわかる. この結果から, スタガード格子上において 4 次精度差分法を適用する場合には, $\alpha=0.92$ を採用すると良いことが結論される.

(oe0.82)$—(40^{*}40^{*}40)$$-(30^{*}30^{*}30)$—–$\cdot$$(50^{*}50^{*}50)$$(30*30*30)$$-(40^{*}40^{*}40)$$—(40^{*}40^{*}40)$–—–(a) $\alpha=0.92$(oeO.98346)$-(30^{*}30^{*}30)$$\cdot$$(50^{*}50^{*}50)$—–$\cdot$$(50*50*50)$(b) $\alpha=0.82$ and $\alpha=0.98346$

Fig. 4 $Re=400$ に対する速度分布の拡大図

4.

自然対流の計算

キャビティー流の計算結果を踏まえて, 同じ計算方法で立方体領域内の自然対流の計算をプラントル数$Pr=7.0$ とおいて行った. 立方体領域は下壁面 $(y=0)$ から加熱されるとし, $y=0$ での温度境界条件は$T=\Delta T$

そと固定し, 上壁面 $y=1$ では $T=0$ と固定した. 側面は断熱壁, すなわち法線方向座標を$/\partial n=0$$n$としと仮定した.

Rayleigh 数 $Ra$ は$\Delta T$$Ra=\beta gL^{3}\Delta T/\nu\kappa$で与えられる. ここで$\beta$を用いて... .. . .. .. .. ........... ................... . . .. ........ ... .. .. ..... (16)は流体の体積膨張率;$g$は重力加速度である. 温度場は$\Delta T$で無次元化し, 熱伝導50による平衡温度場と撹乱温度$T=1-y+\tau$.

$\tau$の和に分解した.

. ........ .... ........... .. .... ........ . ........... .. .. ...... .. ...... . (17)$\cdot$初期には流体は静止しているとし, 温度場にのみ $x,$$\tau(t=0)=\tau_{m}y(1-y)$.

$\cdot$$z$方向に一様で$y$方向に変化する撹乱を加えた.

............. .. ..... .... ...... .. ................. ..... ...... . (18)

ここで $\tau_{m}/4$ は初期撹乱温度の領域内の初期時刻における最大値で, 今回の計算では全て$\tau_{m}=\mathrm{O}\mathrm{J}$とした. 時間刻み $\Delta t=2\mathrm{x}10^{-6}$ とし, 格子点数は $(N_{l}, N_{y}, N_{z})=(30,30,30)$ として計算を行った. なお CPUが Alpha 21264 の$5W\mathrm{M}\mathrm{H}\mathrm{z}$の計算機を使用した場合, 1 ステップを計算するのに約 15 秒かかった. & を変化させ, 時間の経過と共に撹乱が発達, あるいは減衰するかを調べた. 流れ場の撹乱の大きさの指標として, 撹乱温度の絶対値と速度場の大きさの領域内の最大値の時間発展を調べた.

Ra=10,$3\alpha 10,50\mathrm{m},$Fig 5, 6 に低 Rayleigh 数$60W$ の場合の結果を示す.

$–_{9}--\cdot$$-^{\mathrm{v}}---\cdot$$\mathrm{t}$$\mathrm{t}$Fig. 5 小さい & に対する${\rm Max}|\tau|$

Fig. 6 小さい & に対するの時間発展$\mathrm{M}\alpha|u|$の時間発展

これらの図では, 時間の経過と共に撹乱温度も速度も共にゼロに近づく様子が見られる. つまり, これらの Rayleigh 数では,自然対流は初期に一時的存在するが発達することはない.

次に, Fig.7, 8 に $Ra=6050$ , $7W0$ ,75,$80\mathrm{N},$$85W$,匍 ,$1\mathrm{m}\infty$

での撹乱の時間発展の様子を示す.

51$\overline{\cong\underline{-}.}$$\mathrm{t}$$\mathrm{t}$Fig. 7 大きい & に対する${\rm Max}|\tau|$

Fig. 8 大きい & に対する ${\rm Max}|u|$ の時間発展の時間発展

これらの図より, $Ra=6050$ を境として撹乱の時間発展の様子が大きく変化していることがわかる.

$Ra<6050$ では撹乱が減衰し, 撹乱温度, 速度場が共にゼロに近づく様子が見られるのに対し, $Ra>6050$

では撹乱温度, 速度場共にたとえ一度減少してもその後また増加し, 最終的には有限の値を保つ. この結果から対流の起こる臨界レイリー数はおよそ$Ra_{\mathrm{c}}\approx 6050$と推定される. なお, この値は Rayleigh 数を 10 刻みで変化させることによって求めた. 側壁が完全熱伝導境界条件では臨界 Rayleigh 数は $Ra_{c}=6798$ であった (1) が, 本研究では側壁が断熱境界条件なので流れがより不安定になりこのような差が生じたと思われる.

0.

–0.

い—クー0.

$\mathrm{t}$Fig. 9$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$に対する${\rm Max}|\tau|$の時間発展

52
次に,$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$の場合の流れの時間変化の様子を, 時間ステップ $t=\Delta t,$ $0.1$ ,0.2134, 0.6 で平面 $z=0.5$に投影される流線を観測することによって調べる. なおこの可視化には AVS を用いた, Fig 10 に点$\mathrm{b}(t=0.1),$$\mathrm{c}(t=0.2134)$, $\mathrm{d}(t=0.6)$ で示される時間での流線を示す.

(a) 流線 $(t=\Delta t)$(b) 流線 $(t=0.1)$(c) 流線 $(t=0.2134)$(d) 流線 $(t=0.6)$

Fig. 10$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$に対する断面 $z=0.5$ での流線$\mathrm{a}(t=\Delta t)$,

53
なお横方向は$x$軸, 縦方向は$y$軸である. Fig 10 より初期時刻$t\approx \mathrm{O}$では直線 $x=0.5$ に関して対称な対流が見られる. 実際これは, Fig 11 (a) のように水平面中心線 $(x=z=0.5)$ に関して対象な渦構造の断面を見ている. しかし, 撹乱の最大値が増加するにっれ, Fig 9 の点$\mathrm{b}$, っまり $t=0.1$ あたりでは既に初期の対称性が崩れ, Fig 10 (b) のように対流の様子が変化している. 対称性が崩れた結果, Fig 10 (d) のような$x-y$ 面で一周する渦の存在する流れになる. 実際これは, Fig 11 (b) のように垂直面中心線 $(x=y=0.5)$を中心に持つ一つの大きな渦の断面であって, これから立方体領域での自然対流の超臨界での平衡状態は全領域を一つのロールが占める様な流れであることがわかった.

(a) 初期(b) 完全発達期

Fig. 11$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$(a) $t=0.64$ における 4 次精度の流線

Fig. 12$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$に対する流線のスケッチ(b)$t=1.66$ における 2 次精度の流線における 4 次精度、 2 次精度の流線の比較

54
最後に自然対流の遷移を数値的に調べるには高精度の計算が必要であることを示すため, 全く同じ初期条件から出発した 2 次精度での計算結果と比較してみた. まず, Fig 12 に$Ra=1\mathrm{O}\mathrm{O}\mathrm{O}\mathrm{O}$に対してほぼ定常状態に達したと考えられる時刻, 4 次精度では t=0.64(図 (a)), 2 次精度では $t=1.66$ (図 (b)) の 3 次元的な流線図を比較した. これらは横方向はでは$z$$x$軸, 縦方向は$z$軸で立方体箱を$y$の正方向から眺めた図で, Fig .12(a)方向の下半分 $(0.5\leq z\leq 1)$ のみの流線を示している. これらの図からわかるように, 4 次精度では流線が直線 $x=y=0.5$ を取り巻く軌跡を描いているのに対し, 2 次精度では直線 $x+z=1,$ $y=0.5$ を取り巻く軌跡を描いている. 一般的に, 熱対流の超臨界状態で成長する渦パターンは, 渦軸が側面に垂直となる傾向があることが知られている (10) ので, 4 次精度の結果がより信頼性が高いと結論される.

次に, 4 次精度と 2 次精度の撹乱温度$\tau$の最大値の時間発展を比較した (Fig 13). $Ra\leq 6000$ では両者はある程度一致しているが, &\geq 6050 となって 4 次精度の結果における臨界点を越えると, 両者の不一致は急激に大きくなる. 一つは臨界 Rayleigh 数の値で, 2 次精度では &。は約 80 である. さらに,4次精度の計算で見られる撹乱のオーバーシュートは, 2 次精度計算では全く見られない. キャビテイー計算では 2 次精度と 4 次精度のこの様な定性的な違いが見られなかったばかりでなく, 2 次精度でもほぼ定量的に満足できる結果が得られたのとは全く異なる様相で, 自然対流の計算が高精度を要求していることを示して
いる.

5.

結論

最初に, 適当な座標変換を伴う空間 4 次精度の差分法によってキャビテイー流の計算を行$\mathrm{a}$$\mathrm{A}$, スペクトル法およひ空間 6 次精度による方法での計算と比べて遜色のない程度に精度の高い解を得ることができた. 計

算精度の向上には激しい変化の起こる境界付近で格子間隔を小さくとる, 不等間隔格子を用いていることが重要であることは言うまでもないが, 座標変換パラメータの最適値は 4 次精度と 2 次精度では異なる. 4 次精度では $\alpha=0.92$ 付近が最適値であり, 2 次精度では $\alpha=0.98346$ 付近が最適値であることがわかった.

$\mathrm{t}$(a)& $\leq 6000$$(\mathrm{b})Ra\geq 8000$

Fig. 13 2 次精度と 4 次精度差分による撹乱の時間発展比較

55
次に, キャビティー流の計算結果から得られた最適パラメータを用いて, 立方体領域での自然対流の数値計算を行った. 初期時刻には水平面中心線 $(x=z=0.5)$ に対して対称な対流が成長するが, 時間の経過と共に, 低 Rayleigh 数では対流は減衰し, 高 Rayleigh 数ではおおよそ垂直面中心線 $(x=y=0.5)$ を中心とする一つの渦が成長し, 平衡状態へ近づく. この変化の発生する臨界 Rayleigh 数は $Ra_{c}\approx 6050$ であることがわかった.

最後に, 空間精度の計算に対する影響を見るために, 2 次精度の計算との比較を行った. その結果, 自然対流の計算には 2 次精度では全く精度が足りないことが示された. 本研究ではただ一つの初期条件に対して撹乱の成長・減衰を調べたため, 今回得られた臨界 Rayleigh 数は真の臨界 Rayleigh 数の上限である. さらに,数値計算による誤差は人工粘性として働き臨界 Rayleigh 数を引き上げる影響があり, どちらにしろ真の臨
界 Rayleigh 数は今回得られた数値より低い可能性が高い. 従って, 臨界 Rayleigh 数のより正確な値を求めるために, さらに精密な計算が必要である.

文献(1) J. Mizushima and O. Matsuda, J. Phys. Soc.

$\mathrm{J}\mathrm{p}\mathrm{n}.,$$66(1997),2337$-2341.

(2) S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (1961) Oxford Univ. Press.
(3) 梶島岳夫, 機論,60-574,B(l994),2058-206.3.
(4) H.-C. Ku,$\mathrm{R}.\mathrm{S}$. Hirsh and$\mathrm{T}.\mathrm{D}$. Taylor, J. Comput. Phys. 70 (1987), 439-462.

(5) 西田秀利, 博士論文第 4 章 (1992).
(6) J. Kim and P. Moin, J. Comput. Phys. 59(1985), 308323.

(7) C.A.J. Fletcher, Computational Techniques for Fluid Dynamics$\mathrm{V}\mathrm{o}\mathrm{l}.2(1991)$, Springer-Verlag.

(8) P.Moin and J. Kim, J.Fluid Mech. 118(1982), 341-377.
(9) 神部勉, 岩津玲磨, 流体のカオス的ふるまい, パリティ別冊 No .10(1994), 79-82.
(10) J.P.Gollub, A.R.McCarriar and J.F.Steinman, J.Fluid Mech. 125 (1982), 259281.