発表内容
OpenFOAMによる大気拡散計算
１．OpenFOAMによる拡散計算のためのソルバー変更
２．煙突ダウンウォッシュへの適用

電力中央研究所

瀧本 浩史
NEE研究会2015/1/27

12スカラー輸送方程式の導入（１）（バージョン2.1.1に対して）
１．ベースとなるソルバを選択（simpleFoam, pimpleFoam, buoyantBoussinesqPimpleFoamなど）→以降はpimpleFoamをベースに話をすすめる。
２．pimpleFoamの構成ファイルは以下の通り。
$FOAM_SOLVERS/incompressible/pimpleFoamにおいて./pimpleFoam.C ・・・ メインルーチン./UEqn.H ・・・ 運動方程式./pEqn.H ・・・ 圧力方程式./createFields.H ・・・ 変数設定、基準点圧力の読み込みなど./Make/options ・・・ 必要なヘッダファイルの場所、ライブラリへのリンクを指定./Make/files ・・・ コンパイルするファイルとソルバ名の指定OPENFOAMによる拡散計算

その他のファイルはビルド時に作成される。

34

スカラー輸送方程式の導入（２）スカラー輸送方程式の導入（３）新しいソルバを作成するため、ベースとなるソルバをディレクトリごとコピーし、別の名前をつける。
$ cd $FOAM_SOLVERS/incompressible$ cp -r pimpleFoam pimpleFoamC

スカラー輸送方程式を解くCEqn.Hを作成。
volScalarField kappact(“kappact”, turbulence->nut()/Sct);
kappact.correctBoundaryConditions();//渦拡散係数、Sctの設定は後述

volScalarField kappacEff(“kappacEff”, turbulence->nu()/Sc + kappact); //拡散係数＋渦拡散係数
fvScalarMatrix CEqn
(fvm::ddt(C)+ fvm::div(phi, C)- fvm::laplacian(kappacEff, C));

メインルーチンのファイル名を変更。
$ cd pimpleFoamC$ mv pimpleFoam.C pimpleFoamC.C

CEqn.relax();
CEqn.solve();./Make/filesのソルバ名を書き換え。
pimpleFoam.C
EXE = $(FOAM_APPBIN)/pimpleFoam↓
pimpleFoamC.C
EXE = $(FOAM_USER_APPBIN)/pimpleFoamC+−+=0//緩和係数の設定

56

スカラー輸送方程式の導入（４）スカラー輸送方程式の導入（５）createFields.Hに以下を追加。

pimpleFoam.Cから読み込む。

