三菱重工技報 Vol.49 No.1 (2012) 新製品・新技術特集

技 術 論 文75

複雑地形へ適用可能な非定常乱流拡散
数値シミュレーション技術の開発－風洞実験代替に向けて－Development of Unsteady Numerical Simulation Modelfor Gas Diffusion Over Complicated Terrains by Using Large-Eddy Simulation-Numerical Model for Replacement of Wind Tunnel Experiments-

岡 林 一 木 *1

永 山 圭 憲 *2

Kazuki Okabayashi

Yoshinori Nagayama

原智 宏 *3

堀英 一 *4

Tomohiro Hara

Eiichi Hori

大 屋 裕 二 *5

内 田 孝 紀 *6

Yuji Ohya

Takanori Uchida

今般，原子力学会標準委員会では，原子力関連施設の安全審査における周辺被ばく評価に関して，従来行われてきた風洞実験の代わりに数値シミュレーションを用いるための基準作りが始められている．これに対応するために，乱流モデルの１つであるラージ・エディ・シミュレーション（LES）を用いた非定常乱流計算を複雑地形上の拡散場へ適用して，これまでの風洞実験結果との比較検証を行った．その結果，実用的な計算時間内で，精度の良い予測結果が得られ，本解析技術が風洞実験に代替可能な非定常乱流拡散計算法となることが確認された．

|1. はじめに近年の計算機の性能向上及び高速計算技術の進展により，これまでかなりの時間を要していた大規模な非定常乱流計算が，実用的な計算時間で可能となりつつある．こうした中，現在，原子力学会では，周辺被ばく評価のために実施している風洞実験(1)を数値シミュレーションで代用するための標準制定に向けた取組みを開始している(2)．
現在，乱流モデルとして計算時間が比較的少ないレイノルズ平均応力モデル（RANS）が広く一般に利用されており，既に風洞実験を代用する目的で，比較的平坦な地形や単純な建物を対象にした RANS による検証もなされつつある(3)．しかし，複雑地形を対象とした場合には検証例は少なく，計算精度の検証が必要であり，複雑な流れ場及び拡散場への適用には精度的な限界もある(4)(5)．そこで本研究では，乱流モデルの１つであるラージ・エディ・シミュレーション（LES）を複雑な地形上の非定常拡散場へ適用することを試みた． LES は計算格子よりも小さなスケールの乱れをモデル化して計算格子よりも大きなスケールの乱れを直接解くものであり，RANS よりも計算時間を要する．これまで，LES については多くの研究例が見られるが，大気拡散分野ではほとんどが単純な形状を対象とした基礎研究に限られており(6)，複雑地形を対象にした実用的な大規模計算例は見られない．そこで，複雑地形を対象に LES による非定常拡散計算を実施し，別途当社ほかにて実施された既存の風洞実験結果と比較検討を行うことで計算手法の有効性を検証した．

*1 技術統括本部長崎研究所主席研究員 博（工）*2 技術統括本部長崎研究所*3 技術統括本部長崎研究所 博（工）*4 原子力事業本部原子力プラント技術総括部原子炉安全技術部*5 九州大学応用力学研究所 新エネルギー力学部門 風工学分野 教授 工博*6 九州大学応用力学研究所 新エネルギー力学部門 風工学分野 准教授 博（工）

三菱重工技報 Vol.49 No.1 (2012)76|2. 大気拡散数値モデルの概要
2.1 気流計算
流れの支配方程式は，式(1)及び式(2)に示すとおり，非圧縮性流体の連続の式とナビエ・ストークス方程式である．これらの解析にあたっては九州大学応用力学研究所で開発された気流解
析モデル RIAM-COMPACT(7)を用いた．
u i
0
xi

(1)

 ij
u i
u
p
1  2uiuj i t
x j
x i Re x j x j x j

