Silviase/jgov_v1.4_full_ann huggingface dataset:
$$ \theta=\left(\mathrm{T}_{\mathrm{r}}-\mathrm{T}_{0}\right)/\left(\mathrm{T}_{\mathrm{h}}-\mathrm{T}_{0}\right) $$

(6) 式を満足する無次元流れ関数Φを用い、無次元速度 U、V を記述する。

$$ \mathrm{U}=\partial\Phi/\partial\mathrm{Y} $$

$$ \mathrm{V}=-\partial\Phi/\partial\mathrm{X} $$

(15)、(16)式から無次元渦度Ωは次式で記述される。

$$ \Omega=-\left(\partial^{2}/\partial\mathrm{X}^{2}+\partial^{2}/\partial\mathrm{Y}^{2}\right)\Phi $$

(7) 式の密度差の項を (5) 式で置き換え、(7) 式を無次元渦度関数 Ω で記述する。

$$ \partial\Omega/\partial\mathrm{T}+\partial\left(\mathrm{U}\Omega\right)/\partial\mathrm{X}+\partial\left(\mathrm{V}\Omega\right)/\partial\mathrm{Y}= $$

$$ \mathrm{Gr}\partial\theta/\partial\mathrm{X}+\left(\partial^{2}/\partial\mathrm{X}^{2}+\partial^{2}/\partial\mathrm{Y}^{2}\right)\Omega $$

$$ \partial\theta/\partial\mathrm{T}+\partial\left(\mathrm{U}\theta\right)/\partial\mathrm{X}+\partial\left(\mathrm{V}\theta\right)/\partial\mathrm{Y}= $$

$$ \mathrm{P r}^{-1}\left(\partial^{2}/\partial\mathrm{X}^{2}+\partial^{2}/\partial\mathrm{Y}^{2}\right)\theta $$

$$ \mathrm{Gr}=\alpha\mathrm{g}\mathrm{\boldmath~a}^{3}\left(\mathrm{T}_{\mathrm{h}}-\mathrm{T}_{0}\right)/\nu^{2} $$

$$ \mathrm{P r}=\left(\mu/\rho_{0}\right)/\kappa $$

### 2.3 差分化と数値計算

区画の高さと幅を等しいとし、熱源の幅を区画幅の1/5にした。Fig.2に示すようなN×Mの空間要素に分割し、初期条件と境界条件を与えた。軸対称を仮定して計算を行うことが計算機資源節約の目的で行われているが、ここでは軸対称を仮定せずに計算を行って対称性を

Fig.2 Schematic of nondimensional representation of the flow field

調べた。今回の計算では垂直方向と水平方向の区間分割数と空間刻み幅は同一に設定した。

区画の中心 P では、熱源からの気流の垂直方向速度の時間変化と熱源との温度差との関係を調べた。区画中心軸 I=M/2 を対称軸として速度の非対称性を空間対称位置の速度差の絶対値を速度の絶対値で規格化、二乘積算し流れ場の非対称性指数に用いた。

(18)、(19)式の空間微分を中心差分、時間微分を前進差分で置き換えた差分方程式を用い計算を行った。(18)、(19)式の左辺第2項と第3項には、速度と渦度、速度と温度の積の微分が含まれている。そのため速度の小さい領域では空間微分を中心差分で置き換えることは問題がないが速度が大きくなり運動量や熱の拡散速度を上回る時には下流側の条件が物理的に計算点に作用しないにもかかわらず中心差分で置き換えることにより運動量や熱を輸送する形になる。この非現実的な差分化を回避するためには、計算時に流れの方向を調べ、上流側の条件を反映し計算を行った。このため、速度零では、厳密解と一致するが速度の絶対値が増加すると厳密解との差が増加する。さらに速度の絶対値が十分大きくなると厳密解との差が10%程度 $ ^{10)} $ になることが知られている。

漸進加速緩和法（successive overrelaxation method）を用い、渦度分布から流れ関数を求めた。計算の安定性を考慮し、Lax-Richtmyerの安定条件 $ ^{5)} $ を満たす最大時間刻みのおよそ90%を超えないように時間刻み幅を調整した。Lax-Richtmyerの安定条件は、新たな時間における渦度と温度の値を過去の渦度と温度から計算する過程で過去の渦度と温度の線型結合の係数の符号がすべて正でなければならない。この計算では、温度の計算における制限を満たすことが最も厳しい条件となる。無次元の時間刻み幅を無次元空間刻み幅を無次元速度の最大値で割った値以下に保つことで達成される $ ^{5)} $ 。

計算開始時には、最大速度がグラスホフ数の平方根の0.3倍に比例することと開始時の速度は、最大値の0.1倍程度であることが過去の研究 $ ^{5)} $ で示されているので、無次元空間刻み1/300をこの最大速度で除したグラスホフ数の平方根を基準に計算を行い、計算開始時の速度は、最大値の0.01倍よりも小さいことを確認し、計算を迅速に行う目的で無次元の時間刻み幅をグラスホフ数の平方根の10倍とした。計算実行時には無次元速度の最大値を用い、無次元時間刻み幅を調整した。

数値計算には4台のワークステーションを使用し、3種類のフォートランコンパイラーで同一のソースプログラムから実行プログラムを作成した。代表的な初期条件と境界条件で同一の結果が得られるかを異なるワークステーション上で確認を行った。
### 3. 数值计算結果

#### 3.1 グラスホフ数 $ 1.0 \times 10^{6} $ の計算結果

