WindPerfect 建築ＣＦＤ 概説
２０１７年５月
株式会社環境シミュレーション

１．WindPerfect における気流解析シミュレーション

ここでは、（株）環境シミュレーションの開発した３次元熱流体解析シミュレーションプログラ
ム WindPerfect の利用に関して建築ＣＦＤの理論を含めた概要について述べる。

一般的な風環境解析での CFD（計算流体力学）の概要

解析ソフト

WindPerfect2015

乱流モデル

標準 k-εモデル, LES（標準 Smagorinsky モデル）

解析方法

擬定常解析（TimeMarching 法）

移流項スキ-ム

ハイブリッド中心差分法

流入境界条件

別途記載

流出境界条件

自由流出入または規定流速

乱流擾乱：25%

地表面,建物境界条件

対数則

側面,上空境界条件

Free-slip

初期値

べき乗則による鉛直風速プロファイルを全域に設定

1２．基礎方程式

流体解析で取り扱う流体は、通常，水，油，溶融金属などの非圧縮粘性流体です。 また、空気，水蒸気，ガスなどもマッハ数（音速との速度比）が０．３以下ならば非圧縮性流体として取り扱っても問題無いとされています。 このような流体の挙動は、次のような物理方程式で記述されます。

◆運動方程式（Navier-Stokes 方程式）：∂u/∂ｔ ＋ ∇(uu) ＝ －∇(P/ρ) ＋ (μ/ρ)∇２u ＋ Ｆ・・ (1)◆連続の式：∇･u ＝0・・ (2)◆エネルギー方程式：∂E/∂ｔ ＋ ∇(Eu) ＝ －(k/Cρ)∇２E ＋ q

ここでは、
u：速度ベクトル，P：圧力，ρ：流体の密度， μ：粘性係数，
Ｆ：体積力， E：エネルギー密度（温度)， k：熱伝導率，
C：比熱， q：発熱量 ∂/∂ｔ：時間に対する偏微分，∇：nabra（∂/∂x+∂/∂Y+∂/∂z）と、しています。

2・・ (3)

(1) 式で示される非定常 Navier-Stokes 方程式は、流体の移流項を含めた挙動が du/dt＝∂u/∂ｔ＋∇(uu)なる実質微分(Substantial Derivative) で表され、それが圧力の勾配と速度の発散に比例する粘性力に比例する事を示しています。 数学的な性質としては双曲型方程式と楕円方方程式の混合形であり、これが双曲型方程式と法物方方程式の混合形である
非定常圧縮性流体と分けて論議される所以です。 (2) に示される連続の式は質量保存側を表現しており、Navier-Stokes 式と連成して解かれます。 (3)式のエネルギー方程式は、直接これら方程式の解法に影響する訳ではなく、例えば浮力項（(1) のＦ）を介して温度分布の影響が流れ場の評価に反映されます。

3

３．数値解析方法

ここでは、前項で得られた基礎方程式に対して、数値的な解を求める方法と計算のアルゴリズムを述べる。
微分方程式の離散化手法については有限差分法、有限体積法（コントロールボリューム法）、有限要素法などいくつかの種類があるが、ここでは非定常の流体計算において適用例が多い有限差分法を用いる。
まず、前項の式（1）の左辺を速度の時間変化項のみにする。

ui u juitxj

1 PxixjSGSuixjuj（1）xi

ここで、式（1）の右辺を u と P の関数 f u, P とすると

uitnnf ui , P n（2）これを離散化（時間に関して前進差分化）する。その際、圧力と速度に関して陽的（explicit）に取り扱う。

uin 1uinnt f ui , P n（3）

式（3）は full-explicit である。

また、速度場の積分に関しては、SMAC（Simplified Marker And Cell）法を採用する。この方法は速度場を陽的に解いた後、得られた中間速度を予測子として、次の時間ステップで連続の式を満たすように圧力場を求めようとするものである。
なお、移流項の積分については、ハイブリッド中心差分法を用いた。その他の空間差分につ い て は 、 1 次 精 度 の 中 心 差 分 を 採 用 し た 。 計 算 時 間 間 隔 は 、CFL(Courant-Friedrichs1-Lewy)条件の 80％である。（クーラン数 C4utx0 .8 ）

上述の数値解析フロー図を図 1 に示す。

①渦粘性係数の算出
n

速度分布 u をもとに前述のモデルより算出

