Article
Uncertainty Analysis of a Numerical Model for Calculating an Effective Source Height
2022 Volume 21 Issue 3 Pages 127-143
Details
2022 Volume 21 Issue 3 Pages 127-143
The verification and validation processes including uncertainty analysis are now mandatory when utilizing a numerical model in safety analysis. The uncertainty analysis of a numerical model for calculating an effective source height has not been carried out. One of the reasons is a lack of experiment data with uncertainty obtained by the repetition of experiments. In this study, we repeated wind tunnel experiments to obtain uncertainty-quantified experiment data to validate the numerical model under various terrain conditions. Using the obtained wind tunnel data, the uncertainty analysis of the numerical model was carried out to quantify an estimated error of the numerical model. We suggested an adjustment method for the practical application of the numerical model using the estimated error. This method allows us to evaluate an effective source height without the need for wind tunnel results under the same experimental conditions.
放出源の有効高さとは，ガウスプルームモデルを用いた大気拡散計算の際に建屋や地形による拡散への影響を考慮するために用いられるパラメータであり，「発電用原子炉施設の安全解析に関する気象指針1）」（以下，気象指針）では風洞実験によってその妥当性を評価することとされている。風洞実験では，まず建屋や地形のない条件（平地条件）において複数の放出源高度を設定して地表濃度分布の測定を行う。次に建屋や地形の模型を設置した条件（模型条件）において地表濃度分布の測定を行い，得られた濃度が平地条件のどの放出高度に相当するかを評価して放出源の有効高さとする。これらの手順や，実施時の風洞内気流に対する要求基準などが「発電用原子炉施設の安全解析における放出源の有効高さを求めるための風洞実験実施基準：2019（AESJ-SC-P003：2019）2）」（以下，風洞実験実施基準）にまとめられている。また，数値流体力学（CFD）を用いて風洞実験を数値計算で再現することによって放出源の有効高さを評価する手法も提案されており，その手順や計算コードが満たすべき基準などが「発電用原子炉施設の安全解析における放出源の有効高さを求めるための数値モデル計算実施基準：2011（AESJ-SC-A004：2011）3）」（以下，数値モデル実施基準）にまとめられている。CFDを用いた放出源の有効高さ評価について，佐田ら4）は一般座標系に基づいたReynolds Averaged Navier-Stokes（RANS）モデルを使用した気流計算と，Thomson5）による粒子追跡型拡散モデルを組み合わせた手法を提案しており，個々の風向に対する放出源の有効高さの値は風洞実験の結果と多少の差異がみられるものの，常に風洞実験による値よりも大きく，あるいは小さく等，どちらかに偏った値と評価されることはないと報告しており，その結果は数値モデル実施基準にも引用されている。岡林ら6）はLarge Eddy Simulation（LES）に基づいた気流・拡散計算を用いた数値モデルを提案しており，個々の風向に対して風洞実験結果を精度よく再現できることを示した。また，小野ら7）は地表面付近での拡散性状をよりよく再現するためSGS応力の非等方性を考慮したLESに基づく数値モデルを提案し，個々の風向に対する放出源の有効高さおよび地表濃度を精度よく再現できることを示した。
このように，風洞実験によって行われてきた放出源の有効高さ評価をLESに基づいた気流・拡散計算によって代替できる可能性が示されている。しかし，これらの放出源の有効高さを評価するための数値モデルについて不確かさ評価を行った事例はない。原子力分野における数値シミュレーションの利用に際しては「シミュレーションの信頼性確保に関するガイドライン：2015（AESJ-SC-A008：2015）8）」（以下，信頼性ガイド）に示されているように，不確かさ評価を含めた適切なVerification and Validation（V&V）を行うことが望ましい。数値モデル実施基準においても数値モデルのV&Vについて規定がなされており，平地条件を対象としたVerificationおよび単体建物を対象とした濃度分布に対するValidation，平地条件と模型条件の結果から得られた放出源の有効高さに対するValidationについて言及されているが，不確かさ評価については今後検討すべき事項として残されている。数値モデルの不確かさ評価を行うには，比較対象となる風洞実験についても不確かさを定量化する必要があり，放出源の有効高さ評価に関してはそのような風洞実験データが存在していなかったことも数値モデルの不確かさ評価が行われなかった要因であると考えられる。
そこで本研究では，まず新たに風洞実験を実施して風洞実験結果の不確かさについて定量化する（II）。そして，その風洞実験結果を用いて数値モデル実施基準に示されたV&Vのプロセスを実施して数値モデル精度検証を実施する（III）。さらに，数値モデルに関する不確かさ評価と風洞実験の不確かさ評価を行い，数値モデルの推定誤差を定量化するとともに，推定誤差を活用した放出源の有効高さ評価方法を提案する（IV）。
風洞実験は，(一財)電力中央研究所が保有する乱流輸送モデリング風洞の第1試験セクションを用いて実施した。回流式の大型拡散風洞で，試験セクションの寸法は長さ17 m，幅3 m，高さ1.7 mである。測定の概要をFig. 1に示す。実験は風洞実験実施基準に準拠した内容で実施した。試験セクションの風上側には風洞内気流性状を調節するための粗度要素を配置し，これにより後述する所定の条件を満足する気流に調節した。試験セクション中央部には縮尺1：2,000で再現された地形模型および発電所構内の建屋模型を設置し，実験条件ごとに定められた放出源高度の位置に設置されたΓ（ガンマ）型煙突から空気とエチレンの混合ガスをトレーサーガスとして放出した。トラバース装置に設置した水素炎イオン化検出器（FID）によってエチレン濃度を計測した。気流調整等の際の風速計測にはレーザードップラー流速計（LDV）を用いた。
Overview of a wind tunnel experiment
放出源の有効高さへの影響の大小という観点から地形起伏を定量化し，典型的な地形および極力影響の大きい地形を不確かさ評価の対象とするため，日本全国の原子力立地地点の地形について感度解析を実施した。風向は通常の放出源の有効高さ評価と同様に，16方位のうち風下側に陸地が存在する方位とした。各地点・各風向において風上2.5 kmから放出源位置までの風上側領域と，放出源位置から風下5 kmまでの風下側領域における標高，標高の勾配，標高のラプラシアンを算出し，それらの最大値，最小値，平均値，標準偏差などの統計量を求めた。また，数値モデル（ただし計算負荷削減のため水平方向解像度を2倍とした）によって各地点・各風向の地形を再現した計算を行い，放出源の有効高さ（本来は風下の評価範囲内における最低値をとるが，ここでは距離ごとに算出した値）やその風下方向の勾配，濃度分布中心軸のずれについて最大値，最小値，平均値，標準偏差などの統計量を求めた。ここでは，地形による影響に主眼をおき，建屋は全地点共通で70 m角の立方体1つのみとした。放出源高度は後述する実験条件とおおむね同程度となる80 mおよび200 mとした。これらの計算結果から地形に関連する諸量と放出源の有効高さに関連する諸量との相関関係を調べた。ここでは相関係数の絶対値が0.75以上の場合を強い相関があるとした。放出源から風下側の領域を対象とした場合に強い相関のあった組をFig. 2に，放出源から風上側の領域を対象とした場合に強い相関のあった組をFig. 3に示す。なお，H80は放出源高度80 mの場合，H200は放出源高度200 mの場合の放出源の有効高さを表し，MAXは最大値，MINは最小値，MEANは平均値，STDは標準偏差を意味する。図中のRは相関係数を表す。
Correlation between effective source height and topography of downwind side
Correlation between effective source height and topography of upwind side
地形の影響が全くない場合には放出源の有効高さは風下方向に変化しないため，放出源の有効高さの標準偏差や，放出源の有効高さの勾配の標準偏差は地形による影響を示すよい指標となると考えられる。風下側地形については標高の最大値・標準偏差や地形勾配の平均値がこれらと相関が強く，地形選定の際に重視する項目とした。風上側については地形勾配の標準偏差と放出源の有効高さの最小値に強い負の相関がみられた。地形勾配の標準偏差は地形起伏の複雑度を表現する指標となり，風上側に複雑な地形が存在することで放出後すぐから拡散が促進され，放出源の有効高さの最小値の低下につながっていると考えられる。したがって風上側地形については地形勾配の標準偏差を重視する項目とした。Fig. 4にこれらの指標のヒストグラムおよび累積度数分布を示す。
Histgrams of topographical statistics
まず，典型的な地形を想定して各項目がおおむね最大頻度となる実サイトを選定し，そのうち海岸線におおむね直交する風向をCase A，海岸線におおむね平行する風向をCase Bとした。次に，風下1 kmまでの地形および建屋をCase A，Bと共通とした上で，各項目がおおむね上位5%以内となるような仮想地形を風上側・風下側のそれぞれで作成した。この仮想地形を風下側のみ配置したケースをCase C，風上側・風下側の両方配置したケースをCase Dとした。各ケースの地形をFig. 5に示す。本論文における風洞実験および数値モデルでは風下7 kmまでの地形を再現しているが，放出源の有効高さ評価は風下5 kmまでで行われる2,3）。評価対象範囲内におけるCase A，Bの風下側地形はCase C，Dに比べるとほぼ平坦な形状となっている。
Terrain surface models
各ケースの放出源周辺の地形および建屋をFig. 6に示す。Case C，Dでは風下1 km付近においてCase Aの地形と滑らかに接続するように風下側仮想地形を配置した。Case Dの風上側地形は，Case Aの陸上部分にあたる箇所はCase Aと同一とし，Case Aの海上部分にあたる箇所に新たに陸地を配置した。なお，Case A，Bの風洞実験とCase C，Dの風洞実験は別々に実施したため，気流調整および平地条件の濃度測定もそれぞれに対応して行った。
Wind direction and alignment of buildings
風洞実験実施基準では，地形模型および建屋模型を設置しない平地条件において，風洞内気流について次の要件を満たすことを規定している。
また，この気流条件のもとで拡散実験を行い，水平方向の拡がりのパラメータσyおよび鉛直方向の拡がりのパラメータσzが，気象指針の第1図y方向の拡がりのパラメータおよび第2図z方向の拡がりのパラメータに示される大気安定度C～Dの値になっていることを確認することを規定している。
後述するように，各濃度測定は9回の繰り返し測定を行っており，3回に1回は試験セクション内の模型や粗度要素すべてを取り除き，再度配置した。気流調整は粗度要素の再配置のたびに実施したため，1組の実験に対して3回行った。Fig. 7にCase A，BとCase C，Dのそれぞれについて風速分布に関する結果を示す。平均風速は高度800 mにおける平均風速で規格化しており，乱流強度は風速変動の標準偏差をその位置における平均風速で除した値である。平均風速の鉛直分布はべき指数1/7のべき乗則によく従っており，境界層高さも400 m以上となった。主流方向の乱流強度についても高度30 mにおいて13～16%の範囲内に収まる結果となった。
Wind profiles at source location
Figure 8にCase A，BとCase C，Dのそれぞれについてトレーサーガスの水平方向の拡がりのパラメータσyおよび鉛直方向の拡がりのパラメータσzを示す。図中のA～Fは気象指針に示された大気安定度階級別の拡がりのパラメータである。いずれの結果も，大気安定度CとDの間に収まっていることを確認した。
Horizontal and vertical spread of plume concentration
地形模型および建屋模型のない平地条件において，放出源高度HSを0 mから280 mまで変化させてトレーサーガスを放出し，その地表煙軸濃度を風下方向に測定した。複数回の濃度測定におけるばらつきによる不確かさを定量化するため9回の繰り返し測定を行い，測定された濃度の平均値Cと平均値に付随する不確かさu(C)を算出した。放出源の有効高さを評価する風洞実験では1～3回程度の測定を実施して代表値や平均値を用いることが一般的であるが，ここでは（1）式のように9回実施した測定結果の平均値（標本平均）を風洞実験値として用いることとした。したがって，（2）式で計算される標本平均の標準誤差が風洞実験結果（平均値）の不確かさとなる。なお，9回の測定回数で評価される不確かさの不確かさは約25%であり9），後述する不確かさの評価結果の大きさを踏まえると十分な繰り返し回数であると判断できる。また，試験セクション内に設置する粗度要素配置のばらつきも不確かさ要因として考慮し，9回の繰り返しのうち，3回ごとに粗度要素を含むすべてを撤去して再配置を行った。
| \begin{equation} C = \frac{1}{9}\sum_{n = 1}^{9}c_{n} \end{equation} | (1) |
| \begin{equation} u(C) = \sqrt{\frac{1}{9(9 - 1)}\sum\nolimits_{n = 1}^{9}(c_{n} - C)^{2}} \end{equation} | (2) |
ここで，cnはそれぞれの繰り返しn回目における濃度である。また，Uを境界層上端平均風速，Qをトレーサーガス流量，rSを縮率として，基準化濃度CNormを次式で定義した。以降，濃度を示す場合ことわりがなければ基準化濃度を表すものとする。
| \begin{equation} C_{\textit{Norm}} = \frac{CU}{Q}r_{S}^{2} \end{equation} | (3) |
Figure 9にCase A，BとCase C，Dのそれぞれについて平地条件の地表煙軸濃度分布を示す。エラーバーは2u(C)の範囲を示す。放出源高度0 mの場合が最も濃度が高く，放出源高度があがるにつれて濃度は減少する。放出源高度40 m以下では測定対象の風下距離範囲内では排気筒位置から遠ざかるにつれて単調に減少する分布となり，放出源高度60 m以上では測定対象範囲内にピークをもつ分布となった。
Maximum ground level concentration in flat condition
地形模型および建屋模型を設置した模型条件において，放出源高度HSを73 m，150 m，200 mの3通りに変化させてトレーサーガスを放出し，その地表煙軸濃度を風下方向に測定した。放出源高度73 mは想定事故時，150 mおよび200 mは平常運転時をそれぞれ想定しているa）。平地実験と同様にそれぞれの風下距離において9回繰り返し測定を行い，基準化濃度の平均値と平均値に付随する不確かさを算出した。Fig. 10に模型条件の地表煙軸濃度分布を示す。エラーバーは2u(C)の範囲を示す。放出源高度73 mの場合は，すべてのケースにおいて測定範囲内で風下距離とともに単調減少する分布となり，地形起伏の小さいCase A，Case Bはほぼ一致した分布で濃度が高く，Case C，Case Dの順に地形起伏の増加に伴い濃度が低くなる結果となった。放出源高度が低いため建屋による巻き込みの影響を受けて煙軸高度はほぼ地表面と同程度になると考えられ，地形起伏の増加による乱流拡散の増加が煙軸濃度の低下に寄与したと考えられる。放出源高度150 m，200 mの場合は測定対象範囲内にピークをもつ分布となり，Case A，Case Bはほぼ一致し，Case C，Case Dの順に地形起伏の増加に伴ってピーク濃度は高く，ピーク位置は排気筒設置位置に近くなる結果となった。放出源高度が比較的高いため建屋による巻き込みの影響が限定的となり，煙軸高度が維持されて上空を移流拡散していると考えられ，地形起伏の増加による乱流拡散の増加がプルームの着地を早めたと考えられる。
Maximum ground level concentration in terrain and buildings condition
ある風下距離において模型条件で得られた濃度CMが，同じ風下距離における平地条件で得られた2つの放出源高度HSL，HSHにおけるそれぞれの濃度CL，CHの間にあるとすると，その風下距離における放出源の有効高さHは以下の式より評価できる。
| \begin{align} H &= f(C_{M},C_{L},C_{H}) \\ &= \sqrt{H_{\textit{SL}}^{2} + \frac{\ln C_{L} - \ln C_{M}}{\ln C_{L} - \ln C_{H}}(H_{\textit{SH}}^{2} - H_{\textit{SL}}^{2})} \end{align} | (4) |
これを風下距離5 kmまで評価し，敷地境界または周辺監視区域以遠における最も小さなHの値を放出源の有効高さとする。なお，本来は最終的に5 m単位に丸めた結果を放出源の有効高さとするが2），本報では（4）式の結果を丸めずに示している。
また，放出源の有効高さの不確かさu(H)については，放出源の有効高さが決定された風下距離における模型条件および平地条件の濃度の不確かさu(CM), u(CL), u(CH)から不確かさの伝播則9）を用いて以下の式より評価できる。
| \begin{align} &u(H) \\ &= \sqrt{\left(\frac{\partial f}{\partial C_{M}}\right)^{2}u^{2}(C_{M}) + \left(\frac{\partial f}{\partial C_{L}}\right)^{2}u^{2}(C_{L}) + \left(\frac{\partial f}{\partial C_{H}}\right)^{2}u^{2}(C_{H})} \cr \end{align} | (5) |
Figure 11に風洞実験による放出源の有効高さ評価結果を，Fig. 12に風洞実験による放出源の有効高さの不確かさの評価結果を示す。図中のCaseはCase A～Dのいずれかを，Hsは放出源高度を示す。放出源高度73 mの場合は地形の起伏が大きくなるにつれて放出源の有効高さが大きく，放出源高度150 mおよび200 mの場合は地形の起伏が大きくなるにつれて放出源の有効高さが小さくなった。これは前述した地表煙軸濃度分布の傾向とも整合している。
Effective source height from wind tunnel experiment
Uncertainty in wind tunnel experiment
放出源の有効高さの不確かさについては最大で約3 mであり，実際の評価では放出源の有効高さを5 m単位で丸めることを考慮すると十分に小さいと考えられる。放出源高度に対する明確な変化はなく，地形に対する変化としては，比較的地形起伏の小さいCase A，Bでは大きく，比較的地形起伏の大きいCase C，Dでは小さい傾向にあった。地形起伏の小さいケースではプルームの拡がり幅が小さいため空間的な濃度勾配が大きく，計測ごとのばらつきが生じやすいと考えられる。
7. 地形による風洞実験の不確かさの変化対象とする地形によって風洞実験の不確かさにどの程度の変化があるのかを明らかにするため，地形模型を一切設置せずに建屋模型についても単独の円柱のみを使用した実験を実施した。円柱模型は直径45 m，高さ65 mとし，放出源高度は平常運転時を想定した162.5 mとした。Case A～Dと同様に9回の繰り返し測定を行い，放出源の有効高さの不確かさを評価した。なお，風洞内の気流はCase A～Dの実施時と同様に粗度要素を配置して乱流境界層を発達させている。円柱模型設置位置における乱流境界層のレイノルズ数は約1.6 × 105であり，拡散実験に対するレイノルズ数依存性は小さいと考えられる10）。Figure 13に不確かさの評価結果をCase A～Dの平常運転時（放出源高度150, 200 m）の結果と合わせて示す。なお，ここでは放出源の有効高さが決定する風下距離以外についても算出している。ほぼすべての風下距離において，地形なしの結果は地形あり（Case A～D）の結果の範囲内に収まっており，大きな差異がないことが確認できる。この結果から，本報で使用した風洞設備および実験条件下において求めた放出源の有効高さの不確かさは，対象とする地形によって大きく変化せず，おおむね0.5～3 m程度であるといえる。
Comparison of uncertainty in wind tunnel experiment
数値モデルは小野ら7）の提案した非構造格子の有限体積法に基づくLES数値モデルを使用した。この数値モデルではサブグリッドスケールモデルとしてAbe11）の非等方効果を考慮した1方程式モデルを採用している。本報では，ベースとなる格子幅を主流方向に5～10 mm（実スケール換算で10～20 m），スパン方向に5 mm（実スケール換算で10 m），鉛直方向に1～20 mm（実スケール換算で2～40 m）とした。地形の再現範囲は実スケール換算で風下7 kmまでとし，風下7 km以降は風下7 km位置の地形を風下方向に0.5 km伸ばすことで流出境界において地形影響による逆流が生じないように配慮した。建屋の近傍では主流方向およびスパン方向にそれぞれ格子を2分割する細分化を2段階施している。時間刻みは0.00025 sとし，平均化時間は120 sとした。その他は文献7）と同様で，計算領域は計算負荷の削減を目的に風洞実験よりやや小さいスパン方向2 m・鉛直方向1 mとした。地表面・建物表面にはSpalding則を利用した対数則を用い，側面および天井面はslip条件とした。始めに粗度要素をすべて再現した計算を実施して煙源の1.5 mおよび2.5 m前方の断面における風速の時系列データを保存し，これを流入条件とした。流出条件は圧力を0とし，その他は勾配0条件とした。
2. モデル検証一般的にモデル検証では，計算コードで使用している方程式の離散化，収束判定や計算アルゴリズムの選択などが適切かどうかを確認する。加えて数値モデル実施基準では，「大気安定度がほぼ中立の気流および拡散状態を再現できることが放出源の有効高さを評価するうえで必須の条件である」とし，平地条件の計算において以下の気流設定条件および拡散設定条件のそれぞれを満たすことを確認することがモデル検証として定めている。
（1） 気流設定条件気流設定に関して，満たすべき条件は以下の3つである。
なお，これらは風洞実験実施基準における気流調整のための要件とほぼ同等のものであるが，前述したように風洞実験実施基準では主流方向の乱流強度は13～16%とやや厳しい条件になっており，ここでは風洞実験実施基準の値に従うこととした。Fig. 14に排気筒位置における平均風速および主流方向の乱流強度の鉛直分布を示す。風洞実験と同様に平均風速は高度800 mにおける平均風速で規格化しており，乱流強度は風速変動の標準偏差をその位置における平均風速で除した値である。上記の気流設定条件を満足していることが確認できる。
Wind profiles at source location
数値モデル実施基準では，（1）の気流条件のもとに計算領域内を平地とした拡散計算を行い，「計算対象範囲における鉛直方向の拡がりのパラメータσzが気象指針の第2図z方向の拡がりのパラメータに示される大気安定度C～Dに対応するほぼ中立の値になっている」ことを確認することが定められている。また，風洞実験実施基準では上記に加えて，水平方向の拡がりのパラメータσyについても大気安定度C～Dに対応することが追加されている。
Figure 14に示した気流条件のもとに拡散計算を行った結果をFig. 15に示す。水平方向，鉛直方向ともに，拡がりのパラメータが計算対象範囲である風下距離5 kmまで大気安定度C～Dに対応する分布となっており，数値モデル実施基準および風洞実験実施基準の要件を満たしていることがわかる。文献7）では再現対象とした風洞実験において，大気安定度C～Dに対応する範囲が風下距離3 km程度までであり，また風洞の試験セクション外の凹凸も試験セクション内の気流形成に寄与していたため数値モデルでの再現が難しい点があった。本論文で実施した風洞実験は風下距離5 kmまで大気安定度C～Dに対応しているうえ，試験セクション外の凹凸を極力小さくし，試験セクション内に影響を及ぼさないように配慮した。そのため，試験セクション内の気流に影響を及ぼす粗度要素をすべて数値モデル上で再現できたことにより，良好な再現結果が得られたと考えられる。
Horizontal and vertical spread of plume concentration
数値モデル実施基準では，数値モデルの妥当性確認として「直方体建屋の風洞実験結果による妥当性確認」および「放出源の有効高さの風洞実験結果による妥当性確認」の2つの手順を規定している。前者は単純な形状の建屋のみ考慮して地形を考慮しないケースにおける濃度予測結果の評価，後者は実際の発電所相当の建屋・地形を考慮したケースにおける放出源の有効高さ予測結果の評価となっている。
（1） 直方体建屋の風洞実験結果による妥当性確認数値モデル実施基準では，「直方体建屋等の単純形状の建屋対象にした風洞実験結果との比較による数値モデルの妥当性確認を行う」として，CEDVAL A1-5b）の風洞実験結果を用いた妥当性確認を参考として掲載しており，全点濃度の比較においてFAC2c）が0.54以上，煙軸濃度の比較においてFAC2が0.89以上という基準が示されている。本報においても同様にCEDVAL A1-5を使用して妥当性確認を行った。CEDVAL A1-5は，主流方向の幅20 m，スパン方向の幅30 m，高さ25 mの直方体建屋を縮尺$1:200$で再現された模型を風洞内に設置し，建屋模型の風下側壁面に設けられた4つの排気口からトレーサーガスを放出させた実験である。
計算領域をFig. 16に示す。CEDVALでは粗度要素としてスパイアとブロックを使用したと記載があるが，その大きさや配置に関しては示されていない。そこで実験写真等を参考にしたうえで，直方体建屋の1 m手前で計測された風速鉛直分布が風洞実験結果と整合するようにスパイアとブロックの大きさ，配置を決定した。
Computational region for validation test of CEDVAL A1-5
Figure 17に風洞実験結果との風速鉛直分布の比較結果を示す。風洞実験では2成分同時測定が可能な風速計を使用し，主流方向（U）および鉛直方向（W）の同時測定と，主流方向およびスパン方向（V）の同時測定が行われており，その両方の結果を示している。主流方向の平均風速については，やや風洞実験結果よりも鉛直勾配が大きいものの，おおむね風洞実験の傾向を再現できている。速度変動（風速変動の標準偏差）については，スパン方向でやや過大評価であるが，主流方向および鉛直方向については風洞実験結果とおおむね一致する結果が得られた。全点濃度の比較結果をFig. 18に示す。なお，無次元濃度Kは以下のように定義されている。ここでは，Cmeasuredは各地点の計測濃度，Csourceは放出濃度，Qは放出強度，Urefは上空風速，Hbは建物高さである。
| \begin{equation} K = \frac{C_{\textit{measured}}}{C_{\textit{source}}Q}U_{\textit{ref}}H_{b}^{2} \end{equation} | (6) |
Wind profiles of approach flow
Point-by-point comparison of concentration
無次元濃度で0.1以下ではFAC2範囲を超えてばらつく点がみられるが，無次元濃度で1以上の点は多くがFAC2範囲内に収まる結果となった。FAC2の値も0.87であり，基準を上回った。地表煙軸濃度の比較結果をFig. 19に示す。すべての点でFAC2範囲内，すなわちFAC2は1であり，基準を上回ることを確認した。
Comparison of maximum ground level concentration
次に，数値モデル実施基準では「数値モデルを適用する発電所において風洞実験および数値モデルから評価された放出源の有効高さを比較し，以下に示すa)およびb)を満足していることを確認する」としており，ここではII章のCase A～Dの風洞実験結果を用いて比較を行う。
a) 数値モデルにより求めた放出源の有効高さと風洞実験により求めた放出源の有効高さの回帰直線の傾きが0.9～1.1倍でできるだけ1.0に近く，相関係数は0.9以上であること。
b) 次式より求めた変動係数dが20%以内であること。
| \begin{equation} \begin{split} &d = \sqrt{\frac{S_{yy} - 2S_{xy} + S_{xx}}{n - 2}} \times \frac{1}{{\overline{Y}}} \times 100\\ &S_{xx} = \sum(X - \overline{X})^{2}\\ &S_{yy} = \sum(Y - \overline{Y})^{2}\\ &S_{xy} = \sum(X - \overline{X})(Y - \overline{Y}) \end{split} \end{equation} | (7) |
ここでは記号の定義は数値モデル実施基準に拠っており，Xは各ケースの数値モデルによる放出源の有効高さ評価値，Yは各ケースの風洞実験による放出源の有効高さ評価値，$\overline{\phantom{\mathrm{M}}}$は全ケースにおける平均操作，nはケース数である。
数値モデルによって平地条件および模型条件の濃度計算を行い，風洞実験と同様に（4）式によって放出源の有効高さを評価した。数値モデルによる，平地条件における地表煙軸濃度分布の計算結果をFig. 20に，模型条件における地表煙軸濃度分布の計算結果をFig. 21に示す。平地条件では風洞実験結果と同様に放出源高度0 mの場合が最も濃度が高く，放出源高度があがるにつれて濃度は減少する傾向が再現できた。模型条件においても，地形起伏の小さいCase AとCase Bがほぼ同程度の濃度を示し，地形起伏の増加によって濃度分布形状が変化するという風洞実験結果の傾向を再現した。なお，数値モデル計算結果では濃度分布に凹凸がみられるが，これは小規模な地形起伏による濃度増減の影響を反映しているためである。模型条件における地表濃度分布の計算結果をFig. 22に示す。数値モデルおよび風洞実験においてそれぞれ評価された放出源の有効高さを比較した結果をFig. 23に示す。図中には回帰直線とその傾き，相関係数，変動係数を合わせて示す。いずれも数値モデル実施基準に定められたa)およびb)の基準を満足していることを確認した。
Maximum ground level concentration in flat condition
Maximum ground level concentration in terrain and buildings condition
Ground level concentration in Hs = 73 m
Comparison of effective source height
信頼性ガイドでは，「測定量またはシミュレーション結果と不可知である真値との差」として誤差を定義し，これを推定する手法，すなわち推定誤差を求める手法について解説している。数値モデルによるシミュレーション結果をyS，実験データをyDとすると，両者の偏差Ecomparison (= yS − yD)を用いて，モデルの推定誤差$\tilde{\delta }_{\textit{model}}$を次式から評価できるとしている。
| \begin{equation} \tilde{\delta }_{\textit{model}} = E_{\textit{comparison}} \pm u_{\textit{val}} \end{equation} | (8) |
ここで，uvalは
| \begin{equation} u_{\textit{val}} = \sqrt{u_{\textit{num}}^{2} + u_{\textit{input}}^{2} + u_{D}^{2}} \end{equation} | (9) |
| \begin{equation} \tilde{\delta }_{\textit{model}} = E_{\textit{comparison}} \pm ku_{\textit{val}} \end{equation} | (10) |
数値モデルによる計算結果は必ず誤差を含んでおり，実際の安全解析で用いる場合にはその誤差の程度を明らかにして必要ならば調整を行うことが望ましい。数値モデル実施基準では，数値モデルが風洞実験を代替するものとして開発されてきたことを鑑み，数値モデルによって計算された放出源の有効高さが，95%の確率で風洞実験による評価結果を上回ることのないように調整する方法が示されている。各ケースの風洞実験による評価結果をY，数値モデルによる評価結果をXとすると，調整後の数値モデルによる評価結果XCSを（11）式から算出する。
| \begin{equation} \begin{split} &X_{\textit{CS}} = X\left(\frac{Y}{aY + b + 2\sigma_{x}}\right)\\ &\sigma_{x} = \sqrt{\frac{\displaystyle\sum(x - \overline{x})}{n - 1}}\\ &x = X - Y \end{split} \end{equation} | (11) |
ここでは記号の定義は数値モデル実施基準に拠っており，$\overline{\phantom{\mathrm{M}}}$は全ケースにおける平均操作，nはケース数，aは回帰直線の傾き，bは回帰直線の切片である。
この方法では，調整を行う際に各ケースの風洞実験結果が必要となり，新設・増設や大幅な敷地内建屋改変等に伴う有効高さ評価への数値モデル適用が困難である。そこで本報では，これに代わる方法として推定誤差を用いた調整方法を次のように提案する。この方法による調整後の評価結果XEは（12）式で算出する。
| \begin{equation} X_{E} = X - \tilde{\delta }_{\textit{model}} \end{equation} | (12) |
これは，推定誤差がシミュレーション結果と真値との差についての統計的な符号付き推定量として定義される8）ものであることから，数値モデルによる計算結果と推定誤差を用いて放出源の有効高さの真値を推定しようとするものである。なお，推定誤差$\tilde{\delta }_{\textit{model}}$は（10）式に示されたように信頼区間に応じた範囲をもつ推定量であるため，（12）式で真値推定結果として評価される調整後の評価結果XEも範囲をもつことになる。安全解析における地表付近の濃度評価を考えると放出源の有効高さが小さい方が濃度は高くなるため，調整後の評価結果XEのうち最も小さい値を最終的な評価結果とすることで保守性を有した評価が可能と考えられる。信頼水準を95%としたときの下限値（95%LCL）は以下となる。
| \begin{align} X_{E,95\% \textit{LCL}} &= X - \tilde{\delta }_{\textit{model},95\% \textit{LCL}} \\ &= X - (E_{\textit{comparison}} + 2u_{\textit{val}}) \end{align} | (13) |
一方で，確率変数としてのuvalがほぼ正規分布に従うとするとその期待値は0となるため，評価結果XEのうち最も真値である確率が高い値（BE）は以下で表される。
| \begin{equation} X_{E,\textit{BE}} = X - \tilde{\delta }_{\textit{model},\textit{BE}} = X - E_{\textit{comparison}} \end{equation} | (14) |
ここでは，III章にて検証・妥当性確認を行った数値モデルについて，II章の風洞実験データを用いて推定誤差の評価を行った。推定誤差を構成する各要素のうち，数値モデル計算結果と風洞実験結果の偏差Ecomparison，風洞実験結果の不確かさuDについてはすでにII章およびIII章にて評価済みである。また，流入風速や放出強度などの入力パラメータについてはすべてのケースで同一で，ケース間で異なるのは地形・建屋配置の違いによる計算メッシュの差異のみであるため，離散化等に伴う不確かさunumに内包されると考えてuinput ≈ 0とした。
離散化等に伴う不確かさunumについては，信頼性ガイドおよびASME V&V20を参考に，Grid Convergence Index（GCI）を用いて評価した。GCIは，同一のケースに対して複数の格子解像度で計算を行い，その計算結果から格子解像度に対する解の収束性を判断するものである。GCIは以下の手順で評価される14,15）。
a) 複数の異なる解像度をもつ計算格子にて計算を行う。ここでは3段階の解像度とし，それぞれの格子幅を細かい方から順にh1，h2，h3とし，計算結果をϕ1，ϕ2，ϕ3とする。3つの解像度は，格子幅の比r21 = h2/h1，r32 = h3/h2がそれぞれ1.3以上となるように決定することが望ましい。
b) 収束の次数pcを以下の関係式から評価する。
| \begin{equation} \begin{split} &p_{c} = | \ln | \varepsilon_{32}/\varepsilon_{21} | + q(p_{c}) |/\ln r_{21}\\ &q(p_{c}) = \ln \left(\frac{r_{21}^{p_{c}} - s}{r_{32}^{p_{c}} - s}\right)\\ &s = 1\cdot \text{sign}(\varepsilon_{32}/\varepsilon_{21}) \end{split} \end{equation} | (15) |
ただし，ε32 = ϕ3 − ϕ2，ε21 = ϕ2 − ϕ1である。格子幅の比r21およびr32が同一であればq(pc) = 0となるが，それ以外の場合は反復法などで（15）式を解く必要がある。なお，GCIを求める対象が局所的な変数である場合，（15）式で求まるpcが空間的に大きく変動する場合がある。Cadafalch et al.16）はpcの空間平均を用いる方法を提案しており，本報では主流方向の平均化を行った。
c) 求めたpcを用いて，最も細かい計算格子による計算結果に対するGCIを次式より評価する。Fsは安全係数であり，Roache14）は3段階以上の計算格子を用いた場合は1.25とすることを推奨している。
| \begin{equation} \text{GCI} = \frac{F_{s}}{r_{21}^{p_{c}} - 1}| \phi_{1} - \phi_{2} | \end{equation} | (16) |
以上のようにして求められたGCIは，離散化に伴う誤差が0である理想解に対する実際の計算結果の不確かさを95%信頼区間の値として表している。一方で（9）式においてuvalの評価に用いられるunumは1σの値であるため，換算が必要である。ASME V&V20では，最も細かい格子における誤差分布がshifted Gaussianにほぼ従うとして1.15で除すことが示されている。
| \begin{equation} u_{\textit{num}} = \text{GCI}/1.15 \end{equation} | (17) |
数値モデル計算においても風洞実験と同様に，平地条件および模型条件の2種類の計算で地表煙軸濃度を算出し，（4）式を用いて放出源の有効高さを求める。したがって，平地条件および模型条件の各計算でのGCI評価により濃度に対する不確かさを求め，それらの値から（5）式の不確かさの伝播則を用いて放出源の有効高さに対する不確かさを評価した。
GCI評価に使用した計算格子をTable 1に示す。FineはIII章において放出源の有効高さの風洞実験結果による妥当性確認を実施した際と同じものであり，MiddleはFineに対して各方向に1.3倍，CoarseはMiddleに対して各方向に1.5倍とした。ただし，鉛直方向の地表面に隣接する第1層目に関してはサブグリッドスケールのモデリングとも密接に関わるため，Fineと同じ厚さとした。Coarseの場合は第1層目をFineと同じ厚さ，第2層目をMiddleと同じ厚さとし，第3層目以降を本来のCoarseの厚さとした。
| Name | Grid width (as real scale [m]) | ||
|---|---|---|---|
| Streamwise | Spanwise | Vertical | |
| Fine | 10 to 20 | 10 | 2 to 40 |
| Middle | 13 to 26 | 13 | 2 to 52 |
| Coarse | 19.5 to 39 | 19.5 | 2 to 78 |
Figure 24に平地条件および模型条件で計算した地表煙軸濃度に対する不確かさを示す。これをもとに，放出源の有効高さが評価される風下距離において，（5）式を用いて放出源の有効高さに対する不確かさを算出した。結果をFig. 25に示す。図中のCaseはCase A～Dのいずれかを，Hsは放出源高度を示す。算出された不確かさは約0.5～7 mの範囲でばらついた。地形に対する変化はさほど明瞭にはみられないが，放出源高度に対しては，想定事故時を想定した73 mでは平常運転時を想定した150 m，200 mの場合に比べて小さい傾向がみられた。
Uncertainty estimated by means of GCI
Uncertainty due to discretization
以上の結果より，推定誤差を構成する各要素および95%信頼区間を想定した場合の推定誤差の範囲について算出した。結果をTable 2に示す。表中のCaseはCase A～Dの各ケースを，Hsは放出源高度を表す。Ecomparison，unum，uDがケースごとに異なった値を取ることから，推定誤差$\tilde{\delta }_{\textit{model}}$についてもケースごとに差異が生じており，その絶対値は約0.5～20の範囲でばらつく結果となった。
| Case, H | yS [m] |
yD [m] |
unum [m] |
uD [m] |
Ecomparison [m] |
uval [m] |
$\tilde{\delta }_{\textit{model}}$ [m] |
|---|---|---|---|---|---|---|---|
| A, 73 | 46.0 | 36.9 | 0.6 | 0.5 | 9.2 | 0.7 | 7.7 ∼ 10.7 |
| A, 150 | 112.9 | 107.6 | 7.1 | 2.5 | 5.3 | 7.5 | −9.7 ∼ 20.3 |
| A, 200 | 156.9 | 161.2 | 1.5 | 1.9 | −4.3 | 2.5 | −9.3 ∼ 0.6 |
| B, 73 | 42.4 | 46.4 | 2.2 | 1.6 | −4.0 | 2.7 | −9.5 ∼ 1.4 |
| B, 150 | 105.2 | 105.0 | 1.4 | 3.0 | 0.2 | 3.3 | −6.4 ∼ 6.9 |
| B, 200 | 159.5 | 153.5 | 3.8 | 1.8 | 6.0 | 4.2 | −2.4 ∼ 14.4 |
| C, 73 | 57.3 | 49.7 | 0.4 | 0.6 | 7.6 | 0.7 | 6.2 ∼ 9.0 |
| C, 150 | 109.9 | 100.7 | 1.4 | 0.7 | 9.2 | 1.5 | 6.1 ∼ 12.3 |
| C, 200 | 144.0 | 140.6 | 1.1 | 0.8 | 3.4 | 1.4 | 0.7 ∼ 6.2 |
| D, 73 | 77.7 | 74.5 | 1.8 | 0.5 | 3.2 | 1.8 | −0.5 ∼ 6.9 |
| D, 150 | 84.6 | 82.7 | 4.0 | 0.6 | 1.8 | 4.1 | −6.3 ∼ 10.0 |
| D, 200 | 103.9 | 100.4 | 4.2 | 0.9 | 3.6 | 4.3 | −5.0 ∼ 12.1 |
前節で述べた，推定誤差を用いて真値を推定することでモデル評価結果を調整する手法を実用するうえでは，地形など個別地点に左右される条件によらずに数値モデルの推定誤差が定量化されていることが望ましい。そこで，Table 2に示す結果をもとに，本論文で用いた数値モデルにおける推定誤差の代表値を定量化することを試みた。まず，地形起伏の大小によって推定誤差の定量結果に有意な差異が生じないことを確かめるため，比較的平坦で起伏の少ない地形であるCase A，Bグループと，国内立地の中でも上位5%程度の起伏の大きさをもつCase C，Dグループにわけ，“両者の間に推定誤差の差異はない”という帰無仮説を立ててT検定を行った。次に，代表値としてTable 2に示す結果の平均値を用いることの妥当性について検討するため，Case A～Dのすべての結果を対象に推定誤差のばらつきの傾向を調べた。なお，推定誤差$\tilde{\delta }_{\textit{model}}$は，符号付きの成分Ecomparisonと，符号なしの成分uvalからなるので，それぞれの成分について検討を行った。
T検定を適用するにはそれぞれの標本が正規分布に従い，さらに等分散である必要がある。したがって，それぞれのグループについてのShapiro-Wilk検定（正規性検定），グループ間のF検定（等分散性検定）も合わせて行った。いずれの検定も有意水準を5%とした。それぞれの検定について検定統計量の実現量が出現する確率であるp値をTable 3に示す。
| Test | p-value | |
|---|---|---|
| Ecomparison | uval | |
| Shapiro-Wilk (Case A, B) | 0.40 | 0.54 |
| Shapiro-Wilk (Case C, D) | 0.16 | 0.14 |
| F | 0.17 | 0.39 |
| T | 0.32 | 0.32 |
Shaprio-Wilk検定，F検定でいずれも検定結果のp値が0.05以上であることから，どちらのグループも正規性があるという仮説およびグループ間で分散が同程度という仮説は棄却されなかったため，T検定を適用した。また，T検定ではいずれも検定結果のp値が0.05以上であり，グループ間でEcomparisonとuvalの平均値に差がないという仮説は棄却されなかった。すなわち，地形起伏の大小によって推定誤差に有意な差異が生じることは確認されなかった。
次に，Case A～Dのすべての結果を標本とし，Ecomparisonおよびuvalのそれぞれの成分についてShapiro-Wilk検定を行った。帰無仮説は“全データのばらつきは正規分布に従う”であり，その場合Case A～Dの平均値を代表値とできると考えた。また，データの偏りを評価するため，それぞれの標本の歪度を評価した。結果をFig. 26に示す。
Histgrams of uncertainty elements
Ecomparisonおよびuvalの双方とも，正規分布に従うという仮説は棄却されなかった。また，Ecomparisonについては平均値と最頻値がほぼ同じ値をとっており，歪度も0.53と比較的小さいことから平均値がある程度代表性をもつと考えられる。一方で，uvalについては歪度が1.14と大きく，平均値よりも小さい側に偏りをもつ分布である。平均値は最頻値よりも大きな値となっている。しかし，数値モデルによる計算結果の補正に用いるための定量化であることを念頭に置くと，（13），（14）式のように確率変数としてのuvalが取り得る値のうち，0以上の正の値を常に使用することになる。したがって，大きめのuvalを用いることは放出源の有効高さをより小さく評価することとなり，被ばく評価の上では安全側であることから代表値として平均値を用いることとした。
以上の結果より，Ecomparisonとuvalについてそれぞれ平均値を用いて推定誤差$\tilde{\delta }_{\textit{model}}$を評価した。また，評価された$\tilde{\delta }_{\textit{model}}$を用いて（11）式によって調整した放出源の有効高さ評価結果をFig. 27に示す。バーの高さが期待値を示し，ひげが95%信頼区間を示す。なお，（11）式による調整結果も合わせて示す。図中のCase, HsはCase A～Dのいずれかおよび放出源高度の組み合わせを示す。例えばB, 150はCase Bで放出源高度が150 mのケースをさす。
Adjusted effective source height
本報で提案した方法で調整された$X_{E,95\% \textit{LCL}}$は，数値モデル実施基準に示された調整結果XCSとおおむね同程度の値となり，さらに全ケースで風洞実験結果を下回る結果となった。これにより，数値モデル実施基準で示された手法とおおむね同等な保守性を有することが示された。さらに，数値モデル実施基準に示された方法では，評価対象となるケースごとに対応する風洞実験結果が必要となるのに対し，本報で提案した手法では調整に用いる推定誤差を代表値としてあらかじめ定量化することで，風洞実験結果が存在しないケースにおいても評価が可能となる。
放出源の有効高さを評価するための数値モデルについて，シミュレーションの信頼性確保に関するガイドライン：2015やASME V&V20に基づいて不確かさ評価を行った。さらに，不確かさ評価によって得られた推定誤差を用いて，数値モデルによる放出源の有効高さ算出結果から真値の存在範囲を推定し，それに基づいて放出源の有効高さを決定する手法について提案した。不確かさ評価においては，国内立地において典型的な地形特徴をもつ地点および上位5%程度の大きな地形起伏をもつ地点を選定し，両者から得られた推定誤差に明確な差異がみられない点を示し，それらの結果の平均値から推定誤差の代表値を定量化した。これにより従来の方法と異なり，新設・増設や大幅な敷地内建屋改変等の際においてもその条件に対応した風洞実験結果を必要とせずに数値モデルによる放出源の有効高さ評価が可能となり，より効率的な実機適用が期待される。また，この手法で推定された真値のうち95%信頼区間の下限値に相当する値を採用した場合，大きく保守性を見込んだ従来の手法と同程度の結果が得られることがわかり，保守性についても十分確保できているものと考えられる。