グラスホフ数を $ 1.0 \times 10^{6} $ の条件で区画高さと幅をそれぞれ 300 分割し、計算した結果を Fig.3 から Fig.6 に示す。Fig.3 では無次元時間 $ T = 1.49 \times 10^{-2} $ の区画内の無次元温度 $ \theta $ 、水平方向無次元速度 U、垂直方向無次元速度 V の分布を濃淡で示してある。無次元変数が零のとき灰色で表示し、値が正の方向へ増加すると白く表示され、負の方向へ減少すると黒く表示される。無次元温度は、ほとんどの部分が零に近いので零から 0.1 の範囲を 0-1 に強調表示してある。無次元速度は、グラスホフ数の平方根程度の大きさなので速度の方向が変化する零近傍を表示し、無次元速度 -8 から +8 の範囲を濃淡の -1 から +1 で表現し、-8 以下や 8 以上はそれぞれ濃淡の -1 と 1 に丸めて表示した。Fig.3 に水平方向無次元速度 U の方向を矢印で示した。この図では水平方向無次元速度 U は、中心軸に対して速度の符号が交代する奇関数である。垂直方向無次元速度 V は、区画の垂直軸に対して速度の符号は交代しない偶関数である。

過去の数値計算結果 $ ^{5)} $ と比較する計算座標系が異なる

Fig.3 Nondimensional temperature, vertical and horizontal velocities

Fig.4 Vertical velocity $ Gr = 1.00 \times 10^{6} $

ものの無次元温度分布は良く似ている。Fig.4は、壁面から最も遠く熱気流の中央に位置する点Pにおける垂直方向無次元速度Vの時間変化である。熱源による加熱を受け、垂直方向無次元速度Vが増加し、 $ T=7.0\times10^{-3} $ 付近で最大値となり、その後、収束している。

Fig.5に点Pでの熱源との温度差DTと垂直方向無次元速度Vを示す。温度勾配が対流の駆動力の源なので温度差が縮小するにつれて無次元速度Vが増加する $ ^{10} $ 。温度差が最小値になった後は、温度差が拡大しながら無次元速度Vが増加し、無次元速度Vが最大値から少し減少した点に収束している。

Fig.6 に区画内の対称位置での速度差の二乗和を速度の二乗和で規格化した非対称度合いの時間変化を示す。完全に対称であれば非対称度合いは零である。計算開始して熱源から高温部が成長するにつれ非対称度合いは減少している。

#### 3.2 グラスホフ数 $ 1.0 \times 10^{8} $ 付近の計算結果

グラスホフ数を $ 1.0 \times 10^{6} $ の条件で区画高さと幅をそれぞれ 300 分割した場合には、過去の研究で得られた温度分布を再現できたのでグラスホフ数をさらに増加させ計

Fig.5 Vertical velocity and temperature gradient ( $ Gr = 1.00 \times 10^{6} $ )

Fig.6 Asymmetric degree $ (\mathrm{Gr}=1.00\times10^{6}) $
算を行った。グラスホフ数を増加させることは、この研究では区画の高さaを増加させた状況を示している。

区画高さと幅をそれぞれ300分割し、計算した結果をFig.7からFig.15に示す。グラスホフ数は $ 1.27\times10^{8} $ である。無次元時間 $ 4.0\times10^{-3} $ は、代表長さ3mで2400秒（40分）、30mで240000秒（16.7日）である。代表長さ3mでは、Fig.15加熱開始から40分後、代表長さ30mでは、同様

Fig.7 Nondimensional temperature, vertical and horizontal velocities

Fig.8 Nondimensional temperature, vertical and horizontal velocities

Fig.9 Nondimensional temperature, vertical and horizontal velocities

に加熱開始から16.7日後である。

Fig.7の無次元温度θと水平方向無次元速度Uを比べると上昇するドーム状高温部の上下では周囲から空気が中心に向かって流入している。ドーム状高温部の高さでは周囲に向かって空気が流出している。無次元温度θと垂直方向無次元速度Vを比べると上昇するドーム状高温部周辺は上に動いているがその直上部では、下に動いている。Fig.7のドーム状高温部の両側に二股に枝分かれした流れができているが、この部分はドーム状高温部中心部上部に形成されたよどみ点から両側へ拡大する流れと高温部の境界に形成された渦によって形成されたことが水平方向速度成分Uの方向からわかる。

Fig.8 に区画内全域の気体が運動している無次元時間 $ T=8.81 \times 10^{-4} $ の様子を示す。水平方向無次元速度 U は、ドーム状高温部の構造に対応して分布が作られている。垂直壁面近傍には、灰色の速度が零に近い領域が形成されている。境界層の厚みは左右に広がった高温部の高さ近傍で最小となっている。区画の隅角部では厚くなっている。垂直方向無次元速度 V はドーム状高温部の近傍では上向きである。左右に広がった高温部下側では、下向きとなっている。

Fig.9 に熱源とドーム状高温部間柱状の部分の上半分が細くなっている T=1.03×10^{-3} の様子を示す。Fig.10 の T=1.50×10^{-3} の様子を見ると高温部は、区画天井に達している。柱状の部分の太さが中央部で太くなり、その上下で細くなっている。水平方向無次元速度 U は、この太くなった柱状の部分で著しく歪んでいる。

Fig.11 に $ T=1.62\times10^{-3} $ の様子を示す。水平方向無次元速度 U の分布を見ると中心軸に対する分布は奇関数からは大きく異なる状況となっており、柱状高温部上部で大きく歪んでいる。垂直方向無次元速度 V の分布は天井近傍を除くと歪みは小さい。

Fig.12 に $ T=1.68\times10^{-3} $ の様子を示す。天井直下に島状

Fig.10 Nondimensional temperature, vertical and horizontal velocities
高温部が形成されている。Fig.13に $ T=1.88\times10^{-3} $ の様子を示す。無次元温度θの分布を見ると天井直下で左右の非対称性が大きい。中心と左壁面の中間に高温部分が明確に見えるが、右壁面の中間部には比較的温度の低い灰色部分が拡がっている。Fig.14に $ T=3.68\times10^{-3} $ の様子を示す。無次元温度θの分布は均一化が進行しているが、中央部の柱状高温部は緩やかに湾曲している。垂直方

