Studies on an Air Environmental Flow Analysis in Urban Area by Stabilized Finite Element Method Based on Large Eddy Simulation

(1)

修士論文要旨（

2008

年度）

LES に基づく安定化有限要素法による都市の大気環境流れ解析に関する研究

Studies on an Air Environmental Flow Analysis in Urban Area by Stabilized Finite Element Method Based on Large Eddy Simulation

土木工学専攻

33

号八田政知

Masatomo HATTA

1. 研究目的

都市域において，良好な大気環境を創出するためには，

都市計画の段階から，ビル風やヒートアイランド現象などの大気環境の予測・検討を行う必要がある．都市域の大気環境問題で取り扱われる流れは，一般にレイノルズ数が非常に大きく，そして，建物周辺では剥離，再付着，

循環などさまざまな流れの性状を含み複雑な乱流場を形成する．このような現象の検証と評価は模型実験のみならず，近年では数値解析手法の進歩により数値解析によっても行われてきている．また，計算機性能の発展に伴い

DNS

による乱流解析も行われるようになってきているが，計算負荷・容量が膨大なものとなり，実用上は適用が困難な場合が多い．このため，

LES

や

RANS

に代表される乱流モデルが提案され実用的に用いられている．その中でも

LES

は，

RANS

などと比較して，計算量は多くなるが，非定常解析を行うため時間変動する現象を把握するには有効な手段である．また，数値解析により都市の大気環境流れ解析を定量的に行うためには，

構造物や地形などを正確に表現することが望まれる．これらの任意複雑形状の適合性において，有限要素法は優れた手法である．

著者らはこれまで，

Smagorinsky

モデル⁽¹⁾ を用いた安定化有限要素法⁽²⁾ による平行平板間乱流の

LES

を行ってきた⁽³⁾．しかし，建物周辺気流などの複雑な流れ場や非等温場への適用は十分ではなかった．そこで本論文は，都市の大気環境流れ解析手法の構築を目的とし，

複雑乱流場および非等温場における本手法の精度検証を行った．数値解析例として立方体周辺気流解析および地表面高温領域周辺気流解析を取り上げ，実験値⁽⁴⁾⁻⁽⁶⁾ および既存の数値解析解⁽⁴⁾^，⁽⁵⁾と比較を行った．

2. 数値解析手法 2.1

基礎方程式

非等温場における非圧縮性粘性流体を考え，

Boussi- nesq

近似を仮定する．そのときフィルタリングおよび無次元化を施した，

Grid Scale

（

GS

）の運動方程式，連続式，エネルギー方程式はそれぞれ式

(1)

，

(2)

，

(3)

で表される．

運動方程式；

∂u

i

∂t + u

j

∂u

i

∂x

_j

+ ∂p

∂x

_i

− 1 Re

∂

∂x

_j

³ ∂u

i

∂x

_j

+ ∂u

j

∂x

_i

´ + ∂τ

ij

∂x

_j

− Arθδ

i2

= 0 (1)

連続式；

∂u

_i

∂x

i

= 0 (2)

エネルギー方程式；

∂ θ ¯

∂t + u

_j

∂ θ ¯

∂x

j

− 1 P rRe

∂

²

θ ¯

∂x

²_j

+ ∂h

_j

∂x

j

= 0 (3)

τ

_ij

= u

_i

u

_j

− u

_i

u

_j

(4)

h

j

= u

j

θ − u

j

θ (5)

ここで，

u

i，

p

，

θ

はそれぞれフィルタリングを施した流速，圧力，温度である．

Re(= U L/ν)

はレイノルズ数，

P r(= ν/α)

はプラントル数，

Ar(= − gβ(θ

w

− θ

c

)L/U )

はアルキメデス数をそれぞれ表している．但し，

U

は代表流速，

L

は代表長さ，

ν

は動粘性係数，

g

は重力加速度，

α

は温度拡散係数，

β

は体膨張係数，

(θ

w

− θ

c

)