(2)131i j0i j ij  u iu j  u k u k  ij  2 SGS S ij ，  ij   SGS  C s f s   S， S  1  u i  u j ij
2  x j xi 2S  2 S ij S ij

 ， f  1  exp z 25 ，   h h h 
1213



xsyz

ただし，A ：諸量 A に対するフィルタ操作(ここでは，計算格子幅による平均値)
xi ：i 方向座標(x =x：主流方向，x =y：水平方向，x =z：鉛直方向)
12

ui ：速度の i 方向成分ui ：速度の i 方向変動成分

3p ：圧力Re：レイノルズ数

SGS：サブグリッドスケール(計算格子幅以下のスケール)の拡散係数
Cs：標準スマゴリンスキーモデル係数(=0.1)
fｓ：減衰関数z+：壁座標hx，hy，hz：各方向成分のフィルタ幅(ここでは，計算格子幅)
ここで，サブグリッドスケールの渦による応力を表すモデルとして，標準スマゴリンスキーモデル
3を用いた．なお，上式中で重複する添え字には総和規約が適用される．すなわち，a b は  a j b jj jj 1を意味する．境界条件としては，側方境界面と上部境界面に滑り条件，流出境界面は対流型流出条件とし，地面は Non-Slip 条件とした．

2.2 拡散計算
大気拡散計算モデルとして，以下の拡散方程式を用いた．
C
C
C
C   C    C    C  Kuv
w Kq
 K
t
x
y
z x  x  y  y  z  z （3）1K SGSRe Sc Sc turb

ただし，
C：濃度

K：拡散係数

Scturb：乱流シュミット数

Sc：シュミット数

q：煙源におけるガスの放出強度

ここで，乱流シュミット数 Scturb は計算結果に大きな影響を与えるものではないため，実験結果との比較も踏まえ，極力小さめに設定した．なお，シュミット数 Sc は 0.71 とした．
拡散計算モデルは，RIAM-COMPACT と同様に，対流項及び拡散項の差分スキームに，それぞれ３次精度風上差分と２次精度中心差分を用いた．時間進行法はオイラー陽解法である．境界条件は，流入境界面において濃度を０とし，他境界面及び地形境界面ではすべてノイマン条件とした．

三菱重工技報 Vol.49 No.1 (2012)77|3. 検証方法
EU 共同研究(SMEDIS project)(8)では，複数の研究機関が CFD モデル(STAR-CD，FLUENT，CFX ほか)による実験データの再現計算を実施し，その結果に基づいて数値モデルの精度検証ガイドライン(COST732)が作成された (9) ．これに従い，第２章で記述した本計算手法の性能をCOST732 に従って評価する．なお，COST732 による性能基準は，直方体建屋背後の風下地表面から放出されたガスの濃度について評価を行い，良好事例(Best Performance)からその適合性を定めたものである．検証は、COST732 に示される基本的な検証手法例に基づき、以下の３つの評価指標を用いて行う．
①FAC2（Fraction of Xp within a factor of two of Xo）：ファクター2
次式によって定義される FAC2 は，実験値 Xo と計算値 Xp の比が 0.5～2.0 に入る割合を示す．

FAC 2 nN 1
  Nin n i 1Xp
 2.0
1 for 0.5 N Xo
else 0

②FB（Fractional Mean Bias）：バイアス諸量 A の平均値を A で表すと，FB は次式で定義される．これは実験値 Xo と計算値 Xp の差（誤差）の平均値相当のものであり，全体的な偏移量を表すための指標である．

FB  2( X o  X p ) /( X o  X p )③NMSE（Normalized Mean Square Error）：正規化二乗誤差
NMSE は次式で定義され，誤差（データ）の広がりを表すものである．

NMSE  (( X o  X p ) 2 ) /( X o  X p )|4. 計算精度の検証
第２章で記述した大気拡散数値モデルの精度検証のため，第３章に示した COST732 に基づく評価基準(9)を用いて，過去に当社で実施した原子力関連施設の安全審査における周辺被ばく評
価のための風洞実験結果ほかとの比較検討を行った．