Fig.11 Nondimensional temperature, vertical and horizontal velocities

Fig.12 Nondimensional temperature, vertical and horizontal velocities

Fig.13 Nondimensional temperature, vertical and horizontal velocities

向無次元速度 V の垂直中心軸に対しての対称性は Fig.13 に比べ低い。

Fig.15 に $ T=4.01\times10^{-3} $ の様子を示す。無次元温度 θ の分布を見ると垂直中心軸に対して対称性が高く見える。ところが無次元速度 V の分布を見ると垂直中心軸に対する対称性は低い。

Fig.16 に区画の中央 Fig.2 の点 P での垂直方向無次元

Fig.14 Nondimensional temperature, vertical and horizontal velocities

Fig.15 Nondimensional temperature, vertical and horizontal velocities

Fig.16 Vertical velocity $ \left(\mathrm{Gr}=1.27\times10^{8}\right) $
速度 V の時間変化を示す。Fig.7 の $ T=6.26\times10^{-4} $ の時点では、高温部は区画の下半分にあり点 P で観測された速度は小さい。Fig.8 の $ T=8.81\times10^{-4} $ の時点では高温部が点 P に到達し 4000 程度の無次元速度 V となっている。その後、3300 から 4500 程度の無次元速度 V で $ T=2.5\times10^{-3} $ までの期間変動を繰り返している。その後、無次元速度 V は $ 5\times10^{-4} $ 程度の周期で変動している。

Fig.17 Vertical velocity and temperature gradient $ \left(\mathrm{Gr}=1.27\times10^{8}\right) $

Fig.18 Asymmetric degree $ Gr = 1.27 \times 10^{8} $

Fig.19 Vertical velocity (Gr = 1.35x10^{8})

Fig.17に点Pでの熱源との温度差DTと垂直方向無次元速度Vを示す。温度勾配が対流の駆動力の源なので温度差が縮小するにつれて無次元速度Vが増加する $ ^{11)} $ 。温度差0.8から0.92の範囲で右下がりの曲線軌道を作る。この図を見るとDT=1.0から始まる初期の軌道の周囲DT=0.85を中心に（DT、V）の軌跡は振動しているように見える。

Fig.18 に区画内の対称位置での速度差の二乗和を速度の二乗和で規格化した非対称度合いの時間変化を示す。完全に対称であれば非対称度合いは零である。計算開始して熱源から高温部が成長するにつれ非対称度合いは減少し、Fig.7 の $ T=8.81\times10^{-3} $ の時に最小となっている。Fig.10 の $ T=1.50\times10^{-3} $ 以降急激に非対称度合いが増加している。 $ T=2.20\times10^{-3} $ までは一時減少するがその後再び増加してゆく。これは、垂直方向無次元速度の短周期の変動が収束する時点とはぼ一致している。

Fig.19 から Fig.21 にグラスホフ数を $ 1.35 \times 10^{8} $ に増加させた場合の時間変化を示す。Fig.19 の点 P での垂直方向無次元速度 V の時間変化を見ると $ T = 2.5 \times 10^{-3} $ までの挙動は類似しているが、その後無次元速度が 2500 程度に低下している。Fig.20 に点 P での熱源との温度差 DT と垂直方向無次元速度 V を示す。Fig.17 と比べると垂直方向無次元速度 V が小さく熱源との温度差の大きい位置へ無次元速度 $ T = 2.5 \times 10^{-3} $ 以降移動している。Fig.21 と Fig.18 を比較すると急速に非対称度合いが増加するのがわかる。

Fig.22からFig.24にグラスホフ数を $ 2.00\times10^{8} $ に増加させた場合の時間変化を示す。Fig.22の点Pでの垂直方向無次元速度Vの時間変化を見ると $ T=8.0\times10^{-3} $ 以降は無次元速度4800程度で安定している。Fig.23に点Pでの熱源との温度差DTと垂直方向無次元速度Vを示す。Fig.17と比べると単純な軌跡となっている。Fig.24とFig.18を比較すると非対称度合いが小さく、変動も少な

Fig.20 Vertical velocity and temperature gradient $ \left(\mathrm{Gr}=1.35\times10^{8}\right) $
い。

Fig.25からFig.27にグラスホフ数を1.00×10⁸に減少させた場合の時間変化を示す。Fig.16とFig.25を比べると点Pでの垂直方向無次元速度Vの時間変化は小さく短時間で収束しているように見える。Fig.26をFig.17と比べると熱源との温度差DTと垂直方向無次元速度Vともに変動が小さいのがわかる。Fig.27の非対称度合いも

Fig.21 Asymmetric degree (Gr = 1.35x10^{8})

Fig.22 Vertical velocity $ \left(\mathrm{Gr}=2.00\times10^{8}\right) $

Fig.23 Vertical velocity and temperature gradient $ \left(\mathrm{Gr}=2.00\times10^{8}\right) $

$ T=3.0\times10^{-3} $ までは小さい。

Fig.28 にグラスホフ数を $ 1.00 \times 10^{8} $ 、区画高さと幅をそれぞれ 100 分割した場合の点 P での垂直方向無次元速度 V の時間変化を示す。Fig.25 と比較すると速度変動が大きくなっている。Fig.29 にグラスホフ数を $ 1.00 \times 10^{8} $ 、区画高さと幅をそれぞれ 50 分割した場合の点 P での垂直方向無次元速度 V の時間変化を示す。時間刻みの許

Fig.24 Asymmetric degree (Gr = 2.00x10^{8})

Fig.25 Vertical velocity $ \left(\mathrm{Gr}=1.00\times10^{8}\right) $