は代表温度差をそれぞれ表している．また，

τ

ij は

SubGrid Scale

（

SGS

）応力，

h

j は

SGS

熱流束を表す．格子で捉えきれない

SGS

の乱れによる

GS

の流れ場への影響は，

τ

ijを通じて

GS

の運動方程式に組み込まれる．

h

jはエネルギー方程式にフィルタリングを施すことにより，新たに現れる未知数である．したがって，

LES

を用いて非等温場を解析する場合は

τ

ij，

h

jの

2

つに対してモデル化を行う．

2.2 Smagorinsky

モデル

Smagorinsky

モデル⁽¹⁾は，

SGS

応力

τ

ijに対するモデル化の中で最も代表的なものである．乱流エネルギーの収支において生産項と散逸項が釣り合うという局所平衡状態を仮定すると，

SGS

応力は以下のようにモデル化される．

τ

_ij

= − 2ν

_SGS

S

_ij

(6) ν

SGS

= (C

s

∆)

²

q

2S

ij

S

ij

(7) S

_ij

= 1

2 ³ ∂u

_i

∂x

j

+ ∂u

_j

∂x

i

´ (8)

∆ = V

e

1

3

(9)

(2)

h

j

= − α

SGS

∂θ

∂x

j

= − ν

SGS

P r

SGS

∂θ

∂x

j

(10) ν

SGS は渦粘性係数，

S

ij は

GS

の変形速度テンソル，

C

sは

Smagorinsky

定数

, ∆

はフィルター幅であり

4

面体要素の体積

V

eの

3

乗根で定義する．また，

α

SGSは

SGS

温度拡散係数，

P r

SGSは

SGS

プラントル数である．

Smagorinsky

モデルの場合，式

(10)

の

ν

SGSは式

(7)

から算出し，

SGS

プラントル数

P r

SGSに対しては

0.5

程度の一定値を与える．

2.3

境界条件

Dirichlet

型，

Neumann

型の境界条件は，以下のように与えられる．

u

_i

= g

_i

on Γ

_g

(11)

½

− pδ

_ij

+ 2 ³ 1

Re + ν

_SGS

´ S

_ij

¾

n

_j

= h

_i

on Γ

_h

(12)

θ = t on Γ

t

(13)

³ 1

P rRe + α

SGS

´ ∂θ

∂x

_j

= k on Γ

k

(14)

ここで，

Γ

_g，

Γ

_tは

Dirichlet

境界，

Γ

_h，

Γ

_kは

Neumann

境界を表している．

g

_i，

h

_iはそれぞれ流速，トラクションの既知量，

t

は壁面既知温度，

k

は壁面上の法線方向温度勾配，

n

_jは外向き単位法線ベクトルを表す．

2.4

安定化有限要素法

基礎方程式

(1)

，

(2)

に

SUPG/PSPG

法，式

(3)

に

SUPG

法に基づく安定化有限要素法²⁾を適用すると以下の弱形式が得られる．なお，

w

i，

q

は重み関数である．

運動方程式，連続式；

Z

Ω

w

i

³ ∂u

i

∂t + u

j

∂u

i

∂x

_j

´ dΩ −

Z

Ω

∂w

i

∂x

_i

pdΩ +

Z

Ω

³ 1

Re + ν

SGS

´ ∂w

i

∂x

_j

³ ∂u

i

∂x

_j

+ ∂u

j

∂x

_i

´ dΩ

− Z

Ω

w

i

Arθδ

i2

dΩ Z

Ω

+q ∂u

i

∂x

_i

dΩ +

nel

X

e=1

Z

Ω_e

n τ

_supg1

u

_k

∂w

_i

∂x

k

+ τ

_pspg

∂q

∂x

i

o

· n ∂u

i

∂t + u

j

∂u

i

∂x

j

+ ∂p

∂x

i

− ∂

∂x

j

³ 1

Re + ν

SGS

´³ ∂u

i

∂x

j

+ ∂u

j