n② f u ,Pnの移流項積分値算出
n

速度分布 u をもとにハイブリッド中心差分法より算出*③速度場の予測子 u を算出
第１ステップでは、式（22）より陽的に算出④正規化圧力2n 1、un 1）を算出

u * 、 s を修正

u*s⑤uPss（n 1を算出
su*s、n 1snsst

速度の残差が誤差範囲に収束

ＹＥＳ
次のタイムステップ

図 1 数値計算フロー図5ＮＯ４．アルゴリズム

前項で示した方程式の離散化は次のような手順で行われます。
まず、(a)式の左辺を速度の時間変化項のみにします。

∂u/∂ｔ ＝ －∇(uu)－ ∇(P/ρ) ＋ (μ/ρ)∇２u ＋ Ｆ・・(1)

ここで、(1)式の右辺を u の関数なるｆ(u)とすると、∂u/∂ｔ ＝ ｆ(u)・・(2)

これを離散化（時間に関して前進差分化）すると次の式を得ます。

un+1 ＝ un ＋ △ｔ・ｆｎ・・(3)ここで、t＝n+△ｔ，ｆｎ＝ｆ(un, Pn) であり、この差分化は Euler 法と呼ばれます。 これは時間に対して１次精度の差分であり、また陽解法(Explicit Method)であるため、時間刻み△ｔを CFL条件を満たすように十分小さくとる必要があります。
このように離散化した Navier-Stokes 式を、連続の式と直接連成させて解く事は非常に困難です。 これは Navier-Stokes 式の非線形性が強いためと、連続の式との組み合わせではNavier-Stokes が剛な方程式になり、安定な収束が得られにくいためです。 このため、Navier-Stokes 式と連続の式を直接連成させるのではなく、弱い連成条件の下に解く手法が一6

般的となっています。 その中の１つがここで採用する SMAC 法(Simplified Marker And CellMethod)です。 この方法は速度場を陽的に解いた後、得られた中間速度を予測子として、次の時間ステップで連続の式を満たすように圧力場を求めようとするものです。
SMAC 法の手順は次のようである。

・u' ＝ un ＋ △t・ｆ(un, Pn)
により速度場の予測子 u' を求めます。

・∇2φs ＝ ー∇・u'よりスカラポテンシャルφs（＝Ps／ρ）を求めます。
（∂φs/∂n＝0 on Neumann 境界，φs＝0 on Diricret 境界）

・un+1 ＝ u' ＋ ∇2φs ， φn+1 ＝ φn － φs／△tを用いて un+1，φn+1 を求める。 ・～・を繰り返します。

7

5．乱流の取り扱い

風の現象は、いわゆる乱流と呼ばれる空気の流れであり、微少地形の影響を受けた広域にわたる風の場を再現するためには、DNS（Direct Numerical Simulation）と呼ばれる手法によって、連続の式と運動方程式を直接数値的に解き、流れに含まれる渦を完全に計算することが望ましい。 しかし、高レイノルズ数の乱流現象を解析するには多くの計算格子が必要となり、現状ではコンピュータの性能からいっても困難なため、強い乱流の流れの場を対象に計算を行うには乱れをモデル化して取り扱う必要がある。 これまで提案されている代表的な乱流モデルとしては、図 2 に示すように、平均操作を時間について行うものと空間的行うものに大別できる。

Direct NumericalSimulation

フィルタリング（空間平均）

乱流数値シミュレーション
N-S方程式

LES

非線形項（2次モーメント）

ＤＳＭ（応力方程式モデル）

ui u j2-方程式
モデル
レイノルズ平均（時間平均）
1-方程式
モデル

0-方程式
モデル
渦粘性型モデル
ＲＡＮＳモデル図 2 乱流モデルの分類

ここでは、普遍性が期待できるミクロなスケール（格子スケール）以下の渦だけをモデル化し、流れの場の影響を強くうけるマクロなスケール以上の渦についてはモデルに頼らず直接計算することで、DNS より粗い格子でも理論的破綻が少ない LES を適用するのが最近の趨勢であるし、実用的には標準 k-εモデルのような RANS（レイノルズ平均ナビエストークス）モデルが用いられる事も多い、

8６．LES 乱流モデル

非圧縮性流体の流れの場を支配する基礎方程式は、以下に示す質量保存の連続式と運動量保存の Navier-Stokes の式である。