Fig.26 Vertical velocity and temperature gradient $ \left(\mathrm{Gr}=1.00\times10^{8}\right) $
容幅が大きくなっているので長時間の計算結果となっている。Fig.25やFig.28と比較すると初期に無次元速度4500程度のオーバーシュートが生じている。

### 3.3 グラスホフ数と垂直方向無次元速度の関係

Fig.30 に区画高さと幅をそれぞれ 300 分割した場合グラスホフ数を変化させた場合の垂直方向無次元速度 V の変化を示す。図中の破線は、グラスホフ数の平方根に 0.5

Fig.27 Asymmetric degree $ Gr = 1.00 \times 10^{8} $

Fig.28 Vertical velocity with reduced grid N=100 (Gr=1.00x10^{8})

Fig.29 Vertical velocity with reduced grid N=50 (Gr=1.00x10^{8})

を掛けた垂直方向無次元速度の近似曲線である。グラスホフ数が $ 1.00\times10^{6} $ 以下の領域では、粘性力が浮力に比べ無視できないために近似曲線よりも低い無次元速度となっている $ ^{5,8)} $ 。グラスホフ数が $ 1.00\times10^{7} $ 以上の領域では、近似曲線との間隔がグラスホフ数の増加にしたがって変動を繰り返す。

Fig.31 に区画高さと幅をそれぞれ 100 分割した場合グ

Fig.30 Vertical velocity with Grashof number N=300

Fig.31 Vertical velocity with Grashof number N=100

Fig.32 Vertical velocity with Grashof number N=50
ラスホフ数を変化させた場合の垂直方向無次元速度 V の変化を示す。グラスホフ数が $ 1.00 \times 10^{6} $ 以下の領域では、Fig.30と同様に近似曲線よりも低い無次元速度となっている。分割数が Fig.30より小さくなっているため粘性力の影響が小さくなり近似曲線との差は小さくなっている $ ^{5,8)} $ 。Fig.30と比べるとグラスホフ数が $ 1.00 \times 10^{7} $ 以上の領域では、近似曲線との間隔のグラスホフ数増加にもなう変動は小さくなっている。

Fig.32 に区画高さと幅をそれぞれ 50 分割した場合グラスホフ数を変化させた場合の垂直方向無次元速度 V の変化を示す。Fig.31 と比べるとグラスホフ数が $ 1.00 \times 10^{6} $ 以下の領域では、ほぼ同一の挙動である。 $ 1.00 \times 10^{8} $ 付近の領域では、Fig.31 よりも近似曲線に近くなっている。ところが $ 1.00 \times 10^{10} $ 以上の領域では、Fig.31 よりも近似曲線との差が大きくなっている。

Fig.33 に区画高さと幅をそれぞれ 512 分割し、グラスホフ数を $ 2.0 \times 10^{8} $ に選び計算を行った時の点 P における垂直方向無次元速度 V の時間変化を示す。300 分割の場合には、無次元速度 V の変動は小さかったが分割数を増加させると再び変動が増加している。無次元速度 V の平均値はほぼ同じ 4800 程度である。

### 4. 考察

代表長さが大きな空間内の熱対流を数値計算で求める場合、高グラスホフ数の条件で計算を行うことになる。代表長さが大きい場合、特性時間が代表長さの二乗に比例して増大する。数値計算で定常になった流れ場は、通常の火災では想定できないほど長時間継続しなければならない。今回計算を行ったグラスホフ数1.0×10⁸程度の領域で代表長さを30mとした場合、10日程度経過しても流れ場の変動が見られる。代表長さが大きな空間内の熱対流を数値計算で求める場合、定常解よりも非定常解を用いて実時間程度の流れ場を検討する必要がある。

Fig.33 Vertical velocity $ (\mathrm{Gr}=2.00\times10^{8},\mathrm{N}=512) $

グラスホフ数が $ 1.0 \times 10^{6} $ 程度であれば空間を 300 分割した計算を行っても緩やかに収束する。この条件であれば速度変動が小さくなったことを収束判定にすることが可能である。

速度変動の大きさや平均速度からの差を収束判定の条件に使用する場合、グラスホフ数を $ 1.0\times10^{8} $ 程度にすると20%程度の速度変動を許容しなければ収束と判断されないことになる。

軸対称の条件で解析を行った研究 $ ^{2,5)} $ と比較するとグラスホフ数 $ 1.0\times10^{6} $ では対称な温度場と流れ場が観測されたが、グラスホフ数 $ 1.0\times10^{8} $ 付近では非対称な温度場と速度場の変位が観測された。本来、支配方程式、境界条件、初期条件の何れもが中心軸に対して対称性を有していれば、数値計算結果も対称となるはずである。数値解析法の数値誤差の微小な非対称性から大規模な非対称性が発生していると推定されるが、数値誤差の計算過程での積算解析は行わなかった。数値計算では、均一な初期条件を与えることは容易であるが現実の問題では均一な初期条件を与えることは容易ではない。数値解析法の影響と初期条件のゆらぎについては今後別に検討することにした。

なお、類似の二次元熱伝達現象を解析的に検討した結果によると火源と壁面までの距離と区画高さの比を2とした条件でレイレー数2.0×10⁶以上で伝ば型の擾乱が成長する $ ^{12)} $ となっている。プラントル数0.72の空気の場合、グラスホフ数は、2.7×10⁶となる。今回の数値計算結果を見ると伝ば型の擾乱が成長する場合には、定在波となるまでは非対称の擾乱が存在すると考えられる。解析は、変位が小さい条件で行われているので今回の数値計算で見られた変動のどの程度が伝ば型の擾乱の寄与であるかは今後の検討が必要である。また、この解析では温度場からエネルギーが擾乱に供給される場合、位置エネルギーから擾乱に供給される場合、主流から擾乱にエネルギーが供給される場合があるとしており、今回の数値計算結果に見られるように温度場と速度場が変化している場合には擾乱の成長も変化するとおもわれる。