4.1 平地実験との比較
はじめに，平板上で風洞実験により得られた気流場及び所定の煙源高さに応じた濃度場が数値シミュレーションにより再現できるかを検討した．計算では風洞実験と同様に，風路前流側に地表面粗度としておかれた境界層制御用のスパイヤー（高さ 520ｍｍの二等辺三角形）を再現して，風下方向へ発達する平板上の乱流境界層の計算を実施した．併せて，その境界層中で異なる高さから放出されるガスの拡散計算も実施した．計算格子幅は，流れ方向に 13ｍｍ～26ｍｍ，流れと直角方向に一定の 6.5ｍｍ，高さ方向には 1.5ｍｍ～130ｍｍとして地表面付近が密となるように格子を寄せ，全格子数としては約 1900 万メッシュとした．計算は 16CPU(CentralProcessing Unit：中央処理装置)の並列機を用いて行い，後述する気流データベースも利用することで，50 時間で計算可能となった．図１に，煙源位置での平均風速及び乱れ強さ３成分の高さ方向分布の計算結果を，実験値と比較して示す．計算結果は 6.5 秒間の平均値であり，平均化時間によってその値がほぼ変化しないことを事前に確認している．図１より計算結果は実験値と非常によい一致を示していることがわかる．なお，計算結果は，風洞実験と同様に縮率 1/2000 としてスケール換算している．図２には，再現された気流場にて煙源高さを変化させ，最大地表濃度の風下方向分布（地表煙軸濃度分布）を実験値と比較した結果を示している．図中の濃度 C は正規化濃度 UC/Q（m-2）（U：一様流風速，Q：放出量）で整理しており，各煙源高さの地表煙軸濃度分布は実験値とよく整合していることがわかる．

三菱重工技報 Vol.49 No.1 (2012)78

さらに図３では，地上点源より放出された場合のガスの鉛直方向拡散幅（濃度分布の標準偏差に相当する長さ）の実験と計算の比較を行い，ほぼ一致していることを確認している．この拡散幅は，濃度分布を正規分布と仮定して，最小二乗法により求めたものである．これらの結果を，大気安定度ごとに拡散幅が定義されている Pasquill-Gifford 線図（1）とともに示すことで，再現目標とするほぼ中立状態の大気安定度Ｃ～Ｄを再現できていることが確認できる．

図１ 平板上での気流分布
平均風速及び乱れ強さ３成分の高さ方向分布に関して実験と計算との比較を示す．
（縮率 1/2 000 で実機スケール換算）

図２ 煙源高さ別の風下方向地表煙軸濃度分布

図３ 鉛直方向拡散幅

平 板 上 で の 各 煙 源 高 さ （ 0 ， 20 ， 40 ， 60 ， 80 ， 100 ，200m）における地表煙軸濃度分布の実験値との比較を示す．（縮率 1/2 000 で実機スケール換算）

平板上での鉛直方向拡散幅 σｚ の実験と計算との比較を示す．（縮率 1/2 000 で実機スケール換算）

なお，次節に示す複雑地形の計算では，図４に示すように，助走区間 A で得られた所定の風下断面での気流計算結果（気流データベース）を，計算領域 B の流入境界条件としても利用することで，計算領域から助走区間 A を削減し，計算時間の大幅な短縮を図っている．

図４ 計算領域の削減
助走区間 A の風下断面の気流データを，複雑地形の計算領域 B の流入条件として利用する．

三菱重工技報 Vol.49 No.1 (2012)

79