連続の式

uixi0（1）Ｎavier-Stokes 方程式（等温状態）

uit

ui u jxj

1 Pxixj

ここで、 u i は各軸方向速度、 x i は座標位置、
：

uixj

uj

（2）

xiは密度、 P は圧力、 は分子動粘性係数（＝は粘性係数）である。

連続の式は流体中における質量保存則を表す。Navier-Stokes 方程式は非圧縮性流体の運動方程式であり、左辺第 1 項を時間変化項、第 2 項を移流（対流）項、右辺第１項を圧力項、第 2 項を粘性項と呼ぶ。速度の時間変化は圧力と粘性力に依存するというのがこの式の物理的意味であり、左辺第 2 項の移流項は流れ場を生成する項で流体の非線形性を表す。
マクロスケールの運動は、一般的にミクロスケールの運動よりもはるかに多くのエネルギーを保有しており、その大きさや強さは流体中の保存量の輸送に大きく影響する。一方、ミクロスケールの運動は通常、マクロスケールの運動に比べてはるかに弱く、物理量の輸送にはほとんど寄与しない。そこで、マクロスケールの渦をより正確に取り扱う手法である LES（Large EddySimulation）を採用し、全ての流れの運動の中で、マクロスケールの成分のみを含む速度場の解を計算対象とする。後述するように、流れの基礎方程式にフィルタリングを施した場合にも、SGS（Sub Grid Scale）の新たな未知数が現れ、方程式系がクローズしなくなる。この未知数を何らかの方法でモデル化して与える方法が LES のモデリングであり、この際に用いるモデルのことを SGS モデルと呼ぶ。
フィルタを施した速度場は次のように定義される。

ui x

G x, x u i xdx（3）

ここで、フィルタ関数 G x, x は局所的な関数とされ、ガウシアンフィルタ、矩形フィルタ、カットオフフィルタが多用される。いずれのフィルタも長さスケール に対して定義され、よりも大きい渦が直接シミュレーションされるマクロスケール成分であり、 よりも小さい渦がミ9クロスケール成分でありモデル化を必要とする。
式（1）と式（2）に示す連続式、Navier-Stokes 方程式に空間フィルタリングを施せば、フィルタリングされた連続式と運動方程式が得られる。

uixi0（4）ui u juitxj

ui u jij

1 Pxiuixjxj

uj

ij

xixi

ui u j（5）（6）

式（5）の右辺第 3 項の ij はフィルタリングにより新たに現れる未知数であり、SGS レイノルズ応力と呼ばれる。SGS の乱れによる流れ場の影響は ij を通じてフィルタリングされた運動方程式に組み込まれる。式（6）の右辺第 1 項の u i u j を u ju j uiu i （ ui ：格子平均からのずれ）と表現して分解すれば、 ij は次式のように表される。

ij

Lij

Cij

Rij

（7）

Lijui u jui u j

（8）

C ij

ui u jui u j

（9）

Rij

ui u j（10）

Lij 、 C ij および Rij はそれぞれ、Leonard 項、Cross 項および Reynolds 項と呼ばれる SGS 項である。これらの量を u i などの GS の量に結びつけ、式（4）、式（5）を閉じた方程式系にするための工夫が LES におけるモデリング（SGS モデリング）である。

SGS モデルとして最も歴史が古く、かつ一般的に使われているのは、スマゴリンスキー（Smagorinsky）が提案したスマゴリンスキーモデルである。スマゴリンスキーモデルでは、 ij を構成する Leonard 項 Lij 、Cross 項 C ij 、Reynolds 項 Rij のうち Lij と C ij はお互いにうち消し合う効果があり、両者を合わせると無視できる（ LijC ij

100 ）と考え、 ij は Rij のみで代表されると仮定する。 Rij は SGS の渦粘性係数

SGS を導入し、フィルタリングされた運動方程式の変形速度テンソル（歪み率） S ij に比例するものとして以下のようにモデル化する。

ijS ijRijui u j

12uixj

SGSui
xj

uj
xi

1
ij Rkk3

uj

2 SGS Sij

1ij Rkk
3（12）xi

11（11）

SGS 渦粘性係数は次元解析より、次のように導かれる。