∂x

i

´ − Arθδ

i2

o dΩ

= Z

Γ_h

w

i

h

i

dΓ (15)

エネルギー方程式；

Z

Ω

w ³ ∂θ

∂t + u

_j

∂θ

∂x

j

´ dΩ

+ Z

Ω

³ 1

P rRe + α

_SGS

´ ∂w

∂x

j

∂θ

∂x

j

dΩ

+

n_el

X

e=1

Z

Ωe

τ

supg2

u

k

∂w

∂x

k

³ ∂θ

∂t + u

j

∂θ

∂x

j

− ( 1

P rRe + α

_SGS

) ∂

²

θ

∂x

²_j

´ dΩ =

Z

Γ_k

wkdΓ (16)

式

(15)

の左辺第

1

〜

4

項と右辺項は式

(1)

，

(2)

に対する通常の

Galerkin

項，第

4

項は移流項の卓越に起因する数値不安定性を抑制するための安定化項（

SUPG

項）

および非圧縮条件に起因する圧力振動を回避するための安定化項（

PSPG

項）である．右辺項は，

Galerkin

項の粘性項を部分積分したことにより生じた境界積分項である．式

(16)

の左辺第

1,2

項と右辺項は式

(3)

に対する通常の

Galerkin

項，第

4

項は

SUPG

項である．

τ

_supg1，

τ

_pspg，

τ

_supg2は安定化パラメータである．

Ω

は解析領

域，

Ω

_eは要素領域を表している．

そして，式

(15)

，

(16)

に対して

P1/P1

（流速・圧力

1

次）要素を用いて補間を行うことにより，以下の有限要素方程式を得る．

(M + M

_δ

) ∂u

i

∂t + n

K(u

_j

) + K

_δ

(u

_j

) o

u

_i

− (G − G

_δ

)p + 1

Re Su

_i

− (M + M

_δ

)Arθδ

_i2

= 0 (17)

Cu

_i

+ M

_ϵ

∂u

_i

∂t + K

_ϵ

(u

_j

)u

_i

+ G

_ϵ

p − M

_ϵ

Arθδ

_i2

= 0 (18) (M + M

_δ

) ∂θ

∂t + n

K(u

_j

) + K

_δ

(u

_j

) + 1 P rRe S o

θ = 0 (19)

ここで，

M

，

K

，

G

，

S

，

C

はそれぞれ時間，移流，圧力，

粘性，連続の各項の係数行列であり，添字

δ

，

ϵ

はそれぞれ

SUPG

項，

PSPG

項に起因するものを表す．

時間方向の離散化において時間微分項は次式のように表される．

∂u

i

∂t = u

ⁿ⁺¹_i

− u

ⁿ_i

∆t (20)

時間方向の離散化には

2

次精度を有する

Crank − Nikolson

法を用いた．なお，連続式，圧力は陰的に取り扱っている．移流項における移流速度

u

i は

2

次精度

Adams − Bashforth

法により次式のように近似した．

u

_i

= 3

2 u

_iⁿ

− 1

2 u

_iⁿ⁻¹

(21)

連立一次方程式の解法には，

Element-by-Element Bi-

CGSTAB2

法を用いた．

(3)

15.7L 5.2L

10.0L 9.7L 4.35L

Outlet Inlet

y

z x

図

– 1

解析領域

LineB

x y

z 0 Central axis

1.0L 1.0L

1.0L

(z = 0.0L) LineA

図

– 2

立方体周辺拡大図

2.5

壁関数

建物周辺気流のようなレイノルズ数の大きい流れ場では，壁面境界条件に

No-slip

条件を適用すると計算負荷が増大する．そのため，壁面第

1

節点を全領域にわたって粘性底層内に位置させるのではなく，乱流域に位置することを前提として，壁関数を用いることが多い．瞬時値に対して壁関数を適用する場合として，

Werner

らにより提案された

Linear-1/7 power law

型

2

層モデル⁷⁾ を以下に示す．