4.2 複雑地形実験との比較
次に，複雑地形上での計算手法の有効性を確認するため，図５に示す複雑地形を対象に計算と実験の比較を行った．放出源位置は，図中の A，B の２箇所とし，風向はⅠ～Ⅴの５風向とした．今回実施した計算ケースは，表１に示すとおり，放出源位置，風向，放出源高さを変化させた合計７ケースである(10)．地表煙軸濃度の計算結果を，風上側が海となる平坦な場合（風向Ⅰ）と複雑な地形となる場合（風向Ⅴ）の代表的な２ケースについて図６に示す．両者共に実験値と非常に良い一致を示している．なお，第３章で示した COST732 に基づき評価した３つの評価指標（FAC2，FB，NMSE）の結果を全７ケースについて表２に示す．表２に示すとおり，いずれのケースも，各評価指標の基準値を満足しており，十分な予測精度を有していることが確認できた．
表１ 計算条件一覧
ケース No.

風向

放出源位置

①

放出源高さ(m)

AI

②③144

A58B63④Ⅱ

A58⑤ⅢA144⑥Ⅳ

B63⑦ⅤB63

各計算ケースにおける放出源位置と風向，放出源高さを示す．

図５ 計算対象の複雑地形
計算対象の風向（Ⅰ～Ⅴ）及び２箇所の放出地点（A，B）を示す．なお，図中の線は対象とした地形の等高線を表す．

(a) 風向Ⅰ，放出源位置 A，放出源高さ 144m(b) 風向Ⅴ，放出源位置 B，放出源高さ 63m

図６ 風下方向地表濃度分布の比較
風下方向の地表煙軸濃度分布について実験と計算の比較結果を示す．（縮率 1/2 000 で実機スケール換算）

表２ 本計算手法の性能評価結果
ケース No.

FAC2
煙軸データ
全地表面データ(基準：0.89 以上)(基準：0.54 以上)

煙軸データ(基準：0.5 以下)

全地表面データ(基準：0.3 以下)

NMSE(全地表面データ)(基準：4 以下)

FB

①

0.95

0.82

0.210.210.33②1.000.71

0.32

0.11

1.19

③1.000.720.15

0.12

0.71

④0.890.65

0.120.07

0.52

⑤0.95

0.60

0.20

0.27

2.80

⑥1.00

0.730.160.08

2.57

⑦1.00
0.54
0.16
0.24
3.29
３つの評価指標（FAC2，FB，NMSE）について，計算ケースごとに，地表煙軸濃度データ及び全地表濃度データについて，評価基準値との比較を示す．ここで、FAC2 の基準値については COST732 に記載された値を使用し、FB と NMS については COST732 で参照されている Chang らの結果(11)を使用した．

三菱重工技報 Vol.49 No.1 (2012)

80

4.3 単純建屋実験との比較
最後に，COST732 に示される風洞実験結果を用いて検証計算を行った．その際，数値計算の検証用として公開されているハンブルグ大学で実施された単純建屋(長さ 100mm×幅 150mm×
高さ 125mm)周りの風洞実験データ(12)を用いた．ただし，この実験データは前流気流としての乱流境界層を生成する地表面粗度条件が公開されていないため，非定常計算においては厳密な前流側気流条件の設定が困難である．しかし，建屋近傍後流側では，建屋による乱れの影響が支配的で前流気流の影響は小さいと考えられることから，検証計算では２次元ブロックにより生成される乱流境界層によって風洞実験での前流気流を模擬し，単純建屋周りの拡散計算を実施した．
その結果，図７に示すとおり，建屋後流側の風下方向地表煙軸濃度は風洞実験結果と非常に良く一致していることを確認した．また，地表煙軸濃度に対する評価指標 FAC2 についても，基準値0.89 に対して 1.0 であり，申し分ない予測精度を有していることがわかる．なお，全地表濃度及び空間内濃度（全点濃度）に対しては，評価指標 FAC2 が 0.57 と，基準値 0.54 を若干上回る程度に留まったが，これは上述したとおり，今回使用したデータベースの地表面粗度条件が不明なために，風洞実験で生成された乱流境界層を数値計算で十分に再現できなかったことが原因と考えられる．

比較対象データ

FAC2 （評価基準）

地表煙軸濃度1.00 （0.89）