グラスホフ数 $ 1.0 \times 10^{8} $ 付近では数値計算を行って初期の擾乱が収束し、一時的に非対称性が小さくなる現象が観察された。ところが時間の経過と共に非対称性が増大する現象も観察された。非対称性が小さくなる段階と増大する段階で異なっているのは流れ場の速度である。このことから数値計算の収束条件とは別に摂乱の成長に関するグラスホフ数と速度上限が存在すると推定できる。今回の数値計算ではどのような形の速度上限が存在するかについては明らかではない。

グラスホフ数 $ 1.0 \times 10^{8} $ 付近で軸対称の条件では生成しないはずの非対称な熱気流の蛇行が観察されたこと
は、軸対称の条件を導入することにより蛇行現象が別な現象として表現される可能性が示された。

区画高さと幅を分割するときの分割数を変化させた場合、過去の研究で指摘された $ ^{5,8)} $ 偽粘性と偽熱拡散の影響と思われる流れ場の過剰な安定化は、分割数を300とした数値計算でグラスホフ数が $ 1.00\times10^{6} $ 以下の領域で明確に垂直方向無次元速度の低下として観察された。グラスホフ数がさらに大きくなった $ 1.00\times10^{10} $ 以上の領域では、分割数を300とした場合よりも分割数を50とした場合の方が垂直方向無次元速度の低下は大きかった。分割数が少ない場合には、速度変動成分のうち空間周波数の高い成分が平均化されていることになる。分割数を512として分割数300で垂直方向無次元速度Vの変動の小さくなるグラスホフ数 $ 2.0\times10^{8} $ の条件で計算を行うと無次元速度Vの変動が大きくなった。空間分割の細分化による偽粘性や偽熱拡散による分散と基礎方程式の非線形項とのつりあい $ ^{13)} $ により一種の波動が伝ばしていると考えられる。グラスホフ数 $ 1.0\times10^{8} $ 付近で観測された無次元速度Vの変動も同様なつりあいによる一種の波動によると考えられる。数値計算で観測された波動の原因となる偽粘性と今回モデル化で無視した粘性項の大きさについて評価を行っていないので、この波動が実験的に観測されるかについては現在不明である。速度変動成分を考慮した乱流モデルを導入する場合や速度変動を考慮する場合には空間分割の影響を評価する必要が確認できた。

熱対流現象の理論解析では、グラスホフ数と不安定化する渦スケールの関係が与えられている $ ^{14} $ 。ところが有限の流れ場の乱れと温度場の乱れが成長した段階でのグラスホフ数と不安定化する渦スケールの関係は数値計算でしか与えられない $ ^{2} $ 。理論解析の結果を見ると高グラスホフ数では波数の高い成分に不安定な領域が拡大するので少なくとも数値計算の初期段階では波数の高い成分がある程度成長すると考えられる。このような段階で、波数の高い成分を平均化しながら計算を行うと速度変動の数値計算結果等が見かけ上平滑化したり、変動したりする原因となると推定される。

熱対流の問題では、決定論的非周期的流れ場 $ ^{11)} $ が知られているので支配方程式と数値計算の特性によっては見かけの異なる結果が同一の初期条件から見かけ上ランダムな挙動が生じる。今回の計算を行った分割数を300、グラスホフ数 $ 1.0\times10^{8} $ 程度条件でグラスホフ数を10%程度変化させた場合、点Pでの熱源との温度差DTと垂直方向無次元速度Vの描く軌道は敏感に変化した。このグラフで用いた温度差と垂直方向速度は、決定論的非周期的流れ場 $ ^{11)} $ の研究でも流れ場の特性パラメータとして用いられている。今回の数値計算では、空間分割数の影響も受けているので確定できないが、駆動力となる温度差と垂直方向速度は単純な釣り合いではなく、流れ場の特性に応じた準安定点の周囲で変動しているようにも見える。

### 5 \. 結論

古典的な数値計算法を用い単純な二次元正方形の空間を対象に数値計算を行った。数値計算の安定性を調べる目的でグラスホフ数を増加させ、流れ場を計算した。空間分割数300、グラスホフ数 $ 1.0\times10^{8} $ 程度の条件でグラスホフ数に敏感な結果が得られた。空間分割数を減らし数値計算を行い空間分割数の影響を低グラスホフ数側で観察したが、高グラスホフ数側では不明確であった。

## 参考文献

1) D. Makarov, O. Sugawa, A. Karpov and T. Tsuruda: "CFD Modelling of Garage Fire and Gaseous Fire Suppression", Fire Science & Technology 20-1, pp.37-49 (2000)

2) C.-Y. Chow: "An Introduction to Computational Fluid Mechanics", pp.359-374, JOHN WILEY&SONS (1979)

3）水谷幸夫・香月正司共訳 スハス V.パタンカー："コンピュータによる熱移動と流れの数値解析"、pp.11-41、森北出版株式会社（1985）

4）高洲 達監訳 J. ウオルシュ編："数値解析概論〝、pp.1-171、日本評論社（1970）

5) K. E. TORRANCE and J. A. ROCKETT: "Numerical study of natural convection in an enclosure with localized heating from below-creeping flow to the onset of laminar instability", J. Fluid Mech. 36-1, pp.33-54 (1969)

6) T. HANDA, M. MORITA, K. HAYASHI, and T. ISHII: "COMPUTER SIMULATION OF THE MOTION OF HEAT FLOW INDUCED BY THE FIRE IN THE LONG-SCALE TUNNEL", Fire Science & Technology 1, pp.91-110 (1981)

7) K. SATOH: "A NUMERICAL FINITE DIFFERENCE STUDY AND EXPERIMENTAL VERIFICATION OF AIR ENTRAINMENT BY FIRE PLUME BETWEEN TWO HORIZONTAL PARALLEL PLATES", Fire Science & Technology 2-2, pp.127-136 (1982)

