第 24 回数値流体力学シンポジウム	 
D5-6

オープンソース CFD ツールキット OpenFOAM を用いた市街地風環境予測
CFD Prediction of Pedestrian Wind Environment Around Building
Using Open Source CFD Toolkit OpenFOAM○ 今野 雅, 東大院工, 東京都文京区本郷 7-3-1, imano@arch.t.u-tokyo.ac.jp
大嶋 拓也, 新潟大工, 新潟市西区五十嵐２の町 8050, oshima@eng.niigata-u.ac.jp
挾間 貴雅, 鹿島技研, 調布市飛田給 2-19-1, takamasa@kajima.com
柴田貴裕, MTS 雪氷研, 千代田区神田須田町 1-5, shibata@mtsnow.co.jp
小縣 信也, 森村設計, 東京都目黒区中目黒 1-8-8, ogata@ptmtokyo.co.jp
Masashi Imano, Univ. of Tokyo, 7-3-1 Hongo Bunkyo-ku Tokyo, Japan
Takuya Oshima, Niigata Univ., 8050 Ikarashi-Ninocho, Nishi-ku, Niigata, Japan
Takamasa Hasama, Kajima Corp., 19-1, Tobitakyu 2-chome, Chofu-shi, Tokyo, Japan
Takahiro Shibata, MTS Institute Inc, 1-5 Kandasuda-cho, Chiyoda-ku, Tokyo, Japan
Shinya Ogata, P.T.M LTD, 8-8 Nakameguro 1-chome, Meguro-ku, Tokyo, Japan
Various benchmark tests and guidebook for CFD prediction of pedestrian wind environment around building havebeen published by the working group of the architectural institute of Japan (AIJ). In this paper six of thesebenchmark calculation are carried out using Open Source CFD Toolkit OpenFOAM. In addition to that crosscomparisons among the present CFD result, CFD results of other commercial codes and measurement results areconducted. As a result it has been revealed that OpenFOAM has almost same capability of accuracy in prediction ascompared with other CFD codes. Furthermore we have released OpenFOAM data files for these test cases at WEB siteof the Open CAE Society of Japan.
１．はじめに	 	 	  日本建築学会の流体数値計算による風環境評価ガイドライン作
成 WG が作成した風環境の CFD 解析検証用データベースでは、6ケースに及ぶ様々な状況での建物周りの風環境の実測と各種商用・自作コードによる計算結果の比較を網羅的に行っており、その結果は文献(1)等にまとめているほか、WG の WEB ページ(2)において各種データが公開されている。本報ではオープンソース CFDツールキット OpenFOAM を用いて上記 6 ケースの計算を行い、実測や既往の CFD 解析結果との比較を行なった。また、OpenFOAM の計算入力ファイルをまとめ、オープン CAE 学会のWEB ページ(3)上で公開したので報告する。

大評価は他のソフトウェアと比べて同程度と言える．
Fig.2-5 に，角柱後流域の地表付近での過小評価を見るために、x/b=2.0 での<u>の鉛直分布を示す．Fig.2-4 と同様のソフトウェアを用いた計算結果も合わせて示す．建物高さより上（z/b>2.0）では，OpenFOAM は他のソフトウェア同様，計算結果は実験結果と一致しているが，建物高さより下（z/b<2.0）では，KE1 のと近い結果ではあるが，OpenFOAM の結果は風洞実験と比べて最も過小に評価する結果となっている．これは主に移流項の離散化スキームの差が原因と思われる。
Table 2-1	  Calculation condition
流入境界条件ベンチマーク指定値(1)を使用
上面・側面境界条件 Slip 壁流出境界条件
勾配ゼロ
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル
移流項スキーム
TVD（Limited linear）スキーム
計算アルゴリズムSIMPLE，定常解法
計算格子分割法
ベンチマーク指定の格子分割(1)を使用