全点濃度

0.57 （0.54）

図７ 地表煙軸濃度分布の比較
無次元風下距離 X/Hb に対する無次元化濃度 K＝UCHb2/Q の分布を実験と計算で比較している．
ただし，U：風速，C：地表煙軸濃度，Hb；建屋高さ，Q：放出量である．

|5. まとめ
複雑地形を対象に，乱流モデルの１つの手法であるラージ・エディ・シミュレーション（LES）を用いた非定常気流拡散計算を行い，その計算結果の精度検証を，欧州にて検討されている COST732 の評価手法に従って実施した．最初に，平地条件にて風洞実験結果を十分に再現できることを確認し，次に複雑地形を対象に，地表煙軸濃度データ及び全地表濃度データに対して検証を行った．その際，３つの評価指標(FAC2，FB，NMSE)が，COST732 で推奨されている評価基準を十分に満足していることを示すとともに，別の研究機関で実施された実験データに対しても，評価基準を満足していることを示した．この結果，今回開発した解析技術が十分な予測精度を有しており，風洞実験を代用できるものであることが確認できた．今後は，風洞実験の代替ツールとして，本解析技術の活用によって，検討期間が今までの風洞実験に比べて大幅に短縮できることが期待できる．
謝辞
研究実施にあたり，風洞実験データをご提供いただきました日本原子力発電（株）殿に感謝の意を表します．

三菱重工技報 Vol.49 No.1 (2012)

81

参考文献
(1) （社）日本原子力学会，日本原子力学会標準 発電用原子炉施設の安全解析における放出源の有効高さを求めるための風洞実験実施基準：2009（2010）(2) （社）日本原子力学会 標準委員会 公衆審査 発電用原子炉施設の安全解析における放出源の有効高さを求めるための数値モデル計算実施基準（案）20XX
http://www.aesj.or.jp/sc/public/public-60.html(3) 佐田幸一ほか，安全解析のための数値モデルによる排ガス大気拡散予測及び被ばく線量評価－風洞実
験 と 数 値 モ デ ル に よ る 放 出 源 の 有 効 高 さ ・ 被 ば く 線 量 － ， 日 本 原 子 力 学 会 和 文 論 文 誌 vol.8
No.2（2009） p.184～196(4) 佐藤ほか，低層市街地における建物屋上煙源からの排ガス拡散予測（その２）RANS による濃度予測，第49 回大気環境学会年会 講演要旨集（2008）(5) 道岡，佐藤：低層市街地における建物屋上煙源からの排ガス拡散予測（その３）LES による濃度予測，第49 回大気環境学会年会 講演要旨集（2008）(6) Nakayama，H． et al．，Development of Local-scale High-Resolution Atmospheric Dispersion Model
Using Large-Eddy Simulation Part2: Turbulent Flow and Plume Dispersion around a Cubical Building，Journal of Nuclear Science and Technology，vol.48 No.3 (2011) p.374～383(7) Uchida, T. et al.，Large-eddy Simulation of Turbulent Airflow over Complex Terrain, Journal of Wind
Engineering Industrial Aerodynamics, vol.91 (2003) p.219～229(8) Carissimo, B. et al., The SMEDIS Database and Validation Exercise, International Journal of Environmentand Pollution, vol.16 (2001) p. 614-629(9) Britter, R. et al., Model Evaluation Guidance and Protocol Document, COST Action 732 ，(2007) http://www.cost.esf.org(10) 永山ほか，複雑地形上へ適用した大気拡散数値シミュレーション（LES）結果の風洞実験との比較検証について，日本航空宇宙学会西部支部講演会（2011）講演集 （2011） p.39～42(11) Chang, J. et al., Air quality model performance evaluation, Meteorology and Atmospheric Physics, vol.87(2007) p.167-196(12) ハンブルク大学, 公開データベース http://www.mi.uni-hamburg.de/Data-Sets.432.0.html