8) J. W. DEARDORFF: "A Numerical Study of Two-Dimensional Parallel-Plate Convection", J. Atmos. Sci. 21, pp.419-438 (1964)

9) 高野 暲："流体力学〝、pp.101-156、岩波書店(1975)

10) Anil W. Date: "Introduction to Computational Fluid Dynamics", pp.59–61, CAMBRIDGE (2005)
11) E. N. LORENZ: "Deterministic Nonperiodic Flow", J. Atmos. Sci. 20, pp.130-141 (1963)

12) A. E. GILL and A. DAVEY: "Instabilities of a buoyancy-driven system", J. Fluid Mech. 35-4, pp.775-798 (1969)

13）日本流体力学学会編："流体力学の世界〟、pp.38-39、朝倉書店（1990）

14) S. Chandrasekhar: "Hydrodynamic and Hydromagnetic Stability", pp.1-75, DOBER PUBLICATIONS, INC. (1961)
# A Study on Laminar Numerical Simulation of an Enclosure Fire

(Abstract)

Takashi Tsuruda

(Accepted February 17, 2010)

Numerical techniques are used to analyze fire phenomena for developing extinguishing systems. With an appropriate set of initial and boundary conditions, a unique solution is obtained with the help of numerical techniques. For analyzing fire phenomena, the computational stability of the numerical scheme should be examined as well as the solution of the steady state. In this study, a series of numerical computations has been carried out using a classical numerical scheme for various values of Grashof number. During this series of numerical computations, several cases show that the no dimensional vertical velocity at the center of the flow field is strongly dependent on the value of the Grashof number around $ 1.0 \times 10^{8} $ with 300 grids. For the given boundary conditions, the solution behaves like an initial value problem. The effect of the grid size on the no dimensional vertical velocity at the center of the flow field is observed with lower Grashof number.
論文

# 微細水粒子による対向流拡散火炎の消炎

内藤浩由 * 萩原由梨 ** 上戸徹 ** 佐宗祐子 *** 小谷良信 ** 吉田亮 **

(平成22年2月22日受理)

従来から水は火災の消火に使われてきたが、火災形態に対する水の必要量が科学的に分析された例は少なく、経験的に必要以上な大量の水を散布するため、大部分の水が蒸発しきれず、火災現場等に残り、これが二次的被害をもたらしている。また現在でも、特定防火対象物に設置されているスプリンクラー等は主に経験的に設計されており、効率的な消火、または水損防止の観点からも、科学的根拠のある最適設計が望まれている。

微細化された水粒子による拡散火炎の消炎現象を支配する物理的な要因としては、運動量、物質および熱の輸送現象と水の相変化（蒸発）が考えられ、また化学的要因としては、化学反応速度を考慮する必要がある。微細水粒子による消火研究は、これまで数多く行われているが、しかしながら、様々なパラメータが複合的に関与するため、微細水粒子による拡散火炎の消炎現象の物理的・熱的・化学的機構の詳細はいまだに不明な点が多い。

そこで本研究では、従来必ずしも明確ではなかった微細水粒子の特性を正確に把握すると同時に、それを制御する技術を確立する事に特徴がある。また、流体力学的特性が明確に定義される対向流拡散火炎を使用し、この様に微細水粒子による消炎機構を定量化する事に独創性がある。

微細水粒子の粒径、数密度および粒径分布が消炎におよぼす影響を実験的に検討した結果、微細水粒子の火炎帯滞留時間と水粒子消失時間の考察より、各よどみ速度勾配に対する消失粒径を導き、対向流拡散火炎の消炎に寄与する微細水粒子の最適粒径を示した。

また PDPA (位相ドップラー式レーザ粒子分析計) で計測された、火炎帯近傍もしくは火炎帯中における水粒子径の分布傾向は、微細水粒子の火炎帯滞留時間と水粒子消失時間の考察結果と良い一致を示した。更に、微細水粒子による消炎特性は、水粒子の蒸発面積の増加に伴い消炎時のよどみ速度勾配が減少する事が明らかとなった。

### 1. はじめに

水は熱容量と蒸発潜熱が大きいため、火炎から吸収する熱量は大きい。また、水が蒸発する時、比容積は約1700倍になり、その結果、火炎近傍において酸素は希釈される。さらに環境負荷が小さいため、火災抑制に対して水は好ましい物理的性質を持つ。水粒径を微細にすれば、表面積が十分に大きくなり、熱の吸収や蒸発が促進され、火災抑制に対する水の効果が増加する事が期待されるが、微細水粒子による消炎機構の詳細については不明な点が多い。本研究で使用した対向流拡散火炎は、よどみ速度勾配を制御する事が出来るため、よどみ速度勾配による消炎機構と微細水粒子による消炎機構を分ける事が出来る。本研究では伸長を受けた拡散火炎の微細水粒子による消炎機構を、物理的影響と熱的影響を考慮し定量的に明らかにするため、位相ドップラー式レーザ粒子分析計（Phase Doppler Particle Analyzer: 以下PDPA）を用いて微細水粒子特性を定量化し、微細水粒子の粒子特性が消炎効果におよぼす影響を詳細に検討する事を目的とする。

### 2. 実験装置および方法

### 2.1 対向流拡散火炎バーナー及びPDPAシステム