２．2:2:1 角柱周辺流れを対象とした例（Case	 A	 ）	 ２．１	  計算概要	 	  本章では 2:2:1 角柱周辺流れを対象としたベンチマークケース（Case A）に基づいて作成・公開した計算ケースの結果を示す．
解析領域を Fig.2-1 に，計算条件を Table 2-1 に示す．
移流項に TVD（Limited linear）スキームを，地表面境界に滑面用一般化対数則を適用している以外は学会指定のベンチマーク条件に準じている．
２．２	  計算結果	 
Fig.2-2，2-3 に，主流方向角柱中心線断面における風速の主流方向成分<u>及び乱流エネルギーk の鉛直プロファイルを示す．<u>については、文献(1)の他のソフトでも見られるものがあるが，角柱後流域の地表付近での過小評価が OpenFOAM の結果では出ている．k については、k-εの計算結果の特徴である角柱風上コーナー部分の過大評価が、OpenFOAM でも見られる．
Fig.2-4 に，k の角柱風上コーナー部分の過大評価を見るために、x/b=-0.25 での角柱上面上の k の鉛直分布を示す．文献(1)の他のソフトウェアでの計算結果も合わせて示す（ただし、標準 k-εのみに限定）．他のソフトウェアの計算条件については文献(1)を一部抜粋し，Table 2-2 に再掲する．他のソフトでの計算結果から見て、OpenFOAM の結果は平均的で、k の角柱風上コーナー部分での過Fig.2-1	  Calculation area1Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム	 
D5-6	 
３．4:4:1 角柱周辺流れを対象とした例（Case	 B）	 ３．１	  計算概要	 	  本章では 4:4:1 角柱周辺流れを対象としたベンチマークケース(Case B) に基づいて作成・公開した計算ケースの結果を示す．解析領域を Fig. 3-1 に，計算条件を Table 3-1 に示す．乱流モデルは標準 k-εモデル，Shih モデル，RNG モデル，Realizable モデルの4 ケースを検討した．地表面および角柱壁面に滑面用一般化対数則を適用している以外はメッシュ分割も含め，学会指定のベンチマーク条件に準じている．なお本計算は RANS による定常計算であるため，計算結果の対称性を仮定して主流方向角柱中心線断面で対称境界条件を課し，領域半分のみの解析を行った．計算時間
は Nehalem Xeon プロセッサ (2.93 GHz) 4 Cores を用いて各ケース11 分であった．
Fig. 2-2 Vertical profiles of <u> at the center section of the block.
３．２	  計算結果	 
Table 3-2 に，角柱後方壁面からの再付着距離 XF を示す．標準k-εモデルの結果は実験値より大きく評価されているものの，文献(1)の各種コードにおける同モデルによる結果 0.46 m~0.58 m の範囲内であり，コード間の比較としては妥当な結果といえる．ただし，Shih モデル及び RNG モデルは文献(1)の 0.66 m~0.68 m と比較してやや大きめであり，Realizable モデルは文献(1)の 0.62 m，さらには実験値より小さくなっている．
Fig.3-1，3-2 に，標準 k-εモデルでの主流方向角柱中心線断面における<u>及び k の鉛直プロファイルを示す．<u>は，文献(1)の結果と比較しほぼ同様である．k は角柱近傍部で実験値より過大評価である点は類似しているが，角柱後流域では文献(1)の結果よりやや小さめの傾向がある．
Fig.3-3 に，z = 0.0125 m 面上の各測定点 (位置は文献(1)参照) における風速比を示す．標準 k-εモデルは文献(1)とほぼ同様である．
Shih モデル及び RNG モデルは，文献(1)に見られる領域 3 (角柱後流域) で風速比がそれぞれ最大0.8 及び1.0 程度に過大評価される傾向がより顕著で，それぞれ 1.0 および 1.2 程度となっている．一方，Realizable モデルでは後流域で実験値により近い値が得られているものの，領域 5 (後流側方) では過大評価である．

