簡易圧縮性
k-ε モデルを用いた単体建物周辺における
高温排気ガスの拡散予測に関する研究
STUDY OF DISPERSION PREDICTION OF HIGH TEMPERATURE EXHAUST GAS AROUND ISOLATED BUILDING USING WEAKLY COMPRESSIBLE k-ε MODEL
林 超１） 大岡 龍三２） 菊本 英紀３） 佐藤 大樹４） Chao LIN1), Ryozo OOKA2) , Hideki KIKUMOTO3) and Taiki SATO４)
ABSTRACT
Dispersion of high temperature exhaust gas around an isolated cubic building was predicted by using the weakly compressible k-ε model. The density of the exhaust gas 𝜌𝜌� was set as
0.3 and 0.1 times of that of environmental air 𝜌𝜌�. In the condition of 𝜌𝜌� = 0.3𝜌𝜌�, the results
of concentration were almost the same with the incompressible model and larger than the experimental data. However, in the condition of 𝜌𝜌� = 0.1𝜌𝜌� which lacks experimental data,
the weakly compressible k-ε model showed larger concentration on leeward building wall than incompressible model. The difference between the two models are discussed in the aspect of buoyant acceleration and production of turbulence kinetic energy and its dissipation. Key Words: CFD, Dispersion, Compressibility, Low Mach Approximation
１．はじめに
近年、計算機性能の向上に伴い、数値流体力学(Computational Fluid Dynamics: CFD)に基づく環境シミュレーシ
ョンは、実際の市街地を対象とした都市風環境の予測や評価に利用可能な段階に到達しつつある 1)。ボイラーや発
電機の運転時には、建物周辺で極めて高温な排気ガスが放出される。人々の安全性や建物への影響を解析し適切
な対応策を立案するためにも、CFD を用いた拡散予測が有効であると考える。
CFD を用いた高温排気ガスの拡散予測に関しても先行研究として、例えば、有波ら(2013)2)は可搬型ガスタービン
発電機の高温排熱の拡散性状を実測および数値解析をし、外気温の相違によって高温排熱の拡散性状にほとんど
変化がないことを報告した。Tong and Zhang(2015)3)は都市環境の範囲でディーゼルバックアップ発電機の排出口の
位置が高温ガスの拡散性状に大きく影響を与えることを指摘した。Tominaga and Stathopoulos(2018)4)は地表面付近
１） 東京大学大学院工学系研究科 大学院生 （〒153-8505 目黒区駒場４－６－１） ２） 東京大学生産技術研究所 教授 （〒153-8505 目黒区駒場４－６－１）
３） 東京大学生産技術研究所 講師 （〒153-8505 目黒区駒場４－６－１）
における異なる浮力を伴うガスの拡散性状を解析し、乱流エネルギー消散率の輸送方程式における浮力モデルが結 果にほとんど影響しないことを報告した。 従来の CFD 解析では、以上のような浮力を伴う拡散現象の予測を行う場合でも、非圧縮性の仮定およびブジネス ク近似によって密度の変化は直接扱わず、運動方程式での浮力効果のみを考慮することが多い。しかし、流体密度 の変化が無視できない高温ガスの拡散においては、従来手法の解析精度が十分に検討されているとは言えない。そ こで本研究では、流体の圧縮性を考慮するk-ε モデルを導入し、非圧縮性の予測モデルとの比較を行う。まず、予測 モデルの概要を説明し、既往の浮力を考慮した拡散実験 5)の再現解析を実施する。さらに、実験値がない高浮力条 件に対しても解析を実施し、平均濃度分布の結果について予測モデルと浮力効果の相違に関する考察を行う。 2．解析モデル 本研究では、定常状態の平均速度場および濃度場を対象に Realizable k-ε (RKE)6)モデルを基に高温排気ガスの 拡散予測を行う。ガスの温度だけでなく、それに含まれる汚染物質の濃度も予測対象とする。基礎方程式では、連続 の式および平均流速・温度・濃度、乱流エネルギー、乱流エネルギー消散率の輸送方程式を解く。 流体の圧縮性を厳密に考慮する場合には、気体の状態方程式によって、流体密度を温度および圧力の関数とし て求める必要がある。その際には、圧力波の伝搬過程も捕捉するため、時間離散化幅Δt を極めて小さくする必要が 生じ、計算実行上の困難が伴う 7)。しかし、本研究で対象とする問題は、流体速度が圧力波の伝搬速度よりも十分小 さい。したがって、流体密度𝜌𝜌を求める際に、圧力を基準大気圧𝑃𝑃�で一定とし、温度 T のみの影響を考慮する低マッ ハ数近似（式(1)）の適用が妥当であると考える。 𝑃𝑃�� 𝜌𝜌RT (1) ここで、R は気体定数である。本研究では、上記の近似を用いる手法を簡易圧縮性 k-ε モデルと呼ぶ。簡易圧縮性 モデルでは浮力𝐟𝐟�および体積膨張率𝛽𝛽は、基準流体密度𝜌𝜌�、重力加速度𝐠𝐠を用いて次のように表される。 𝐟𝐟� � �𝜌𝜌 � 𝜌𝜌��𝐠𝐠𝐠 𝛽𝛽 � �1𝜌𝜌𝜕𝜕𝜌𝜌𝜕𝜕𝜕𝜕 �1𝜕𝜕 (2) 一方、従来から用いられる非圧縮性の仮定およびブジネスク近似を用いたモデルにおいては、浮力および体積膨 張率𝛽𝛽は、基準温度𝜕𝜕�を用いて次のように与えられる。 𝐟𝐟� �𝜌𝜌��𝜕𝜕�𝜕𝜕� 𝜕𝜕�𝒈𝒈 � � 𝜌𝜌�𝛽𝛽��𝜕𝜕�� 𝜕𝜕�𝐠𝐠𝐠 𝛽𝛽�� 1 𝜕𝜕� (3) 3．解析概要 3.1． 解析対象 本研究では、富永ら(1992)5)の風洞実験を再現したCFD 解析を行なう。図 1 に解析領域を示す。一辺 Hb = 0.2 m の立方体建物が領域内に配置され、座標原点は地表面における建物背面の中心点である。この建物の後流域内の 床面（立方体建物風下端から 0.5Hb風下）にガスの排出口（排出口形状は一辺 0.025Hbの正方形）が設置されてい る。流入部での軒高平均風速UHは 0.4 m/s であり、環境空気と異なる密度のガスが、上記の排出口から 0.5UHの鉛 直方向速度をもって定常的に放出される。 図1 解析対象 図 2 流入面における無次元速度および乱流強度のプロファイル 0 1 2 3 4 5 0 0.2 0.4 z/H IU 0 1 2 3 4 5 0 1 2 z/H U/UH 1/4 power law
3.2． 解析条件 解析にはSTAR-CCM+ V12.02.010 を用いた。解析条件を表 1 に示す。流入境界は図 2 に示しているように、風洞 実験でデータが得られているものは実験データをべき関数などで近似し与えた。建物前方の地面の境界条件には𝑧𝑧� 型の対数則を用い、流入プロファイルの形状を保持するため、試行錯誤的に𝑧𝑧�� ���������とした。その他の量に 関しては、日本建築学会の流体数値解析ガイドライン1)を参考に推定し与えた。 実験 5)では、密度の異なるガスを用いて浮力を生じさせるが、本研究では等価な浮力を生み出す温度を排出ガス に与え解析する。ただし、排出ガスにはトレーサーを模擬した濃度を与え、濃度の結果を実験と比較する。乱流プラン トル数および乱流シュミット数は、実験が物質拡散に基づくためともに0.7 とした。 体積濃度に基づく基準濃度を𝐶𝐶�� 𝑞𝑞��𝑈𝑈�𝐻𝐻���で定義する（𝑞𝑞はガス発生流量[m3/s]）。本解析条件では𝐶𝐶� = 312.5 ppm である。結果は𝐶𝐶�で無次元した平均濃度𝐶𝐶∗で示す。 表 1 解析条件 解析モデル 定常RANS, Realizable k-ε 2 層モデル 解析領域 16 Hb (x) ×9 Hb (y) ×6 Hb (z) 解析格子 1,260,622 (正六面体格子) アルゴリズム SIMPLE 移流項の差分スキーム 2 次精度風上差分 拡散項の差分スキーム 2 次精度中心差分 流入条件 実験値から近似 地面 (x=-6 Hb ~-2 Hb) z0型境界(z0=0.00183 m) 上空, 側面, 建物表面, 地面(x=-2 Hb~ 10Hb) 滑面型壁関数, 断熱 3.3． 排出ガスの密度条件と解析ケース 実スケールと風洞実験スケールでの現象の相似性は、現象を支配する特徴量から構成される無次元数を用いて 議論されることが多い。本研究では、式(4)に定義する浮力無次元数によってその相似性を考慮する。 𝜌𝜌�𝑈𝑈�� �𝜌𝜌�� 𝜌𝜌��𝑔𝑔𝐻𝐻� (4) ここで𝑔𝑔は重力加速度の大きさである。この無次元数は、圧縮性流体の運動方程式における慣性力と浮力の比の 大きさに対応している。ただし、例えば富永ら(1992)5)の論文で用いている密度フルード数とは異なり、慣性力を表す 項に基準密度ではなく排出ガスの密度を使用していることに注意されたい。風洞実験と実現象との相似性（浮力無次 元数の一致）を仮定すると、実験条件は浮力無次元数を介して実スケールの対応するパラメータの値に変換できる。 排出ガスの密度は、既出の風洞実験では𝜌𝜌� = 0.3𝜌𝜌�であった。しかし、東京都の年平均風速� � �����・年平均気 温𝑇𝑇�� �����������・建物高さ 40 m を想定する実スケールの条件とし、標準大気圧下での空気の気温と密度の関 係も用いると、この風洞実験は約 200℃のガスが排出される状況を再現しているものと解釈できる。しかし、非常用発 電機の排出ガスの温度は 600℃ (873 K)を超える場合も多く、その再現としては浮力の効果が不足している。そこで 本研究では、𝜌𝜌� = 0.1𝜌𝜌�の排出ガス密度条件を設定し、CFD による数値解析を実施した。これは前述の各種条件で 相似性を考慮すると、実スケールで約 740℃の排出ガス温度の状況に相当する。なお、風洞実験において低密度ガ スを生成するのに多く使用されるヘリウムでも空気に対して約 0.14 の比重をもつ。したがって、上記の密度条件を実 験的に実行するのは難しいため、本条件はあくまで数値実験として設定した。 表2 に示す 4 ケースの解析を行った。Case1 と Case2 はそれぞれ𝜌𝜌� = 0.3𝜌𝜌�と𝜌𝜌� = 0.1𝜌𝜌�の排出ガス密度条件で非 圧縮性と簡易圧縮性モデルを用いるケースである。環境温度は 300 K に設定した。ガス出口の温度については、密 度条件からそれぞれブジネスク近似および理想気体状態式から計算して与えた。 4．解析結果 4.1． 平均濃度分布の比較 表 2 解析ケース ケース名 ρs/ρ0 浮力無次元数 圧縮条件 ガス出口の 設定温度 (K) 環境 温度 (K) Case1-1 0.3 -0.116 非圧縮性 512 300 Case1-2 簡易圧縮性 1000 Case2-1 0.1 -0.091 非圧縮性 570 Case2-2 簡易圧縮性 3000
まずは既往実験値がある Case1 の結果を示す。図 3 は、建物風下側での平均濃度の鉛直分布と水平分布 である。実験値およびTominaga ら4)のCFD 解析値（非圧縮性・ブジネスク近似・RKE6）を使用）をともに載せる。 平均濃度は、非圧縮性・圧縮性いずれにおいても実験より大きく評価する傾向が見られる。その原因は k-ε モ デルでは建物周辺におけるスパン方向の運動量およびスカラーの拡散が小さく評価されていることによるものと 考えられる 8) 。圧縮性と非圧縮性のモデリングによる明確な差異はこれらの位置でのプロファイルでは見られな かった。また、解析条件に若干の違いはあるものの既往のCFD 解析値4)とも本解析結果は概ね一致した。 図3 無次元化平均濃度のプロファイル 上段：鉛直断面 (y/Hb = 0), 下段：水平断面 (z/Hb = 0.5) 図4 無次元化平均濃度のプロファイル 上段：鉛直断面 (y/Hb=0), 下段：水平断面 (z/Hb=0.5) 次に高浮力のCase2 に着目する。図 4 は、𝜌𝜌� = 0.1𝜌𝜌�の建物風下側での濃度の鉛直および水平プロファイルであ る。比較のため𝜌𝜌� = 0.3𝜌𝜌�の解析結果をともに載せる。また図5 は、解析領域のスパン方向中心断面内における平均 濃度の分布である。Case2-1(𝜌𝜌�/𝜌𝜌�� ���)では、Case1-1(𝜌𝜌�/𝜌𝜌�� ���)と比較して濃度分布にほぼ差が現れなかった。 一方、Case2-2(𝜌𝜌�/𝜌𝜌�� ���)では、Case1-2(𝜌𝜌�/𝜌𝜌�� ���)より上方に濃度が輸送される傾向が見られる。その反面、ス パン方向(y 方向)への濃度の広がりはやや狭まっている。非圧縮性と簡易圧縮性モデルとの差は、排出ガスがより低 密度の条件で顕著になった。図 6 は、建物風下側鉛直面における平均濃度分布である。それぞれのケースにお いて面内で最大の濃度を示す位置は、Case2-1 と比べ Case2-2 でより高い位置に生じている。また、その最大値 もCase2-1 よりも Case2-2 において大きな値を示すような結果となった。 Case2 の条件においては、実験値が存在しないため CFD 解析値の精度評価を行うことはできない。しかし、 簡易圧縮性モデルの方がより精緻に物理現象をモデリングしていることを考慮すると、従来の非圧縮性モデル による予測では、建物表面に到達する排出ガス濃度（あるいは温度）を過小評価する可能性が示された。
(a) Case2-1：非圧縮性モデル (b) Case2-2：簡易圧縮性モデル
(a) Case2-1 (b) Case2-2 図6 建物風下側鉛直面 (x/Hb =0)における無次元濃度の分布 図7 排出ガスの密度と浮力による加速度との関係 4.2． 予測モデルと浮力による加速度の傾向 前節までの結果を浮力による加速度の面から考察を行う。改めて浮力は次式で表現される。 𝐟𝐟�� �𝜌𝜌 � 𝜌𝜌��𝐠𝐠 (5) 非圧縮性の仮定およびブジネスク近似を用いたモデルでは、式(2)を近似した表現で置き換えるが、ガス排出 口においては等価な浮力を与える温度を設定するため、少なくとも排出口近傍では𝐟𝐟�の大きさは、2 つの予測モ デルにおいて大きくは異ならないと考えられる。 しかし、流体運動の観点からは浮力そのものよりもそれをその場所の密度で割った浮力による加速度𝐚𝐚�への 方 が重 要 である。𝐚𝐚�は、圧 縮 性 を考 慮 する場 合 には、式(5)の両辺をその場の流体密度で割ることによって式 (6)のように表される。一方、非圧縮性を仮定する場合には、密度変化による浮力は考慮するが、密度の変化自 体は考慮しない。つまり非圧縮性のモデルにおいては、式(5)の両辺を基準密度で割ることで式(7)のような加速 度を評価していることを相当する。 𝐚𝐚��𝐟𝐟𝜌𝜌 � �1 �� 𝜌𝜌𝜌𝜌 � 𝐠𝐠� (6) 𝐚𝐚��𝜌𝜌𝐟𝐟� �� � 𝜌𝜌 𝜌𝜌�� 1� 𝐠𝐠 (7) この 2 つの予測モデルによる加速度の違いをグラフにしたのが図7 である。𝑔𝑔 = 9.81 m/s2とした。排出ガスの 密度が小さくなるほどに両モデルにおける加速度の乖離が拡がっていき、流体の圧縮性を考慮したモデルほど より強い浮力が働く傾向が強まる（密度が大きいガスの場合は逆だが、その差はそれほど大きくない）。そして、 この傾向は図6 で示すように Case2-2 における最大の濃度がより高い位置に生じた結果とも整合している。 4.3 予測モデルと乱流拡散 𝐺𝐺�� 𝛽𝛽𝑃𝑃𝑃𝑃𝜇𝜇� ��∇𝑇𝑇� ∙ 𝐠𝐠� (8) 𝐶𝐶��� ���� |𝑣𝑣�| |𝑢𝑢�| (9) 𝛽𝛽 � � 1 𝜌𝜌 𝜕𝜕𝜌𝜌 𝜕𝜕𝑇𝑇 � 1 𝑇𝑇 (2) ��𝛽𝛽 � 𝛽𝛽�� 1 𝑇𝑇� (3) 本節では乱流拡散の面から考察を行う。Realizable k-ε 2 層モデルでは、乱流エネルギーk の輸送方程式にお ける浮力生成項𝐺𝐺�は式(8)で示される。2 節で述べたように、体積膨張率𝛽𝛽は簡易圧縮性と非圧縮性モデルに おいて、それぞれ平均温度𝑇𝑇と基準温度𝑇𝑇�を用いて評価される。 乱流エネルギー消散率ε の輸送方程式では、浮力生成項は k の浮力生成項𝐺𝐺�に係数𝐶𝐶��(式(9))を乗じて表 される。𝐶𝐶��においては、𝑣𝑣�は重力加速度と同じ方向の速度、𝑢𝑢�は重力加速度と垂直な方向の速度である。つま り、流れ場が重力加速度の方向と完全に一致する場合は1 で、垂直な場合は 0 になる9)。 図8 は、建物風下側水平面における渦粘性係数の分布である。建物近傍において、非圧縮性モデルより、簡 易圧縮性モデルの値が小さい傾向が見られる。乱流エネルギーk の浮力生成項𝐺𝐺�は体積膨張率𝛽𝛽に比例する。 𝛽𝛽は、非圧縮性モデルでは常に環境基準温度𝑇𝑇�� �����の逆数で、簡易圧縮性モデルではその場所の平均温 度の逆数であるため排出口近傍の高温域ではより小さい値となる。したがって、簡易圧縮性モデルにおいて 乱流エネルギーk の浮力生成項𝐺𝐺�は非圧縮性モデルより小さくなる傾向にある。 一方、前節で述べた浮力加速度𝐚𝐚�により、簡易圧縮性モデルの方がより大きなｚ方向速度が生成される。した がって、式(9)より𝐶𝐶��は全体的により大きい値になると考えられる。その結果、乱流エネルギー消散率ε の浮力生 ‐100 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 ab (m /s 2) ρ/ρ0 ⾮圧縮性 圧縮性
成項𝐶𝐶��𝐺𝐺�も非圧縮性より簡易圧縮性モデルの方がより大きくなると考えられる。 k-ε モデルでは渦粘度𝜈𝜈�がk の二乗に比例し、 ε に反比例する(𝜈𝜈�� 𝐶𝐶�� � �)ため、渦粘度𝜈𝜈�は非圧 縮性モデルより簡易圧縮性モデルで小さくなり、 乱流拡散が抑制されるものと予想される。そして、 この傾向は図 5 や図 6 で示されるように、建物近 傍での濃度がCase2-1 より Case2-2 で高い値を示 す結果とも整合している。 5． まとめ (1) 浮力を伴う汚染物質の拡散予測に関して、簡易圧縮性 k-ε モデルを用い、既往の風洞実験の再現解析お よび高浮力の予測解析を行なった。 (2) 排出ガスの密度が環境空気の 0.3 倍の条件においては、非圧縮性モデルと概ね同様な結果を与え、実験 値より大きな濃度予測値を与えた。ただし、排出口近傍の高 濃度域においては、予測モデルによる差異も 見られ、簡易圧縮性モデルにおいて非圧縮性モデルよりも高い濃度が建物に到達する結果を与えた。 (3) 排出ガスの密度が環境空気の 0.1 倍の条件においては、簡易圧縮性モデルの方が浮力による加速度がよ り大きく働き、排出ガスの濃度がより上方に輸送されるような分布となった。 (4) 予測モデルと乱流拡散の傾向について考察を行った。非圧縮性モデルに比べて簡易圧縮性モデルにお いては、乱流エネルギーk の浮力生成項𝐺𝐺�が小さくなると同時に乱流エネルギー消散率 ε の浮力生成項 𝐶𝐶��𝐺𝐺�が大きな値になる。その結果、渦粘度𝜈𝜈�が小さくなり乱流拡散が抑制されることを示した。 参考文献
1) Y. Tominaga et al.: AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings, J. Wind Eng. Ind. Aerodyn., vol. 96, no.10–11, pp.1749–1761, 2008.
2) 有波裕貴, 赤林伸一, 石川真也, 他:可搬型ガスタービン発電機からの高温排熱の拡散性状に関する研 究, 日本ガスタービン学会誌, vol. 41, no. 3, pp. 246–253, 2013.
3) Z. Tong and K. M. Zhang: The near-source impacts of diesel backup generators in urban environments, Atmos. Environ., vol. 109, pp. 262–271, 2015.
4) Y. Tominaga and T. Stathopoulos: CFD simulations of near-field pollutant dispersion with different plume buoyancies, Build. Environ., vol. 131, no. October 2017, pp. 128–139, 2018.
5) 富永禎秀, 村上周三, 持田灯, 他:浮力のあるガスが排出された場合の建物周辺の濃度変動、乱流拡散 構造に関する風洞実験, 風工学シンポジウム論文集, vol. 12, pp. 119–124, 1992.
6) T.H. Shih et al.: A NEW k-ε EDDY VISCOSITY MODEL FOR HIGH :REYNOLDS NUMBER TURBULENT FLOWS, Computers Fluids, vol. 24, no. 3, pp. 227–238, 1995.
7) 白石靖幸、加藤信介、石田義洋:低マッハ数近似式との比較による Boussinesq 近似式の予測精度の検 討,日本建築学会環境系論文集, vol. 577, pp. 13–18, 2004.
8) B. Blocken et al.:Numerical evaluation of pollutant dispersion in the built environment: Comparisons between models and experiments, J. Wind Eng. Ind. Aerodyn., vol. 96, no. 10–11, pp. 1817–1831, 2008. 9) R. A. W. M. Henkes et al.,:Naturalconvection flow in a square cavity calculated with lowReynolds
-number turbulence models, Int. J. Heat Mass Transf., vol. 34, no. 2, pp. 377–388, 1991. 謝辞
本研究の実施にあたり、解析対象とした風洞実験を実施した渋谷亜紀子氏（当時東京大学大学院大学院生）より、 修士論文の提供を受けた。記してここに謝意を表する。
(a) Case2-1：非圧縮性モデル (b) Case2-2：簡易圧縮性モデル 図8 建物風下側水平面 (z/Hb =0.25)における渦粘性係数𝜈𝜈�の分 布