V

p

τ

w1/2

= y

⁺

(y

⁺

≤ 11.81) (22) V

_p

τ

w1/2

= 8.3(y

⁺

)

^1/7

(y

⁺

>11.81) (23) V

pは壁面第

1

節点の接線方向速度成分の絶対値

(

底面の場合は

V

p

= p

u

²

+ w

²

)

である．

Linear-1/7 power law

型

2

層モデルは，壁面第

1

節点の瞬時速度の値のみで瞬時の壁面せん断応力を求めることができ，またこれに関する計算負荷をほとんど必要としないという利点がある．

3. 数値解析例 3.1

立方体周辺気流解析

図−

1

に解析領域，図−

2

に立方体周辺拡大図，表−

1

に解析ケースを示す．境界条件は，

Case1

，

Case2

に共通して，流入境界は

u = y

^(1/4) ，

v = w = 0.0.

流出境界は自由流出境界．側面および上端面は

Slip

条件とした．なお，立方体の高さ及び立方体の高さでの流入風速で無次元化している．数値解析に用いたメッシュ

表

– 1

解析ケース

底面および立方体壁面の境界条件

Case1 Linear-1/7 power law

型

2

層モデル

Case2 No-slip

0 1

1 1.1 1.2 1.3 1.4 1.5

ታ㛎୯

᧛਄䉌䋨LES) Case1 Case2

䋼u䋾

y

（

a

）

LineA

0 1

0 0.5 1 1.5

ታ㛎୯

᧛਄䉌䋨LES) Case1 Case2

䋼u䋾

y

（

b

）

LineB

図

– 3

流速

u

の時間平均分布

表

– 2

再付着距離

は，

x

，

y

，

z

方向にそれぞれ

63

×

34

×

48

分割し，節点数は

104,630

，要素数は

582,696

，最小メッシュ幅は

4.16

×

10

⁻²

L

の不等分割メッシュを用いた．そして，

レイノルズ数は

84,000

，微小時間増分量は

1.0

×

10

⁻³，

Smagorinsky

定数は村上らの数値解析^(4),(5) を参考に

0.12

とした．図−

2

より，流速

u

の測定地点は立方体中心の底面を原点として，

LineA

は

x = 0.2L

，

y = 1.0L

〜

1.5L

，

z = 0.0L

，

LineB

は

x = 1.0L

，

y = 0.0L

〜

1.5L

，

z = 0.0L

とした．

解析結果として，図−

3

（

a

），（

b

）にそれぞれ

LineA

，

LineB

での流速

u

の時間平均，表−

2

に屋上面および立方体後方の再付着距離の比較を示す．

図−

3

（

a

）より，

Case1

は実験値および村上らの数値解析解より屋上面付近で流速を過大に評価しているが，

定性的に一致していることが確認できる．一方，

Case2

では屋上面での逆流域を広く評価しているために，屋上面付近の流速が負となっている．図−

3

（

b

）より，

Case1

，

2

ともに，立方体の高さ（

y = 1.0L

）付近までは，実験値および村上らの数値解析解と良い一致を示している．

y = 1.0L

より上部においては，

Case1

は良い一致を示し，

Case2

は流速を小さく評価している．次に，屋上面および立方体後方の再付着距離の比較を行

(4)

8.0L 1.25L

2.0L 2.0L

x y

z

A B C D

0.875L 0.875L

1.3L θ = 1.0

L

図

– 4

解析領域

う．表−

2

より，

X

R，

X

F はそれぞれ屋上面における逆流による再付着距離と立方体後方における循環流による再付着距離である．

Case1

は，

X

R，

X

Fともに村上らの数値解析とほぼ同等の値を示している．

Case2

では，

X

R，

X

F ともに実験値に対して再付着距離を大きく評価している．

3.2

地表面高温領域周辺気流解析

図−

4

に解析領域を示す．境界条件は，流入境界は

u = y

^(1/2)，

v = w = θ = 0.0.