Fig. 2-3 Vertical profiles of k at the center section of the block.
Table 2-2	  Calculation condition of other case
研究機関

ソフトウェア

移流項スキーム

計算アルゴリズム

ケース名A

STREAM ver2.10

QUICK

SIMPLE，定常解法

KE1B

STAR-LT ver2.0

QUICKSIMPLE，定常解法

KE2

FLUENT ver5.0

中心差分

SIMPLE，定常解法

KE3

自作

QUICK

HSMAC，非定常解法

KE4D

Table 3-1	  Calculation conditions.
流入境界条件ベンチマーク指定値(1)を使用
上面境界条件slip 壁側面境界条件slip 壁。ただし対称面は対称境界条件
流出境界条件
勾配ゼロ
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル，Shih モデル，RNG k-εモデル，Realizable k-εモデル
移流項スキーム QUICK
計算格子分割
ベンチマーク指定の格子分割(1)を使用

Fig. 2-4 Vertical profiles of k above the block.

Fig. 3-1 Calculation area.
Fig. 2-5 Vertical profiles of <u> behind the block.
2Copyright © 2010 by JSFM

Table 3-2	  Reattachment distances.
モデル 標準 k-ε
ShihRNG
Realizable
XF [m]
0.55
0.78
0.730.31

第 24 回数値流体力学シンポジウム	 
D5-6
Fig.4-1 に示す．移流項に TVD（Limited linear）スキームを、地表面境界に滑面用一般化対数則を適用している以外は学会指定のベンチマーク条件に準じている．計算格子は、OpenFOAM 付属の自動格子生成ユーティリティと自作の STL データを使用して生成した．
計算モデルを Fig.4-2 に示す．単純建物群のうち中心の建築物 1個が対象建築物であり、残り 8 個が周辺建築物である．周辺建築物は 8 個全てを各辺 D(=0.2m)の立方体とし、対象建築物の形状は「D D D」「2D D D」の 2 ケースとした．
風向は各ケース 0 、22.5 、45 の 3 通りとした．各計測点は Fig.4-3 に示す通りである．計算時間は IntelCore2Duo プロセッサ (2.99 GHz) 1 Core を用
いて各ケース約 1 時間 30 分であった．

実験0.38

Table 4-1	  Calculation condition
ベンチマーク指定値(1)を使用
流入境界条件
風向,22.5 及び 45 の Case は、2 面を流入条件とした.
上面・側面境界条件 slip 壁勾配ゼロ. 風向,22.5 及び 45 の Case
流出境界条件は、2 面を流出条件とした.
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル
移流項スキーム
TVD（Limited linear）スキーム

Fig. 3-1 Vertical profiles of <u> at the center section of the block.

計算格子分割法

ベンチマーク指定の格子分割(1)を使用(約 8 万要素)

Fig. 3-2 Vertical profiles of k at the center section of the block.

Fig.4-1	  Calculation area(D D D)
(2D D D)Fig.4-2	  Building Model

Fig. 3-3 Comparison of normalized velocity between models.	 
４．単純建物群モデルを対象とした例（	 Case	 C	 ）	 ４．１	  計算概要	 	  本章では単純建物群モデル（Case C）を対象としたベンチマー
クテストの結果を示す．計算条件を Table 4-1 に、計算領域を

Fig.4-3	  Measurement Points
3

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム	 
D5-6
５．街区内に建つ高層建物モデルを対象とした例（Case	 Ｄ	 ）	 ５．１	  計算概要	 	  本章では街区内に建つ高層建物モデルを対象としたベンチマー
クケース（Case D）に基づいて作成・公開した計算ケースの結果を示す．計算条件を Table 5-1 に、解析領域を Fig.5-1 に示す．縮尺 1/400 の風洞実験を再現して計算した。移流項に TVD（Limitedlinear）スキームを，地表面境界に滑面用一般化対数則を適用している以外は学会指定のベンチマーク条件に準じている．計算格子生成には OpenFOAM 付属の自動格子生成ユーティリティを使用し，学会提供の CAD データを STL 形式に変換したデータを使用して約 34 万要素の格子を自動生成させた．