Fig.1に示すように、空気はロータリーブロワーからサージタンクで脈動を抑えた後に、燃焼風洞に供給され、流量はブロワー回転数で制御される。燃焼風洞は、広開角拡散筒、整流筒、縮流ノズルからなり、縮流ノズルにより出口流速が一様化され、燃焼室へ導かれる。多孔質円筒バーナーは外径30[mm]の焼結金属でできており、燃焼室中心軸上、縮流ノズル出口下流側100[mm]に設置される。燃料に純度99.8%のCH_{4}を用い、流量はマスフローコントローラーにより制御される。燃料は多孔質円筒バーナーから均一に噴射され、空気流との間に対
向流拡散火炎が形成される。後流への燃料噴射を防ぐため、多孔質円筒バーナーの下流側半分に、耐熱塗料が塗布されている。燃焼室断面は、30 [mm]×120 [mm] の矩形であり、燃焼室の両側は、火炎可視化および光学計測のため、石英のガラス窓2枚で構成されている。多孔質円筒バーナーを設置するため、一方の石英ガラス窓には、穴があけられている。

また PDPA (位相ドップラー式レーザ粒子分析計) を使用し、水粒子速度および粒径計測を行なった。光源は、最大出力 4 [W] のアルゴンイオンレーザーを使用し、分光器で 514.5 [nm] の波長帯へ分光し、ビーム間隔を 20 [mm]、焦点距離 500 [mm] で計測を行なった。また、レシーバーにより受光されたバースト信号は、リアルタイム信号処理器により高速離散フーリエ変換が行なわれ、粒子流速および各バースト信号の位相差により粒径計測が行なわれる。PDPA 計測位置は、Fig.1 の座標系へ示す様に、多孔質円筒バーナー表面を原点 0 とし、上流側を y 方向、半径方向を x とし、微細水粒子の粒径および流速は、x=0 [mm] の条件のもと y 方向のよどみ流線に対して計測を行った。

### 2.2 微細水粒子供給系

Fig.1に示すように、広開角拡散筒の上流側に、微細水粒子供給系を設置し、二流体噴射弁および超音波霧化装置を設置し、ブロワーから供給される空気流と均一混合に

[Supply System of Fine-Water Droplets]

Fig.1 Experimental apparatus.

Fig.2 Characteristic of AFAD for fine-water droplet.

Fig.3 Characteristic of AFAL for fine-water droplet.
Fig.4 Characteristic of MGS and MGL for fine-water droplet.

近い状態となる様に設計した。また今回、流量レンジの異なる2種類の二流体噴射弁（AFAS、AFAL）と、霧化能力の異なる2種類の超音波霧化器(MGS、MGL)を採用した。二流体噴射弁へ供給する空気の質量流量 Qc [ml/min] と水の質量流量 Qw [ml/min] は、それぞれマスフローコントローラーで制御し、超音波霧化器との組み合わせにより微細水粒子体積流量および粒径制御を行った。これにより、微細水粒子の平均粒径を一定にして体積流量を可変したり、また体積流量を一定にして平均粒径を可変したりする事ができ、粒子特性に対する対向流拡散火炎の消炎現象を定量化する事が出来る。

また、二流体噴射弁（AFAS、AFAL）および超音波雾化器（MGS、MGL）の微細水粒子の粒子特性は、PDPAにより計測され、それぞれザウター平均粒径および体積流量の特性をFig.2、Fig.3、Fig.4に示す。

### 3. 実験結果および考察

### 3.1 対向流拡散火炎の微細水粒子混入時の直接写真

Fig.5 に対向流拡散火炎の直接写真を示す。Fig.5（a）、（b）と（c）、（d）は、各々、以下に定義するよどみ速度勾配 2V/R（V：空気流速 [m/s]、R：多孔質円筒バーナ半径 [m]）と無次元吹き出し速度 -fw の条件を合わせており、微細水粒子がない場合（Fig.5：（a）、（c））と、ある場合（Fig.5：（b）、（d））の条件のもと撮影を行った。微細水粒子添加時の体積流量は、5.3 [ml/min] 一定とした。

ここで 2V/R [s⁻¹] は、二次元速度ポテンシャル流れを仮定し、静止している円柱を過ぎる一様流から導く事が出来る。速度 V [m/s] の一様流の中に、半径 R [m] の円筒が置かれた場合、速度ポテンシャル ∇φは、

$$ \phi=V\left[r+\frac{R^{2}}{r}\right]\cos\theta $$

で表す事が出来る $ ^{1)} $ 。ここで r [m] は、円筒中心からの距離とする。また、よどみ点から任意の角度 $ \theta $ 位置における円筒表面に沿う速度成分 $ V_{\theta} $ [m/s] は、

$$ v_{_{\theta}}=\frac{1}{r}\frac{\partial\phi}{\partial\theta}=-V\left[1+\frac{R^{2}}{r^{2}}\right]\sin\theta $$

であるので、円筒表面上である r=R では、

$$ v_{_\theta}=-2V\sin\theta $$

となる。また、円筒表面上に沿うようにθ方向と同じ方向にS座標をとると、

$$ S=R\theta $$

と表す事ができ、また角度θが十分に小さい場合は、 $ \sin\theta=\theta $ の関係が成り立つため、円筒表面に沿う速度成分 $ V_{\theta} $ [m/s] は、

$$ v_{_{\theta}}=-2V\theta=-2VS/R $$

と簡易化される。すなわち円筒前方のよどみ領域における円周方向のよどみ速度勾配 $ [s^{-1}] $ は、

$$ \frac{dv_{\theta}}{dS}=\frac{\partial v_{\theta}}{\partial\theta}\frac{\partial\theta}{\partial S}=-2V/R $$

と表され、円筒前方よどみ領域に火炎が形成される場合は、一様流により得られる平均の火炎伸張率として見積もられ、火炎構造や火炎の振る舞いを支配する重要なパラメータ一である。

更に、-fwは、境界層の無次元吹き出し速度を式(7)で定義され、

$$ -f_{_{w}}=v_{_{w}}/V\sqrt{R_{e}/2} $$

となる。ここで Re は空気流の Reynolds 数、 $ v_{w} $ は燃料噴射速度 [m/s] を示す。