Info<< "Reading field C¥n" << endl;
volScalarField C(IOobject(“C",
runTime.timeName(),mesh,
IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);

while (pimple.loop()){#include "UEqn.H"// --- Pressure corrector loop
while (pimple.correct()){#include "pEqn.H"}
if (pimple.turbCorr()){
turbulence->correct();}#include “CEqn.H”

・・・
singlePhaseTransportModel laminarTransport(U,phi);#include “readTransportProperties.H“ //ケースディレクトリのconstant/transportPropertiesの読み込み}78スカラー輸送方程式の導入（６）ケースの設定readTransportProperties.Hを作成
１．ケース/constant/transportPropertiesにシュミット数を追加// Laminar Schmidt number
dimensionedScalar Sc(laminarTransport.lookup("Sc"));Sc
Sct// Turbulent Schmidt number
dimensionedScalar Sct(laminarTransport.lookup("Sct"));Sc [0 0 0 0 0 0 0] 0.9
Sct [0 0 0 0 0 0 0] 0.7

２．ケース/0ディレクトリ内にCの初期条件、境界条件設定ファイルを作成。
ビルド$ wmake

３．ケース/system/fvSchemesで、スカラー輸送方程式の離散化スキームを設定。
４．ケース/system/fvSolutionsで、スカラー輸送に対する代数方程式の解法を設定。

910

その他
・煙突からの放出を考える場合は、煙突の吐出面を固定端境界として濃度を与えればよい。 （濃度×流速×面積が放出量となる）
・空中放出の場合、以下を参考にソース項を実装する必要がある。
今野先生の「輸送方程式のソース項の実装」を参照。
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&cad=rja&uact=8&ved=0CC
MQFjAB&url=http%3A%2F%2Fwww.opencae.jp%2Fraw-attachment%2Fwiki%2FOpenCAE-UserGrouparchive%2F20120421OFBeginner.pdf&ei=OWG_VPL3Ko738QW2gYHYCw&usg=AFQjCNEfKFpQ
O8OBHS8POjxHw339NGHYUg&sig2=dpVxgJZDRpNVfdAaQ8sitQ

煙突ダウンウォッシュ問題への適用1112

背景問題点

発電所からの排ガス拡散の実験・計算
煙突を配置せず水平放出
or煙突を配置し上昇過程を再現→これまでの主流→上昇過程が複雑地形や建物影響を受けるときには，上昇過程を再現する必要がある。火力発電所の低煙突化や山間地に建設される地熱発電所など。

Lateb et al. (2010)

Breuer (2000)
有効高さ・煙突（円筒）周りを真面目に解くのは難しい
・発電所の場合、煙突の径が~10mであるのに対して、計算領域は~10km、粘性底層は~1mm→総格子数、タイムステップの制約から煙突をあまり細かく解像できない。

1314目的煙突ダウンウォッシュ
煙突後背部の低圧部に排ガスが巻き込まれる現象

現状では、煙突周囲の格子サイズ、格子タイプに関する指標はない。
→円筒形の煙突を対象に、複数の異なる格子タイプ、格子解像度を用いて煙突付近の流れ場の再現精度が拡散予測に与える影響を把握する。

Tip vortices

Johnston and Wilson (1996)－ レイノルズ数依存（臨界レイノルズ数2x105）－ 実スケールではあまり問題にならないが、風洞実験では重要（風洞実験結果を使ったCFDモデルの精度検証も）－ Ws/Us < 1.5 のときに顕著な影響（亜臨界のレイノルズ数において）
Ws: 排ガス速度, Us: 煙突高さの水平風速

1516計算概要

実験による上昇過程の再現
話は変わりますが・・・

15D

風洞実験においても、発電所周囲の排ガス拡散予測については，有効高さを予め計算し，水平方向に放出する方法が用いられることが多い。
→ 主にレイノルズ数とフルード数による制約のため。
（排煙レイノルズ数，円筒レイノルズ数）

Park and Lee (2002)の実験を参考に設定。。
（実験では吐出速度0）乱流モデル：LES (OpenFOAM v2.1.1)
SGSモデル：1方程式モデル
壁関数：Spalding則
時間項：２次精度の後退差分
対流項：線形補間（２次精度）（濃度の対流項のみTVDスキームを使用）
床面と円筒はno-slip、側面と天井はslip条件

円筒レイノルズ数が臨界レイノルズ数（2x105）を超えると、剥離域は小さくなり、煙突ダウンウォッシュは緩和される。
排ガスに浮力が入るとフルード数の相似が必要。ただし、排ガスの運動量の相似も満たさなければいけないので、密度を変化させることができず、実験風速をかなり遅くする必要がある。

10 m/s32D
10D

煙突径：D=30mm
煙突高さ：Hs=180mm
水平風速U：10m/s
排ガス吐出速度Ws：0 m/s or 10m/s
吐出速度比Ws/Us：１
レイノルズ数（UD/ν）：2 x 104
排ガスの浮力なし

17格子タイプ

18格子タイプ→ 基準ケース
セル数は1600万
O-type Coarse Mesh (OC)

デカルト格子 (CM)O-type Fine Mesh (OF)

非構造格子 (UM)openfoam.org

19

層状格子 (LM)20y+圧力分布（sgsモデル）

CM~40

UM~40

LM~15

Wind

2~40

Smag (Cs=0.05)
Smag (Cs=0.1)
Smag (Cs=0.15)
OneEq
Dynamic Smag
Park and Lee (2002)

1OC

格子タイプ：OF z/Hs=0.806−2OFのみ粘性底層まで解像。

0~3Cp

OF

−1y+

圧力係数

格子タイプ

円筒前面

円筒後背部

円筒側部

−3yPlus.boundaryField()[patchi] = d[patchi]*sqrt(nuEff.boundaryField()[patchi]*mag(U.boundaryField()[patchi].snGrad()))/ nuLam.boundaryField()[patchi]
dは壁面からの距離。

050

100150

Angle (deg.)→LESでvanDriest型の壁面減衰関数を使用する場合、毎ステップy+を計算し、Δを修正する
ので、時間が非常にかかる。
→LESPropertiesでcalcIntervalを設定することが可能。
2122

圧力分布

圧力分布（格子タイプ）
z/Hs=0.806

UM2

OF

UMのケースはほとんど剥離していない。
→OCの風上補間を使ったケースに類似。
OCLM

CM

OC 風上補間全ての格子タイプで剥離領域の大きさはOFと異なる。

−3−2−1

Cp01

Park and Lee (2002)
OF W/U=0
OC W/U=0
CM W/U=0
UM W/U=0
LM W/U=0050100150Angle (deg.)

OF以外は円筒側面で低圧（高風速）に。
特にUMは、円筒後背部での圧力回復が大きい。
剥離位置も格子タイプによって異なる。
23z/Hs=0.824

吐出速度の影響

圧力分布（吐出速度）y2z/Hs=0.806

煙突中心断面の流線

x−2−1Cp01W/U=0W/U=1
Park and Lee (2002)

吐出速度なし−3OF
050100150Angle (deg.)

Tip vorticesとJet wakeの相互作用？

吐出速度ありz/Hs=0.82526

濃度分布

OC

1.0

1.5

OF

0.5

∆H/D

2.0

OF W/U=1
OC W/U=1
CM W/U=1
UM W/U=1
LM W/U=10.0

プルーム上昇高さ

2.5煙軸上昇軌跡

05101520x/D

地表濃度分布は高解像度格子の方が高くなる。

風下距離

煙軸上昇軌跡はUMのみ大きく異なる。
2728

解像度の影響

剥離域の過小評価の要因（１）
tip vorticesの傾き？

OF

1. 壁面せん断応力の過小評価
0.03

OC0.01−0.010.00

τw/U2ref0.02OF
OCCMUMLM

050

100150Angle (deg.)

壁関数による制約（壁面近傍の速度勾配を正確に再現できない）
2930

剥離域の過小評価の要因（２）

剥離域の過小評価の要因（３）3. 数値粘性

2. 数値振動
Turbulence transportLaminar

OFOC-upwind<u’v’>の分布OC<u’v’>の分布
数値振動による擾乱により運動量の乱流輸送？

数値粘性により壁面の速度勾配が過小評価。
せん断応力が小さくなり剥離が起こりにくくなる。
3132

改善策

風洞実験との比較

1. 壁面応力モデルを使用する。

境界層中の拡散煙突内径：6mm煙突外径：8mm
煙突高さ：60mm
吐出速度：1.43m/s
煙突高さ水平風速：1.43m/s

1m

2. 高次の離散化スキームを使用する。
→非構造格子には難しい。

10m
1m

3. レイノルズ数に応じて所定の位置で強制的に剥離させる。

3334

気流場（平地実験との比較）

地表濃度

地表濃度

煙源位置において
OCUM
CM
W/U=1<u>/U

σu/U

σw/U

OFは計算コストの観点から適用できず。
UMの場合最大着地濃度が30%超の過小評価<u’w’>/U23536

その他の煙源の設定

煙源の設定の影響
乱れあり

Case1：基準ケース
Case2：排ガスの乱れあり（乱流強度10%）
Case3：外径考慮せず
Case4：煙突内の流れ再現

基準ケース
基準ケース

煙突内部流

煙突内部流一定流速（基準ケース）

外径＝内径37まとめ
１．煙突周囲の格子タイプや格子解像度により、剥離位置や後流の大きさは大きく変化する。
剥離位置に影響を与える要因は、粗い格子による壁面せん断応力の過小評価、数値粘性による運動量の拡散、数値振動による運動量の乱流輸送である。
２．煙突後方の剥離を正確に再現できなければ、地表濃度を精度よく求めることができない。対策として、壁面応力モデルなど高精度の壁面モデルの利用が考えられる。
３．煙源の設定によって，最大着地濃度に10％～20％程度の差が生じ得る。ただし、濃度を過大評価する要因と過小評価する要因があり、それらが打ち消し合われている。

39

ケース

基準ケースとの比（最大着地濃度）

Case11.00

Case21.00

Case30.93

Case40.83

煙突直上の風速分布38