４．２	  計算結果	 	  Case1 の計算結果を Fig.4-4 に示す．実験値及び機関Ａの計算結果は参考文献 2 から引用した．いずれの値も同じ高さの流入風速で基準化してある．OpenFOAM は実験値に対し、全体としては概ね近い値を示しているが、対象建築物の斜め前方(測定点 9 周辺および測定点 17 周辺)、側面および後流域において乖離が大きくなる傾向が見られる．Case2~Case6 においても側面および後流域において乖離が大きくなる傾向が見られた．OpenFOAM と機関 Aを比較すると、OpenFOAM の方が全体的に大きくなる傾向が見られる．
また、斜め前方は機関Ａの方が実験値によく一致しており、後流域は OpenFOAM の方が実験値に近い値を示している．
計算値と実験値の相関係数及び標準誤差を Table 4-2 に示す．相関係数は概ね 0.8~0.9 の範囲であり、両者の対応は良好である．
公開されている他の CFD 結果(1)と比較すると、相関係数に関しては同等であり、標準誤差に関しては他の CFD 結果(1)よりも小さくなるという結果になった．	 

Table 5-1	  Calculation condition
流入境界条件ベンチマーク指定値(1)を使用
上面境界条件slip 壁側面境界条件slip 壁 (風向 0 時のみ)
流出境界条件
勾配ゼロ
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル
移流項スキーム
TVD（Limited linear）スキーム
計算アルゴリズム SIMPLE，定常解法
計算格子分割法
適合直交格子による 5 段階細分化

Case1

High-risebuilding

Table 4-2	  Correlation factor and standard error
Case	 
対象建築物形状
風向相関係数
標準誤差1	 D D D
00.810.192	 D D D
22.50.880.163	 
D D D
450.850.194	 2D D D
00.880.125	 2D D D
22.50.850.166	 2D D D
450.920.13	 ４．３	  計算結果の補正	 	  風洞実験において、サーミスタ風速計で測定できるのはスカラー風速計の瞬時値の平均値である．一方、k-εモデル等の RANSモデルから求めているのは、各風速成分の平均値から求めたスカラー風速である．そこで Case4 の計算結果を文献(1)の方法により補正し、実験値と比較した． Fig.4-5 に結果を示す．補正後の値は、全体的に風速比が大きくなり、実験値に近づくことがわかる．
補正後の相関係数は 0.89、標準誤差は 0.11 になった．

0.9m

Fig.4-4 Comparison of normalized velocity0o22.5o45oWind

1.8m

1.8m

City blocks

Fig.5-1	  Calculation area５．２	  計算結果	 
Fig.5-2 に示す高さ 5mm(実物換算 2m)での各計測点における風速比に関する CFD 計算結果と実験値の比較を Fig.5-3~5-5 に示す。
ここで M1 や T1 は他の CFD コードによる計算結果であり、M1は構造系直交格子、T1 は非構造系テトラ格子での計算結果である。
なお、風速比の基準風速は高層建物高さにおける流入風速 6.61m/sである。また、実験値はサーミスタ風速計によるスカラー計測値であるため、4 章と同様に CFD 結果は文献(1)による補正を行なっている。どのコードにおいても低風速領域ではバラツキが大きくなる傾向があるが、高速領域は精度良く予測されている。スカラ
ー風速に関する実験値と計算値の標準誤差(RMSE)を Table 5-2 に示す．OpenFOAM の誤差は他のコードよりも小さい結果となったが、この原因は離散化スキームや計算格子の差異によると思われる．

Case4

Fig.4-5	  Comparison of normalized velocity
4Copyright © 2010 by JSFM