流出境界は自由流出境界．側面，上端面は流速は

Slip

条件，温度は断熱．底面は流速は

Linear-1/7 power law

型

2

層モデル，温度は

θ = 1.0

の領域（

Hot Panel

）以外は断熱とした．数値解析に用いたメッシュは，

x

，

y

，

z

方向にそれぞれ

65

×

43

×

31

分割し，節点数は

86,645

，要素数は

483,840

，最小メッシュ幅は

4.44

×

10

⁻¹

L

の不等分割メッシュを用いた．レイノルズ数は

29, 000

，アルキメデス数は

1.21

，プラントル数は

0.71

，

SGS

プラントル数は

0.5

，微小時間増分量は

1.0

×

10

⁻³，

Smagorinsky

定数は

0.10

とした．

解析結果として，図−

5

に

z

軸方向に平均をとった

x − y

断面平均温度分布を示す．

Hot Panel

によって浮力を生じた流体が解析領域後方に移流・拡散している様子が確認できる．また，高温域は底面付近に存在しており，上部は流入温度に近い低温域が存在している．これは，浮力による鉛直流よりも流入流速による水平流が卓越しているため，温度の移流が水平方向に多く行われたためであると考えられる．図−

6

に断面

A

〜断面

D

での平均温度分布を示す．

Hot Panel

直後の断面

A

の底面付近において実験値よりも過小に評価しているが，各断面において実験値と概ね良い一致を示している．

4. 結論と今後の課題

本論文は，都市の大気環境流れ解析手法の構築を目的とし，複雑乱流場および非等温場における

Smagorinsky

モデルを用いた安定化有限要素法の精度検証を行った．

数値解析例として立方体周辺気流解析および地表面高温領域周辺気流解析を取り上げ，以下の結論を得た．

図

– 5 x − y

断面平均温度分布

0 0.1 0.2

0 0.2 0.4 0.6 0.8

1 ᢿ㕙B ᢿ㕙C ᢿ㕙D

0 0.1 0 0.1 0 0.1

y

ᐔဋ᷷ᐲ ታ㛎୯

ᧄ⸃ᨆ ᢿ㕙A

図

– 6

各断面における平均温度分布

•

立方体周辺気流解析において，壁関数を用いた解析結果は，流速分布および再付着距離において，

実験値および既存の数値解析解と概ね良い一致を示し本手法の有効性が確認できた．

•

地表面高温領域周辺気流解析において，解析結果は平均温度分布において，実験値と地表面付近で差異はみられるが，概ね良い一致を示し本手法の有効性が確認できた．

今後の課題としては，地表面の熱収支を考慮した解析，流入境界に関する境界条件の検討等が挙げられる．

参考文献

1) J.Smagorinsky

：

General circulation experiments with the primitive equations ; part1 the basic experiments

，

Monthly Weather Review

，

91

，

pp.99-164

，

1963.

2) T.E.Tezduyar

：

Stabilized ﬁnite element formulations for incompressible ﬂow computations

，

Advance in Ap- plied Mechanics

，

28

，

pp.1-44

，

1991.

3)

仲嶋晋一，倉橋哲弘，樫山和男：安定化有限要素法による

LES

，計算工学講演論文集，日本計算工学会，

10

，

pp.201- 204

，

2005.

4)

村上周三，持田灯，林吉彦：立方体周辺の非等方乱流場の再現に関する

k − ϵ

，

ASM

，

LES

と風洞実験の比較，東京大学生産技術研究所所報，

43(1)

，

pp.27-35

，

1991.

5) S.Murakami

，

A.Mochida and R.Ooka : Numerical sim- ulation of ﬂowﬁeld over surface-mounted cube with var- ious second-moment closure models, Proc. of 9th Sym- posium on Turbulent Shear Flows

，

pp.1-6

，

1993.

6)

村上周三，大場正昭：床面に温度差のある成層流の気流性状並びに拡散に関する風洞実験，日本建築学会大会学術講演梗概集，