SGSCs22S（13）1SS ij S ij 2（14）C s は Smagorinsky 定数と呼ばれ、スマゴリンスキーモデルで唯一のモデル定数であり、一般的には 0.10 である。また、 はフィルタ幅のスケールで、一般に13123（i；i 方向の計算格子幅）である。
式（13）を用いて、式（11）により ij を近似するのが、スマゴリンスキーモデルである。なお、式（11）に現れる Rij の対角成分 Rkk は、

SGS とを用いて次元解析により代数的に与えることもできるが、通常は、式（5）の圧力項に含ませ、陽的に扱わないことが多く、ここでもその扱いを踏襲する。
また、乱流粘性は分子粘性に対して非常に大きい事を考慮すると、式（5）と式（11）から、次の基礎方程式が得られる。

uitui u jxj

ここでは、式（15）の1 PxixjSGSujuixjxi（15）

SGS を壁面近傍以外では上述の式（13）より求め、壁面近傍の箇所については、下記の混合距離モデルにより算出している。
Prandtl の混合距離モデルでは、まず、乱流の特徴的な速度スケール vt とその特徴的な長さスケール l を用いて次元解析により、tvt lで表せるものとする。ここで、（16）t は渦粘性係数（＝t：t は乱流粘性係数）である。次に、特徴的速度 vt は乱流の特徴的的時間スケール t 0 を用いて vtのように表現する。

12l t 0 と定義し、式（16）を以下

l2t0t（17）さらに、壁乱流の性状を考慮して、 t 0 が、 1ux n （ u ：壁面接線方向速度、 x n ：壁面鉛直方向座標）に比例するものとすると、式（17）は次式となる。

l m2tu（18）xn

式（18）中の l m は混合距離である。
壁面近傍で、LES で求めた渦粘性係数（式（13））と混合距離モデルで求めた渦粘性係数（式（18））が対応する必要があり、次式が成り立つ。

Cs22S

l m2u（19）xn

式（19）を満たすように l m を求め、式（15）の渦粘性の効果を評価する。

13７．ｋ－ε乱流モデルの理論

屋外環境の予測については、風、日射、人口排熱（冷房・自動車）など、様々な要素が関連している。 したがって、このような屋外温熱環境を総合的に解析するためには、対流と熱・湿気輸送の連成解析のための基礎方程式が必要である。

① 連続の式

uixi

0② 運動方程式

uiuiujt

1xj

Pxi

xitkxjPk

Gkxjk

uiujxj

xi

C3G k

C2gi

0

③ 乱流エネルギーｋの方程式

ktujkxjtxjk④ 粘性消散率εujttxjxj

C1 Pk⑤ 温度θの方程式

tujxjXjuj⑥ 絶対温度ｑの方程式

qtGkujg3qxjxj

u jqu3の基礎方程式が必要である。これらを連成して解析することで、屋外温熱環境の形成に係る様々な物理現象が明らかとなる。
流体とは、そもそも圧縮性であり、非圧縮性の流体は、圧縮性流体の近似である。非圧縮性
近似が成り立つのは、流速 u が音速ｃと比較して十分に遅く、音速を無限大と近似してもかまわない場合、すなわち、マッハ数（M=u/c）が１よりも十分小さい場合である。したがって、温度14変化などで流体の密度が変化する場合でも流速 u が音速ｃと比較して十分に遅い場合は、非圧縮性の流体と扱っても良いと考えられる。

■非圧縮性流れの支配方程式

建築物周辺の気流を扱う場合、マッハ数が 1 よりも十分に小さいときには、非圧縮が良い近似となる。非圧縮流れでは密度が一定であると仮定する。こうすると、連続の式およびナビエ・ストークス方程式は次のようになる。

uixituit0u j uixj

Pxixj

uixjujxiここで、νは動粘性係数で次のように仮定される。

非圧縮流れでも、移流項は双曲形、粘性項は放物形である。連続の式は、流量の発散が 0 で表される。これを時間で微分し、次にナビエ・ストークス方程式の圧力項と時間微分項だけを抜き出した式に代入すると、
2P
xi xi0が得られ、これは、圧力に関する楕円形方程式である。

■ｋ－ε乱流モデルの適用