1

35 1545 55 65

6 1646 56 66

第 24 回数値流体力学シンポジウム	 
D5-6
会提供の CAD データを STL 形式に変換したデータを使用して約220 万要素の格子を自動生成させた．
Table 6-1	  Calculation condition
流入境界条件ベンチマーク指定値(1)を使用
上面・側面境界条件 slip 壁 (側面は風向 N,E,S,W 時のみ)
流出境界条件
勾配ゼロ
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル
移流項スキーム
TVD（van Leer limiter）スキーム
計算アルゴリズムSIMPLE，定常解法
計算格子分割法
適合直交格子による 5 段階細分化

7 17 25 31 35 39 47 57 67 75 77
8 18 26 32 36 40 48 58 68
9 19 27

41 49 59 69

10 20 28

42 50 60 70

11 21 29 33 37 43 51 61 71
24 12 22 30 34 38 44 52 62 72 76 78
13 2353 63 73

14 24

54 64 74

1.1
1 OpenFOAM
CFD(M1)
0.9
CFD(T1)
0.8
Exp
0.7
0.6
0.5
0.4
0.3
0.2
0.1
005 10 15

Target buildings
B and C

2025

30

35

40

45

50

55

60

65

70

75

1m

Normalized velocity

Fig.5-2	  Measurement Points

80	 

Measurement point numbero

1.1
1 OpenFOAM
CFD(M1)
0.9
CFD(T1)
0.8
Exp
0.7
0.6
0.5
0.4
0.3
0.2
0.1005 10 15

Target building A
Fig.6-1	  Calculation area

202530354045505560657075６．２	  計算結果	 
Fig.6-2 に示す各計測点における風速比の計算値と実測値をFig.6-3,6-4 に示す．なお、計測点 No.47 は、今回作成した格子では建物内に入ってしまったため表示していない。風速比の基準風速は境界層高さ(1m、実スケール 250m)での流入風速 7.8m/s である．
どちらの風向でも CFD の結果は風洞実験の結果と全体的に対応が良好である．いくつかの点(例えば、風向 NNE での No.50, 65,66、風向 W での No.38, 7, 8, 17, 18)では実測値と乖離した傾向を示すが、いずれも対象建物から比較的離れた路地内である。本報では掲載していないが，公開されている他の CFD 結果(1)も概ね同様の傾向を示し、各 CFD 結果の差は小さい．	 
７．新宿副都心高層ビル群を対象とした例（Case	 F	 ）	 ７．１	  計算概要	 	  本章では新宿副都心高層ビル群を対象としたベンチマークケース（Case F）に基づいて作成・公開した計算ケースの結果を示す．
解析領域を Fig.7-1 に，計算条件を Table 7-1 に示す．
移流項に TVD（van Leer limiter）スキームを，地表面境界に滑面用一般化対数則を適用している以外は学会指定のベンチマーク条件に準じている．
計算格子生成には OpenFOAM 付属の自動格子生成ユーティリティを使用し，学会提供の STL データを使用して約 530 万要素の格子を自動生成させた．自動生成に必要な時間は Nehalem Xeon プロセッサ（2.93GHz）１Core を用いて 10 分程度である．全方位を16 風向に分けて計算を実施し，1 方位当たり 6 時間（8 並列）の計算時間を必要とした．	 
７．２	  計算結果	 
各計測点における風向別の風速比を Fig.7-2 に示す．各計測点はFig.7-1 に示す通りである．風洞実験値及び実測値(1)も合わせて示80	 

Measurement point number

oNormalized velocity

Fig.5-4 Comparison of normalized velocity (Wind direction 22.5 )	 
1.1
1 OpenFOAM
CFD(M1)
0.9
CFD(T1)
0.8
Exp
0.7
0.6
0.5
0.4
0.3
0.2
0.1
00
5 10 152m2m
N

Normalized velocity