Fig.5(a)へ示す様に、ある $ v_{w} $ に対し、空気流速Vを増加する場合、すなわちよどみ速度勾配2V/Rが高い場合、火炎の厚さは徐々に減少し、層流二次元火炎が前方よどみ領域からある離れた位置で安定化される。火炎は非常に安定な二次元層流拡散火炎となる。 $ v_{w} $ が減少、もしくはVが増加した場合、火炎は円筒表面に近づき、最終的にはよどみ領域で消炎して、突然wake flame(後流火炎)へ移行する。一方、 $ v_{w} $ を増加しVを減少すれば、円筒からの火炎距離は徐々に増加する。しかしながら、 $ v_{w}/V $ が大きく増加した場合、大きなスケールの乱れが火炎面に現れ、層流二次元火炎を形成する事はできない。更に、Vが非常に大きい場合は、 $ v_{w} $ を十分に大きくしても火炎は消炎する。従って、消炎する限界のよどみ速度勾配(2V/R) $ c[s^{-1}] $ が存在する事になり、その値を越えて火炎を安定化する事はできない。

逆に、Fig.5(c)は、空気流速Vが低く、燃料噴射速度
$ v_{w} $ が大きい場合、すなわちよどみ速度勾配2V/Rが小さい場合、火炎は燃料側に黄色輝炎、酸化剤側に青炎が発生し、火炎の厚さは数mmとなり厚くなる。空気流速を限界まで低くすると黄色輝炎と青炎は完全に別れ、火炎を拡大し観察すると狭いダークゾーンが火炎の中間に現れる。

更に、Fig.5 (a) と (b) の微細水粒子添加時の火炎を比較すれば、微細水粒子を添加した場合、青炎である火炎帯の厚さは、無添加時の火炎帯厚さに比べ約35%薄くなり、発光も弱くなっている。これは、よどみ速度勾配の低い条件であるFig.5 (c) と (d) でも同様な傾向がみられ、火炎の酸化剤側にある青炎部および燃料側の輝炎が微細水粒子を添加する事で、約20%程度火炎帯が薄くなる事が観測された。

### 3.2 微細水粒子非混入時の対向流拡散火炎の安定限界

Fig.6は、-fwに対するよどみ速度勾配2V/Rを示している。また図中には、黄色輝炎が消滅する限界も示している。無次元吹き出し速度-fwの小さい領域では、-fwの増加と共に、2V/Rの値は急激に増加するが、-fwが十分大きくなると2V/Rは一定の値となり、-fwには依存しない消炎限界よどみ速度勾配（2V/R）cが存在する事が分る。-fwが小さい時の消炎は、多孔質円筒バーナー表面近傍での熱的消炎によるものであり、また-fwが大きい場合に観察される、一定の消炎限界よどみ速度勾配（2V/R）cは、化学反应速度の有限性のために生じるもので前者とは区別される。すなわち、後者の消炎は、Damköhler数によって支配されており、この値が臨界値に近づくと火炎が消炎する事が知られている。消炎限界よどみ速度勾配（2V/R）cは、時間の逆数の次元を持ち、流体力学的な特性時間の逆数に比例する。消炎時のDamköhler数が一定であると考えると、流体力学的特性時間は、化学反应の特性時間に比例する。従って、よど

Fig.5 Direct photographs of counterflow diffusion flame without water mist (a, c) and with water mist (b, d)

み速度勾配は燃料と酸化剤の総括的な反応速度の目安と考えられる。

また黄色輝炎消炎限界も同じ傾向を示しており、-fwが比較的小さい場合には、黄色輝炎消炎限界における2V/Rの値は小さい。しかし-fwの増加とともに2V/Rは急激に増加し、十分に大きな-fwでは一定となる。この事は黄色輝炎の消炎に関しても限界のよどみ速度勾配が存在する事を示している。前者の黄色輝炎発生限界においては、燃料流速が遅く火炎が凹筒バーナー表面に近づくため、熱損失の影響が大きくなり、その結果、燃料の熱分解速度が低下するため輝炎が消炎する。また後者の黄色輝炎消炎限界では、火炎帯における固体炭素生成にかかわる総括的な化学反応速度よりも、よどみ速度勾配の方が大きくなると考えられる。

### 3.3 微細水粒子による消炎限界

小さな水粒子は、大きな水粒子より容易に蒸発する事が出来るため、消炎におよぼす効果は、微細な水粒子の方が大きくなる事は予想出来る。また、水の持つ熱容量も消炎に寄与する重要なパラメータであり、空気流中に添加する微細水粒子の質量も重要である。本研究において、微細水粒子の蒸発過程は、対流熱伝達が支配的となり、水粒子の蒸発面積が対向流拡散火炎の消炎において重要である事は予想される。そのため微細水粒子の蒸発面積と質量流量の比である面積パラメータ∑s $$ m^{2}/kg $$ により、対向流拡散火炎の消炎を評価する事が出来る。

式(8)に面積パラメータ∑s $$ m^{2}/kg $$ を示す。

$$ \varSigma_{S}=\frac{\dot{A}_{w}}{\dot{m}_{a}+\dot{m}_{w}}=\frac{6Y_{0}}{\rho_{w}}\frac{\sum\limits_{i=1}^{n}c_{i}D_{i}^{2}}{\sum\limits_{i=1}^{n}c_{i}D_{i}^{3}}=\frac{6}{\rho_{w}}\frac{Y_{0}}{D_{s}} $$

$ Y_{0} $ は水の質量分率、 $ Y_{0}=\dot{m}_{w}/(\dot{m}_{a}+\dot{m}_{w}) $ で定義され、 $ \dot{m}_{a} $ と $ \dot{m}_{w} $ は、それぞれ空気と水の質量流量 [kg/s] である。

Fig.6 Flame stability diagram for methane counterflow diffusion flame.