乱流には、微細な変動から大きな変動まで、様々なスケールの乱れが含まれている。乱流を数値的に解析する場合、直接数値シミュレーション（ＤＮＳ）と呼ばれる手法を用いて「連続の式」「運動方程式」の基礎方程式を解き、流れに含まれるすべてのスケール渦をシミュレートすることが理想的な方法である。しかし、ＤＮＳの場合、すべての渦をシミュレートするのは、十分15に細かい計算格子が必要であったり、計算機能力の制約上現在では不可能であると考える。
そのため、ある程度以上のスケールの渦のみをシミュレートする対象として、それより小さい渦に関しては、モデルにより表現するという手法が用いられている。 これが乱流モデルである。

速度と圧力の瞬時値を平均値と変化量に分ける。

ui

UiuipPpu i : 瞬時値

U i : u iのアンサｈンブル平均

u i : u iの変化量

p : 瞬時値

P : pのアンサンブル平均p : pの変化量

これを連続の式に代入し、アンサンブル平均をとると〈ui’〉=0 を用いて

Ui
xi

ui
xi

ui
xi

ui
xi

Ui

ui

xixi

Uixi

Ui
xi00が得られる。これは、平均流も変化量も連続の式が満足することを示している。
同様に、Ｎ-Ｓ方程式にも代入し、アンサンブル平均をとる。

t

Ui

Uit

uixi

U iU jxj

U iU jPxi

U iu j

U juiui u j
Ujxi

Uixj

xixiPpxjxj

UiuixiUjuj

ui u jが得られる。これは、レイノルズ方程式と呼ばれ u i u j が、変化量の影響がせん断応力を密度で割ったものの形式に変形されることを示しており、レイノルズ応力と呼ばれる。
この方程式は、u i u j が不明であるので解くことができない。そこで、これをモデル化して計算をする。

16

■標準 k-εモデル

u i u j は乱れによって発生するせん断応力である。このせん断応力が平均流の勾配に比例すると仮定するものが渦粘性モデル（EVM）である。すなわち、等方な渦動粘性係数νt を仮定して

ui u jijt

Uixj

Uj2k ij3xi

はクロネッカーのデルタij1(ij ), ij0(ij)と近似する。ここで、k は乱流エネルギーで、次式で定義される。

1ui u j2k12u12u22u32

渦動粘性係数νt は比例定数の形式をとっているが定数ではない。流れの状態により大きく変化する。Νt は渦拡散係数とも呼ばれ、乱れによって行われる拡散を代表する量である。
この点を考慮し、νt は、

tk 1/ 2lとモデル化する。K1/2 は、乱れの平均速度という意味を持ち、ｌは乱れの影響が到達する距離という意味を持つ。
次に粘性消散率εを導入する。εは、分子粘性により乱流エネルギーk が熱エネルギーに変換される割合であり、十分に発達した乱流では次式で評価される。

CD kl k1/ 2CD k 3 / 2l

CD は比例定数である。 l k 1 / 2 は、時間の次元を持っており、乱れの影響が距離ｌまで到達する時間あるいは乱流エネルギーｋが熱エネルギーに変換されるまでの時間と考える。粘性消散率は乱流エネルギーをこの時間で割った程度のものとして評価できる。

17

. k-εモデルの表式は次のようである。 渦粘性係数を算出するに当たり、乱流エネルギーk及び散逸率εの発生と移流拡散を取り扱う。

標準 k-εモデルでの生産項の表式は次のようである。 また各係数の設定を示す。

18

■流入境界条件

流入境界条件に与え乱流エネルギーｋの分布は、例えば次のような風洞実験値を参考に設定することも可能であるが、乱流エネルギーのプロファイルは周辺の地盤や構造物の状況で大きく異なる場合があるので、本解析では便宜上乱流の擾乱を２５％としてｋを設定している。
散逸率εもｋに準じて設定する。

なお、風荷重シミュレーションでは、流入風は実感変動風として設定する必要があり、片岡の方法・小野の方法や合成渦法を適用する。

19

■壁関数

一般に、層流底層から乱流域の主流方向への接線方向風速は壁法則に支配されると考えられる、 本解析では、べき乗則（風工学における鉛直風速プロファイルのべき乗則とは異なる）より比較的精度が良いとされる対数則を壁関数に用い、境界第一セルの風速値を設定し、壁法則を実現している。 壁関数設定の考え方は次図のようである。

境界第一セルの風速値を設定する方法は、対数則から壁面接線応力τw を計算し、次の時間ステップの風速値を推参する方法である。 ただこの方法は発達中の乱流境界層に用いるには問題があり、場合によって適用には注意が必要である。

20