Fig.5-3 Comparison of normalized velocity (Wind direction 0 )	 

2025303540

45505560

65

7075

80

Measurement point number	 Fig.5-5 Comparison of normalized velocity (Wind direction 45o)	 	 	 Table 5-2	  Root mean square error of normalized velocity
コード	 
風向 0
風向22.5
風向 45
OpenFOAM	 
0.121
0.118
0.105
M1	 0.150
0.134
0.147
T1	 0.135
0.124
0.129	 ６．新潟市内低層建物密集地を対象とした例（Case	 Ｅ	 ）	 ６．１	  計算概要	 	  本章では新潟市内低層建物密集地を対象としたベンチマークケース（Case E）に基づいて作成・公開した計算ケースの結果を示す．計算条件を Table 6-1 に，解析領域を Fig.6-1 に示す．縮尺 1/250の風洞実験を再現して計算した。移流項に TVD（van Leer limiter）
スキームを，地表面境界に滑面用一般化対数則を適用している以外は学会指定のベンチマーク条件に準じている．計算格子生成には OpenFOAM 付属の自動格子生成ユーティリティを使用し，学5Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム	 
D5-6
す．何れの結果も，風向 NE〜N〜NW では参照点 1，それ以外の風向は参照点 2 の風速で基準化してある．No.6・No.7・No.13 に関して，CFD 計算は風洞実験及び実測との対応が良好である結果となった．しかしながら No.15 の風向 NE〜SE では実験値・実測値と乖離した傾向を示す．本方向では掲載していないが，公開されている他の CFD 結果(1)も同様の傾向を示す．これは風向 NE を境に基準化する際の速度参照点が変わる事が影響している可能性がある．因みに本報を含めた各 CFD 結果の差は小さい．
Fig.7-3 に風向 S 時の計測点別の風速比を示す．一部の計測点での乖離はあるものの．
全体的には CFD 結果と実験及び実測との対応は概ね良好と言え．他の CFD 結果(1)との差異も小さい．	 
Table 7-1	  Calculation condition流入境界条件
ベンチマーク指定値(1)を使用
上面境界条件slip 壁側面境界条件
風向 N・S・W・E：slip 壁
その他風向：勾配ゼロ
流出境界条件
勾配ゼロ
固体壁面境界滑面用の一般化対数則
乱流モデル
標準 k-εモデル
計算格子分割法 適合直交格子による 3 段階細分化

Fig.7-1	  Calculation area	 
７．まとめ	 	  本報ではオープンソース CFD ツールキット OpenFOAM を用いて上記 6 ケースの計算を行い、実測結果や既往の CFD 解析結果との比較を行なった結果、OpenFOAM は概ね他のコードと同様の解析精度を有することが分かった。OpenFOAM の計算入力ファイルはオープン CAE 学会の WEB サイト(3)で公開しているので、今後OpenFOAM で風環境の計算をされる方は是非参考にされたい。

Fig.7-2 Comparison of normalized velocity at each wind direction

謝辞	 流体数値計算による風環境評価ガイドライン作成 WG(主査：東北大学 持田灯、幹事：新潟工科大学 富永禎秀)の委員の方々から CFD コードでの計算結果を提供頂いた。ここに深く感謝の意を表する。	 	 
参考文献	 
(1) 市街地風環境予測のための数値流体解析ガイドブック ‐ガイドラインと検証用データベース‐，日本建築学会編, 2007(2) http://www.aij.or.jp/jpn/publish/cfdguide/index.htm(3) http://www.opencae.jp/

Fig.7-3	  Comparison of normalized velocity at each measurement point(wind direction: S)
6Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム	 
D5-6

Fig.6-2	  Measurement Points

Fig.6-3 Comparison of normalized velocity (Wind direction NNE)	 	 	 Fig.6-4 Comparison of normalized velocity (Wind direction W)	 	 
7Copyright © 2010 by JSFM