Comments
Description
Transcript
通風気流場を対象とした CFD 解析の精度検証
通風気流場を対象とした CFD 解析の精度検証 Accuracy Verification of CFD for Cross-Ventilated Flow Field 小林知広 * 1 、甲谷寿史 * 2 Tomohiro KOBAYASHI, Hisashi KOTANI *1 立命館大学理工学部 建築都市デザイン学科 Dept. of Architecture and Urban Design, Ritsumeikan University, *2 大阪大学大学院工学研究科 地球総合工学専攻 Dept. of Architectural Eng., Osaka University キーワード:通風(Cross-Ventilation)、風洞実験(Wind Tunnel Test)、標準 k-e モデル (Standard k-e Model)、応力方程式モデル(Reynolds Stress Model)、LES(Large Eddy Simulation) 1. はじめに 自然通風時の自然換気量予測時には、開口の種類で定まる流量係数に基づいた有効開口面積の直列 結合値と、開口を有しない建物模型の実験等から得られる風圧係数を入力値として換気の式を用いる。 この手法は、大開口の通風時には換気量が過小評価されることや、開口面に気流が斜めから流入する 場合には換気量が過大に評価されるという問題点が指摘されている 1, 2) 。これらの問題は、流速が開口 の前後で十分に小さくなることを前提としてベルヌーイの式に基づく換気の式を静圧差と固定値で与 えられる流量係数で表記するために生じ、この手法は実際の通風現象を表していないと言える。開口 を通過する気流の流量を予測するためには、流路である流管内で実際に生じるエネルギー損失に基づ いた予測手法 3) を用いることが適切と考えられるが、流管内部の気流性状の詳細はこれまでに明らか にされていない。流管の性状を定量的に明らかにすることは、風洞実験では困難であるため数値流体 力学(以降、CFD)が有効と考えられるが、その際実験との比較により解析精度を確認する必要がある。 しかし風洞実験により静圧や風速を測定する場合、測定点の気流特性に基づいて適切な測定手法を選 択する必要がある。そこで、これまでに筆者らは比較的簡易な形状を有する室模型を用いて数種の風 洞実験と CFD 解析を行い、風速・静圧・乱流統計量に関する解析精度検証を行ってきた 4, 5) 。本稿で はその一部を抜粋し、風洞実験による測定の手法と CFD 解析の精度検証を紹介する。 2. 風洞実験概要 2.1 流管中心線における流体量測定 実 験 は 測 定 洞 長 さ 9.5m、 高 さ 1.8m、 幅 1.8m の 回 流 式 風 洞 で 行 っ た。 実 験 に 用 い た 模 型 は 図 1 に 示 す 9 室 の 単 室 か ら 構 成 さ れ た 建 物 と し、 中 央 の 1 室 の み 対 面 開 口 が 開 放 さ れ て い る 状 況 を 想 定 し た。 開 口 を 通 り 抜 け る 流 管 の 基 本 的 な 性 状 を 明 ら か に す る た め に、 模 型 は 風 洞 中 央 に 設 置 し て 10m/ s の一様流にさらし、流管内部の風速と静圧の測定を行う。この際、風洞内風速の設定は模型の風上 600mm、風洞壁面より 300mm の位置に設置したピトー管により設定し、風洞内基準静圧もピトー管 位置の静圧を採用した。図 2 に示す通り、測定は比較的容易な開口部の中心線上で行い、トラバーサ Pitot Tube X L 600 200 1,400 200 Static Pressure Tube or Hot Wire Anemometer Model 300 Support (t=10) Model Traverser Pitot Tube Support 200 L Approach Flow (10 m/s) 1,800 Y Leeward 200 4 1,000 d Win 38 384 Z Windward 1,400 .6 181 (1)X-Y 断面 図 1 建物模型 (単位 [mm]) (2)Y-Z 断面 図 2 風洞断面図 に 設 置 し た ア ー ム に 取 り 付 け た I 型 熱 線 流 速 計(0251R-T5 型、φ 0.5μm タ ン グ ス テ ン、KANOMAX 社)または静圧管により X=-150mm から 1,000mm まで 10mm 間隔で風速と静圧の測定を行った。但し、 開口の前後 50mm の範囲は 5mm 間隔で測定した(図 3)。実験条件は開口寸法 L とし、15、30、45、 60、90mm の 5 条件を設定した。図 4 にプローブヘッドの形状を示す。静圧管は自作のものを使用し、 先端から 9.0mm の位置に測定孔を 4 カ所設けている。測定点が流管の中心線であるため、平均風速は プローブに平行となる。図 5 に静圧管を 2 毎に回転させたときの指向特性を示す。ここで、縦軸は各 風向における静圧読み値とプローブが気流に平行な場合(0 )の静圧読み値との差をアプローチフロー の速度圧で無次元化した誤差率であり、静圧 10 の範囲であれば誤差率は概ね 2% 以内に収まる。静 圧測定は微差圧計(MP-45、Validyne 社)を用いて 100Hz、30 秒間を行い、風速は後に乱流エネルギー を算出するために 1.0kHz で 30 秒間測定を行った。 Model 100 Leeward 10.0 9.0 Static Pressure Tube Center Line 4.0 2.0 100 Dimensionless Pressure [-] 5.0 Measured Point 2.11 Windward Measuring Point 181.6 150 1,100 8.0 Hot Wire Anemometer [mm] 図 3 流管中心線上の測定点 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -90 図 4 プローブヘッド -60 -30 0 30 Roration Agle [degree] 60 90 図 5 静圧管の指向特性 2.2 PIV を用いた建物外側の流れ場の測定 室外側の定性的な流れ場を測定するために、前節と同様の模型を用いて Particle Image Velocimetry ( 以 降、PIV) に よ る 測 定 を 行 っ た。PIV と は レ ー ザ ー 光 を 用 い て 気 流 の 可 視 化 を 行 い、 微 小 な 時 間 間 隔 で 2 枚 の 画 像 を 撮 影 し、 ト レ ー サ の 移 動 量 か ら 流 速 を 測 定 す る 手 法 で あ る。 可 視 化 は 煙 発 生 機 (KANOMAX 社、8304)で発生させた煙を風洞上流部に注入することで行った。光源にはダブルパル ス Nd:YAG レーザー(NEWWAVE RESEARCH、DPIV-N50)を用い、風洞側面から中央の水平断面を 投影し、上部に設置した CCD カメラ(LaVision 社、ImagerIntense)により撮影した(図 6)。撮影は 4.0 Hz で 25 秒間行い、1 回の撮影で 200 ms の間隔で 2 枚の画像を撮影し、1 回の測定で 100 組 200 枚の 画像を撮影した。撮影領域は 625 mm 471 mm としたが、模型 の奥行きに対して見付面積が比較的大 きいことから後流の複雑な気流場も広い領域に渡って形成されると考えられるため、図 7 に示すよう に模型周辺と後流域の 2 つの領域において測定を行った。画像 の制御もこれによった。また、風速算定には FFT 相互相関法 6) を用いた。ここでは、統計的信頼性と空間分解能の両者を向上 させるために再帰的相関法 を用いており、瞬時の速度場を計 7) 算するための相互相関係数分布の算出回数は 3 回とした。表 1 に PIV 測定概要をまとめて記す。 CCD Camera 200 Y Y Test Model X Smoke Injector Laser Sheet Double Pulse Nd : YAG Laser Z Program Windward Area 1 Area 2 Leeward Z = 401 mm Wind Z = -69 mm 200 1,800 Approach Flow (10 m/s) CCD Camera Leeward 1,400 Windward 625 mm × 471 mm 1,376 pixel × 1,040 pixel Davis 7.2 FFT Cross-Correlation Method Algorithm (Recursive Correlation Method) Pass 1 : 64 pixel × 64 pixel Interrogation Pass 2 : 32 pixel × 32 pixel Window Size Pass 3 : 32 pixel × 32 pixel Overlap 50 % Total Numver of Vectors 5,590 (86 65) Time Interval of Pulse 200 μs Sampling Frequency 4.0 Hz Sampling Time 25 seconds 50 mJ / Pulse Laser Output 0.3 - 1.0 μm Seeding Size Camera Frame Size 1,800 処理には Davis 7.2(LaVision 社)を使用し、カメラとレーザー 表 1 PIV 測定概要 X = -154 mm [mm] (1)X-Y 断面 (2)Y-Z 断面 図 6 PIV 測定時の風洞断面図 X = 768 mm Test Model 図 7 PIV 測定範囲 2.3 建物周辺の風速及び静圧分布の測定 2.3.1 実験目的 前節では建物周辺の定性的な気流場に着目した測定を目的としたが、本節では建物周辺の圧力分布 とその解析精度に着目した測定を行う。また、風速測定は測定箇所の気流性状に応じて適切な測定手 法を選択する必要があると言え、PIV 測定に関しても測定範囲の中で測定精度が高い箇所と低い箇所 が混在する可能性があるため、測定精度が低下する箇所の確認とそのデータを補う目的で、5 孔ピトー 管を用いて静圧・風速・風向の同時測定を行う。 2.3.2 プローブの較正実験 図 8 に本研究で用いた 5 孔ピトー管のプローブを示す。プローブ先端の球形センサは外径 4.0 mm であり、センサは 5 つの静圧測定孔(以下、静圧孔)を有する。図中の静圧孔 1 から 4 は球の中心か ら測定孔 0 を通る軸に対して 45 に配置されている。5 孔ピトー管の作動原理としては、微差圧計を 用いて静圧孔 1 から 4 における圧力の読み値が同じ値を示すようプローブを 3 次元で回転させた場合 に静圧孔 0 を通る軸が速度ベクトルの方向と一致するというものである 8) 。このときの静圧孔 0 にお ける圧力の読み値に基づいて流速を算出することができる。この手法は測定精度が高く、実験データ の処理も非常に簡易ではあるが、センサをプローブごと回転させる必要があるため、縮小模型や測定 機器に接触する等の問題があり採用は難しい。このような背景の中で Wright 9) はプローブを回転させ ず、かつ比較的容易な較正手法を提案した。本研究ではこの較正手法を用いて測定を行う。 Y Y Z 0.4 90 [deg] 1 0 4 2 4.0 ・頂角φ:風速ベクト ルと X 軸のなす角 0.5 2.5 X 3 2.0 Dimensions in mm (1)Y-Z 立面図 ・夾角δ:風速の YZ 成分と XY 面のなす角 2.0 (2)X-Y 断面図 図 8 センサヘッドの形状 図 9 角度の定義 9) この手法では、まず 5 つの圧力の読み値から風向を算出する。風向は図 9 に示す頂角 f と夾角 d に より定まる。頂角は風向ベクトルと X 軸のなす角であり、夾角は風速ベクトルの YZ 成分と XY 面の なす角である。各静圧孔における圧力の読み値から風向、風速、静圧を求めるにあたり、Wright は以 下の式で表される 3 つのパラメータ(角度係数 K f、速度係数 K V 、圧力係数 K p )を提案した。 1/2 4 ⎡ ⎧ 4 ⎫ ⎤ Kφ = ⎢1 − ∑ ( p0 − pn ) / 2 ⎨∑ ( p0 − pn )2 ⎬ ⎥ ⎩ n =1 ⎭ ⎥⎦ ⎢⎣ n = 1 −1 / 2 ⎡ ⎧4 ⎫ ⎤ KV = ⎢ ρV 2 ⎨∑ ( p0 − pn )2 ⎬ ⎥ ⎩ n =1 ⎭ ⎥⎦ ⎢⎣ 1 / 2 (1) 1 / 2 K p = 2( p0 − pS ) / ρV 2 (2) (3) こ こ で、 p 0 は 静 圧 孔 0 の、 p n は 静 圧 孔 1 か ら 4 に お け る 圧 力 の 読 み 値 で あ る。 V は ス カ ラ ー 流 速 で、 p S はプローブ設置位置の静圧である。5 つの圧力読み値から頂角、夾角、流速、静圧を得るためには 関係式がもう一つ必要となるが、夾角 d が以下の式で与えられる。 tanδ = −( p2 − p4 ) / ( p1 − p3 ) (4) 式(1)から式(3)の 3 つのパラメータは頂角 f により概ね定まり夾角 d によっては大きく変化しな いことが示されており 9) 、あらかじめ較正実験により角度係数、速度係数、圧力係数の較正曲線を f の関数として与えることで風向、風速、静圧の同時測定が可能となる。 こ の プ ロ ー ブ を 風 洞 中 央 で 10 m/s の 一 様 流 に さ ら し、 ア ー ム に 取 り 付 け た 分 度 器 に よ り 水 平 方 向 にのみ管を回転させて d を 90 に固 定 し、f を 90 の範囲で 2 ずつ変化さ せて各静圧孔にお いて圧 力を測定した。圧力の測定には微差圧計(MP45 及び DP103 ( バリダイン社))を用いて 100 Hz で 60 秒間行った。また、圧力係数算出のため、風洞中央に設置した JIS 型ピトー管(LK-1S、岡野製作所) で測定点の静圧も別途測定した。これらの結果から得られた角度係数、速度係数、圧力係数を、それ ぞれ 4 次式、6 次式、4 次式で近似を行った較正曲線と併せて図 10 に示す。この較正曲線を用いて次 節に示す風向、風速、静圧の測定実験を行う。 Pressure Factor [-] Angle Factor [-] 1.5 Kφ = 2.80 × 10 −8 φ 4 − 7.45 × 10 −6 φ 3 + 4.85 × 10 −4 φ 2 + 1.00 × 10 −2 φ + 1.18 × 10 −2 1.2 1 0.8 0.6 0.4 1 0 -0.5 15 30 45 60 Rotation Angle [degree] 75 -1.5 90 1.4 1.35 1.3 1.25 1.2 1.15 1.1 -1 0 KV = 2.44 × 10 −11 φ 6 − 5.56 × 10 −9 φ 5 + 5.20 × 10 −7 φ 4 − 2.60 × 10 −5 φ 3 + 6.91 × 10 −4 φ 2 − 6.97 × 10 −3 φ + 1.12 1.45 0.5 0.2 0 1.5 K p = 6.21 × 10 −8 φ 4 − 2.25 × 10 −6 φ 3 − 5.99 −4 φ 2 + 3.80 × 10 −3 φ + 0.938 Velocity Factor [-] 1.4 1.05 0 15 (1)角度係数 30 45 60 Rotation Angle [degree] 75 1 90 0 15 30 45 60 75 90 Rotation Angle [degree] (2)速度係数 (3)圧力係数 図 10 各係数の較正曲線 2.3.3 5 孔ピトー管を用いた静圧・風向・風速の同時測定 前述の実験と同様に建物モデル(開口寸法 L は 45mm)を風洞の中央に設置し、10 m/s の一様流下 で実験を行った。測定点は図 11 に示す中央の水平断面上に設定した。模型の側面(剥離領域)では センサがアプローチフローに正対するよう設置したが、逆流が生じている可能性もあるため、黒塗り プロット位置では風向 120 に正対する状態でも測定を行った。また、風上側の測定点では壁面に沿う 気流でも測定ができるよう風向 15 に正対するよう回転させて測定を行った。風向算定は角度係数の 較正曲線に基づいて頂角 f が 75 の範囲で行った。図 12 に実験時の風洞断面を示す。風洞内基準静 圧はピトー管とし、測定は較正実験と同様の微差圧計を用いて 100 Hz で 60 秒間行った。 +90 deg Approach Flow (10 m/s) 30 30 30 30 30 30 30 30 図 11 5 孔ピトー管による測定点 1,800 30 30 30 30 30 30 0 deg 1,400 3D Traverser 200 3D Traverser Protractor Protractor Test Model Five-hole Pitot tube Pitot tube Test Model Five-hole Pitot tube Pitot tube Dimensions in [mm] (1)X-Y 断面 1,400 -90 deg 200 200 Ceiling 180 deg 30 30 30 30 30 30 30 30 30 30 Wind direction Leeward 200 Windward (2)Y-Z 断面 図 12 5 孔ピトー管を用いた測定時の風洞断面図 3. 実験を再現した CFD 解析の概要 3.1 解析領域及び解析手法 風洞実験を再現した CFD 解析を行った。ここでは解析結果に大きく影響すると考えられる乱流モデ ルに関して、標準 k-e 2 方程式モデル(以降、SKE)、応力方程式モデル(Reynolds Stress Model、以降、 RSM)、Large Eddy Simulation( 以 降、LES) の 3 種 を 用 い て 解 析 を 行 っ た。 な お、RSM は Launder ら 10) の提案に基づく LRR モデルを使用したが、乱流拡散項は等方性の渦粘性近似により簡易化した モ デ ル を 用 い た。LES で は Smagorinsky-Lilly モ デ ル を 用 い、Smagorinsky 定 数 は 0.1 と し て 与 え た。 解析概要を表 2 及び表 3 に示す。解析対象空間は風洞の断面を再現し、模型の風上側に 1,000mm、後 流域は 4,000mm とした。なお、RANS モデルである SKE 及び RSM では計算負荷を軽減するため対称 面を 2 面設けて空間の 1/4 のみ解析を行った。RANS モデルの計算に用いたメッシュレイアウトを図 13 に示す。また LES を用いた解析においては、SKE の解析結果を初期条件として行った。解析空間 の大きさは RANS と同様としたが、非定常計算であるため、解析領域は設定した解析空間全体とした。 計算時間間隔は 0.0001 秒 (10 kHz) として合計 2.3 秒間の計算を行った。図 14 に全計算期間内の流入 開口、流出開口、後流域の 3 点における主流成分の風速と静圧の変動を示す。この結果から、8,000 タ イムステップ (0.8 秒 ) までの結果を SKE からの移行期間とみなして助走計算とし、本計算は 15,000 タイムステップ(1.5 秒)とした。後述する LES の時間平均値は本計算 1.5 秒間の平均とする。 LES (Large Eddy Simulation) RANS (Reynolds Averaged Navier-Stokes Equation) CFD Code Fluent 6.3 Standard k-ε Model (SKE) Reynolds Stress Model (RSM) ࠉ (SIMPLE) Implicit Method Turbulence Model Algolithm Discretization Scheme for Advection Term CFD Code Turbulence Model Algolithm Discretization Scheme for QUICK Advection Term Time Step Volocity : 10m/s Turbulent Intensity : I=1% Total Calculation Term Turbulent Length Scale : Λ=126mm Inlet Gauge Pressure : 0 [Pa] Boundary Condition Outlet Walls : Generalized Log Low Walls Symmetry : free slip Smagorinsky Coefficient L=15mm : 1,612,128 L=30mm : 1,612,648 Total Number of Cells L=45mm : 1,685,808 L=60mm : 1,685,808 L=90mm : 1,619,688 Inlet Boundary Condition Outlet Walls Total Number of Cells Static Pressure [Pa] 表 3 LES の解析概要 Fluent 6.3 Large Eddy Simulation (Smagorinsky-Lilly Model) ࠉ (SIMPLE) Implicit Method Central Differencing 0.0001s (10kHz) 23,000 time step (2.3sec.) Velocity : 10m/s (Constant) Gauge Pressure : 0 [Pa] Warner & Wengle’s linear-power law 0.1 L=15mm : 1,893,060 L=30mm : 2,050,240 L=45mm : 2,024,400 L=60mm : 2,026,160 L=90mm : 2,028,240 Velocity [m/s] 表 2 RANS の解析概要 80 60 40 20 0 -20 -40 -60 -80 16 14 12 10 8 6 4 2 0 -2 -4 -6 Static Pressure X=0 mm 0 0.2 0.4 0.6 0.8 1 X=181.6 mm X=1,000 mm 1.2 Time [s] 1.6 1.4 1.8 2 2.2 2.4 2 2.2 2.4 X-Velocity X=0 mm X=181.6 mm X=1,000 mm 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] 1.4 1.6 1.8 図 14 LES 解析結果の時変動 181.6 4,000 Symmetry 1,800 1,800 Model Inlet Symmetry Y Z 1,800 900 900 1,000 Model Outlet Wall [mm] X Wall Y X (1)X-Y 断面 [mm] Z (2)Y-Z 断面 図 13 CFD における解析領域とメッシュレイアウト(RANS) 3.2 境界条件 RANS を 用 い た 解 析 に お け る 流 入 境 界 条 件 は、 実 験 に 基 づ き 平 均 風 速 U =10.0m/s、 乱 流 強 度 I =1.0%、乱れの長さスケール L=126mm(風洞の水力直径の 0.07 倍)として与え、SKE において k と e の流入境界条件は以下の式に基づいた値を適用した。 3 k = (U I) 2 2 (5) = CD 3 / 4 k3/2 (6) RSM におけるレイノルズ応力の流入境界条件は以下の式で与えた。 2 ui ui = k 3 (i = 1, 2 3) (7) ui u j = 0 (i, j = 1, 2, 3) (8) 一般に、LES を用いた解析では変動する流入境界条件の作成が重要となるが、本研究で解析対象とし た風洞実験ではラフネスや格子などにより乱れを作成しておらず、アプローチフローは乱れの少ない 一様流であるため、主流方向の流速成分を 10.0 m/s の固定値として与えた。流出境界条件は全解析条 件において圧力規定とし、後流域を長く設けていることから圧力分布は大きくないと判断してゲージ 圧 0Pa を一様に与えた。壁面境界条件は対象面を free-slip 条件とし、RANS 解析において風洞及び模 型の壁面は一般化対数則を適用した。一方、LES を用いた解析では Werner & Wengle 11) による linearpower law 型の 2 層モデルを適用した。 4. 風洞実験と CFD の比較による解析精度検証 4.1 流管中心における速度と圧力の比較 図 15 に風洞実験と CFD から得られた L=15, 45, 90mm 条件における中心線上の静圧、風速、乱流 エネルギーを示す。ここで、静圧と風速はそれぞれ最も風上側の全圧とアプローチフローの風速で基 準化して表している。また、CFD の静圧は、風洞実験でピトー管を設置した位置の静圧を基準として 表記し、基準静圧の位置を一致させている。静圧の分布性状は風上側の衝突域で気流がよどみ、ベルヌー イの式に基づいて速度圧が静圧に変換されて上昇し、開口を通過する際に低下した後に風下側で緩や かに上昇する。風上側ではどの乱流モデルも実験結果と良く一致しているが、室内部においては SKE が静圧を過大評価していることがわかる。これは SKE は流入後に室内で流速があまり上昇しないため 静圧が速度圧に変換されていないことを表し、その原因の一つとして、速度勾配に応じて生成される 乱流エネルギーを過大評価する点があげられる。風下側の静圧は SKE は実験と比較して高い圧力を示 しているが、RSM と LES では比較的実験に近い結果が得られた。 風速に関しては、前述の通り SKE が室内で過小評価されるが、RSM と LES では実験と良く一致し たと言える。風下側に着目すると、実験と CFD で大きな差異が見られるが、これは実験では I 型熱線 流速計を用いているため風速 2 成分の絶対値として出力することに対して、CFD では主流成分のみを 正負で評価するためであり、逆流が存在すると考えられる Wake 内部では厳密に比較することは難し いが、RSM は拡散して流速が低下する点が後側にずれている傾向が見られる。これは RSM が剥離域 や後流域のサイズを過大評価することを示しており、これまでに大岡ら 12) が示した傾向と一致する。 SKE と LES の 後 流 域 の 平 均 風 速 は 概 ね 一 致 し て い る が、 こ れ ら の 解 析 結 果 の 精 度 を 検 証 す る た め、 LES の瞬時値に基づいて熱線流速計で出力される瞬時値を以下の式に基づいて再現する。 vLES hotwire = v 2x + v 2y (9) これは I 型熱線流速計の出力がワイヤーに直交する瞬時風速 2 成分の合成ベクトルと仮定するもので あるが、このように LES で熱線出力を再現した瞬時値に本計算期間内で時間平均を施した結果を図中 LES RSM Measurement 0 1 2 3 4 5 X / Model Length [-] 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 LES RSM Measurement 0 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 1 2 3 4 5 X / Model Length [-] Measurement 0.6 0.4 0.2 SKE 0 LES -0.2 RSM -1 0 1 2 3 4 X / Model Length [-] 1.2 1 0.8 Measurement 0.6 0.4 LES-Hotwire 0.2 0 -1 Measurement 0 1 2 3 4 X / Model Length [-] 1 2 3 4 X / Model Length [-] 5 6 5 6 1.4 1.2 Measurement 1 LES-Hotwire 0.8 0.6 SKE 0.4 0.2 0 LES -0.2 RSM -0.4 0 LES RSM SKE -0.2 -1 0 1 2 3 4 X / Model Length [-] 5 Turbulent Kinetic Energy 18 16 14 6 LES SKE 12 10 8 6 RSM 4 2 LES-Hotwire Measurement 0 6 1.4 6 RSM -1 5 -0.4 SKE LES LES-Hotwire 0.8 6 SKE -1 1 -0.4 Dimensionless Velocity [-] -1 Turbulent Kinetic Energy [m2/s2] SKE 1.2 -1 Turbulent Kinetic Energy [m2/s2] Dimensionless Velocity [-] X-Velocity 1.4 2 2 Turbulent Kinetic Energy [m /s ] Static Pressure 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 Dimensionless Velocity [-] L=15 mm L=45 mm Dimensionless Static Pressure [-] Dimensionless Static Pressure [-] L=90 mm Dimensionless Static Pressure [-] に LES-Hotwire として記載しているが、この結果を見ると LES は実験と良く一致したと言え、主流成 0 1 2 3 4 X / Model Length [-] 5 6 18 16 14 12 LES SKE 10 RSM 8 6 4 2 LES-Hotwire Measurement 0 -1 0 1 2 3 4 X / Model Length [-] 5 6 18 16 14 LES 12 10 8 SKE RSM 6 4 2 Measurement LES-Hotwire 0 -1 0 1 2 3 4 X / Model Length [-] 図 15 風洞実験と CFD から得られた中心線上の静圧・風速、乱流エネルギー 5 6 分の平均風速が似た傾向を示す SKE においても後流域の風速に関しては解析精度は比較的良いと考え られる。また、乱流エネルギーを比較すると、前述の通り SKE では流入前後で過剰に生成されており、 これは一般に言われる通り、通常の渦粘性モデルでは速度勾配が大きい場合に乱流エネルギーの生産 項が大きくなるためである デルも複数提案されており 13) 。また、このような乱流エネルギーの過剰生成を抑制する改良型 k -e モ 14, 15) 、今後それらのモデルを適用した際の解析精度の検証も課題と言える。 RSM では流れ場に対応して乱流エネルギーのピークも後流側にシフトしているが、その値としては実 験値と概ね同程度の結果が得られた。熱線を再現した LES の結果は開口が大きくなった際に若干過小 に評価しているが、実験結果と比較的良く一致する結果が得られた。 4.2 建物周辺の気流場の解析精度検証 図 16 に PIV と CFD に よ り 得 ら れ た 建 物 モ デ ル 周 辺 の 風 速 ベ ク ト ル を 示 す。 こ こ で、PIV の 風 速 ベクトルは縦横ともに 1/4 に省略して示しているため、表示ベクトル数は全体の 1/16 である。また、 CFD 解析結果のベクトル表示点はメッシュレイアウトに関わらず、PIV と合わせて表示する。なお、 PIV は 4.0Hz、25 秒間の時間平均値であり、LES は 10 kHz、1.5 秒間の時間平均値である。PIV の測 定結果において、剥離角度に開口条件による大きな差異は見られないことや、後流域で建物斜め後方 に形成される大きな渦による逆流が生じて室からの流出気流と衝突していることがわかる。また、気 流は衝突後に風上側に向かって流れ、剥離流と合流して風下側または後流の渦へと流れている様子が わかる。室内の PIV 測定結果に関しては、Z 方向に対称となっていないが、これは後流側かつ Z が正 の方向からレーザーが照射され、2 枚のアクリル板を透過していることが影響したと考えられ、この ため室内の PIV 測定に関しては精度が低下すると言える。PIV 測定と CFD 解析の結果を比較すると、 RSM では建物斜め後方の渦が再現されておらず、剥離域の大きさを過大に評価している。これは、流 L=15 mm 350 300 300 250 250 250 200 200 200 150 100 100 50 50 50 0 0 0 -50 -50 -50 0 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 0 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 400 350 350 300 300 300 250 250 250 200 200 200 150 Z [mm] 400 350 Z [mm] 400 150 50 50 0 0 0 -50 -50 -50 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 0 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 400 350 350 300 300 300 250 250 250 200 200 200 Z [mm] 400 350 150 100 100 50 50 0 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 0 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 400 400 350 350 350 300 300 300 250 250 250 200 200 200 100 Z [mm] 400 150 150 100 0 0 0 -50 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 100 -50 -100 -200 -150 -100 -50 0 150 50 50 50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 50 0 -50 0 0 100 0 -100 -200 -150 -100 -50 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 150 -50 -100 -200 -150 -100 -50 0 50 400 150 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] 100 100 0 0 150 100 Z [mm] Z [mm] SKE RSM Z [mm] 150 100 -100 -200 -150 -100 -50 LES Z [mm] 350 300 Z [mm] 350 150 L=90 mm 10 m/s 400 -100 -200 -150 -100 -50 Z [mm] L=45 mm 10 m/s 400 Z [mm] PIV Z [mm] 10 m/s 400 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 X [mm] -100 -200 -150 -100 -50 図 16 PIV と CFD から得られた建物周辺の風速ベクトル 管中心線の風速で見られた後流域の過大評価と一致する。RSM を用いると剥離領域を過大評価すると は必ずしも言えないが、大岡ら 12) や村上 13) も同様の結果を示している通り、現在一般に用いられる RSM のモデリング手法では同様の傾向が得られると考えられる。一方で、RSM は 6 種のレイノルズ 応力の輸送方程式を解いているため、原理的には合理的なモデルと言え、レイノルズ応力の輸送方程 式中のモデリング手法やモデル定数を改良することにより、精度の向上を期待することができると考 えられる。また、SKE で得られた気流場は RSM と比較すると傾向は実験に近いと言えるが、剥離域 の大きさや剥離流の角度は実験とやや異なる。LES の解析結果は L=90 mm の条件でのみ X=500 mm より後流での逆流が生じておらず実験結果と一致しないが、その他の開口条件の気流場では実験結果 と定性的には良く一致したと言える。これらの結果から、建物モデル周辺の気流場の性状としては、 RSM の解析精度が最も低く、SKE、LES では比較的良い結果が得られたと言える。本節では建物モデ ル周辺の定性的な気流のパターンに着目した検証を行ったが、次節では測定手法の信頼性を含め、速度・ 圧力分布に関してより定量的な検証を行う。 4.3 建物周辺の圧力分布及び速度分布の解析精度検証 4.3.1 建物モデル風上側の解析精度 図 17 に L=45mm の 建 物 モ デ ル で 5 孔 ピ ト ー 管 と CFD に よ り 得 ら れ た 建 物 モ デ ル 風 上 側(X=-30, -60, -90 mm) に お け る 静 圧、 ス カ ラ ー 風 速、 風 向 を 示 す。 こ こ で、 ス カ ラ ー 風 速 と 風 向 に 関 し て は PIV 結果も併せて示す。風速はアプローチフローの風速で基準化し、静圧は同時測定したピトー管の 速度圧で無次元化した後に X=0 mm、Z=432 mm(壁面からの距離(Z')が 240 mm)での無次元全圧 で基準化を行った。静圧は SKE のみ過小に評価しているが、RSM と LES の結果は実験結果と良く一 致している。スカラー風速に関しては、3 種 の CFD 結果は概ね一致したが、PIV と 5 孔ピトー管の結 果には差異が見られる。このため、測定箇所により各測定手法の精度が低下すると考えられるが、PIV に関しては模型付近はレーザーシートがアクリル板を透過してトレーサを可視化するため影が発生す ることに加え、模型上部が写り込むという問題点があるため、X=-30mm のライン上では測定精度が低 下したと考えられる。一方、X=-60、-90mm のライン上の PIV 結果は誤差を生む要因が無く信頼性が 高いと考えられ、この領域ではどの CFD 解析結果も良く一致したと言える。5 孔ピトー管のスカラー 0.6 0.4 Five-Hole SKE RSM LES Z [mm] 150 180 1 120 Z [mm] 150 180 0.4 Five-Hole SKE RSM LES 0 0 30 60 30 60 90 120 Z [mm] 150 180 120 Z [mm] 150 180 210 15 0 Five-Hole SKE RSM LES PIV -15 -30 0.8 0 30 60 90 90 120 Z [mm] 150 0.6 0.4 0.2 45 30 15 0 Five-Hole SKE RSM LES PIV -15 -30 -60 30 60 90 120 Z [mm] 150 180 210 1 0 30 60 90 Five-Hole SKE RSM LES PIV 0.8 90 120 Z [mm] 150 0.4 0.2 45 30 15 0 Five-Hole SKE RSM LES PIV -15 -30 -60 60 210 60 -45 30 180 75 0.6 0 210 60 -45 0 180 75 0 90 30 210 Five-Hole SKE RSM LES PIV 1.2 0.6 45 -60 1 210 60 -45 0 0 90 0.8 0.2 0.2 1.2 Five-Hole SKE RSM LES 60 0.4 210 0.4 30 0.6 Wind Direction [Degree] 120 0.6 0 75 0 90 0.8 0.2 1 0.8 Wind Direction 90 Five-Hole SKE RSM LES PIV Wind Direction [Degree] Dimensionless Static Pressure [-] X=-60 mm X=-30 mm 60 1 0 Dimensionless Static Pressure [-] 30 Dimensionless Velocity [-] 0 0 Line where results are indicated Line where measurement with five-hole pitot tube was conducted Dimensionless Velocity [-] 0.8 0.2 Velocity Magnitude 1.2 Wind Direction [Degree] Static Pressure 1 Dimensionless Velocity [-] X=-90 mm Dimensionless Static Pressure [-] 風速は X=-90 mm では CFD より小さくなるが、X=-30 mm では CFD と良く一致し、風向に関しては 90 120 Z [mm] 150 180 210 0 30 60 図 17 PIV と CFD から得られた建物風上側の静圧・スカラー風速・風向 90 120 Z [mm] 150 180 210 模型付近で CFD 結果と異なる。このように 5 孔ピトー管の風向と風速の精度が風上側で悪くなったこ との最も大きな原因は圧力勾配と言える。つまり、角度係数・速度係数・圧力係数の較正曲線は一様 な圧力分布の中で得られたことに対して、衝突域ではそもそも 5 つの測定孔周辺の静圧に分布がある ためである。それにも関わらずこの領域で静圧が CFD 結果と良く一致した理由としては、(3) 式で表 される圧力係数算出に周囲 4 点の静圧孔の測定値は用いられず、中央の 1 点のみの静圧で決定される ため、センサヘッド周辺の圧力分布が重要でなくなるためと言える。このため、衝突域では 5 孔ピトー 管は静圧測定の精度は高く、風向・風速に誤差を生じ、PIV では透過や写り込みの影響が見られる壁 面のごく近傍以外の風向・風速は精度良く測定することができると言える。その上で風上側の CFD 解 析精度に関しては、静圧は SKE のみ精度が少し低く、風向とスカラー風速 はどの乱流モデルでも良い 結果が得られたと言える。 4.3.2 建物モデル側方の解析精度 同様に、図 18 に建物モデル側方の結果を示す。ここで、5 孔ピトー管のプローブ角度で 0 と 120 が存在する点では、予め PIV 結果と比較して適切な測定結果を選択した。静圧に関しては、前述の理 由により 5 孔ピトー管の測定結果は信頼できるものと言え、RSM の結果が実験結果と最も良く一致し た。SKE と LES で は値は やや 異な るも のの、傾向は良く再現できている。スカラー風速に関しては、 X=0mm のラインのみ 5 孔ピトー管と PIV の結果で大きな差異が見られる。ここでは PIV の測定精度 が低下したと考えられ、その理由としては、建物の角部周辺であることから Interrogation Window 内 部で平面的な速度分布が大きくなり空間分解能が不足している可能性があげられる。この領域の PIV 測定精度を向上させるためには、局所的に Interrogation Window のサイズを小さくすることが有効と 考えられる。その他の測定ライン上では 2 つの実験結果は概ね一致しているため測定結果の信頼精度 は高いと考えられる。剥離域内部では乱れの影響で 5 孔ピトー管の風向算定精度がわずかに低下した -1.2 -1.4 -1.6 X=120 mm Dimensionless Static Pressure [-] 60 90 120 Z' [mm] 150 180 210 -0.8 -1 Five-Hole SKE RSM LES -1.2 -1.4 -1.6 0 120 Z' [mm] 150 180 210 -0.8 -1 Five-Hole SKE RSM LES -1.2 -1.4 -1.6 Dimensionless Static Pressure [-] 90 0 120 Z' [mm] 150 180 210 -0.8 Five-Hole SKE RSM LES -1.6 30 60 90 120 Z' [mm] 150 180 210 0.2 30 60 90 120 Z' [mm] 150 180 210 0.6 Five-Hole SKE RSM LES PIV 0.4 0.2 90 120 Z' [mm] 0 150 180 210 Five-Hole SKE RSM LES PIV 0.4 0.2 120 Z' [mm] 150 180 210 120 Z' [mm] 150 180 210 240 Five-Hole SKE RSM LES PIV 150 120 90 60 30 0 30 60 90 120 Z' [mm] 150 180 210 240 Five-Hole SKE RSM LES PIV 180 150 120 90 60 30 0 0 30 60 90 210 0.6 90 90 -30 1 60 60 180 240 0.8 30 30 0 1.2 0 0 210 1 60 30 -30 0.8 30 60 240 1.2 0 90 240 Five-Hole SKE RSM LES PIV 0.4 0 240 210 120 210 0 -0.6 -1.4 180 0.6 240 -0.4 -1.2 150 0.8 1.4 -1 120 Z' [mm] 1 0 -0.2 0 90 1.4 -0.6 60 60 0 -0.4 30 30 Five-Hole SKE RSM LES PIV 150 -30 1.2 240 -0.2 0 0 Dimensionless Velocity [-] -0.6 90 0.2 1.4 -0.4 60 Five-Hole SKE RSM LES PIV 0.4 240 -0.2 30 0.6 0 Dimensionless Velocity [-] X=60 mm Dimensionless Static Pressure [-] 30 0 Wind Direction [Degree] Five-Hole SKE RSM LES 1 0.8 Wind Direction [Degree] -1 1.2 180 Wind Direction [Degree] -0.8 Wind Direction 210 240 Wind Direction [Degree] -0.6 0 X=180 mm Dimensionless Velocity [-] -0.4 0 Line where results are indicated Line where measurement with five-hole pitot tube was conducted Velocity Magnitude 1.4 Dimensionless Velocity [-] X=0 mm Dimensionless Static Pressure [-] Static Pressure 0 -0.2 120 Z' [mm] 150 180 210 240 Five-Hole SKE RSM LES PIV 180 150 120 90 60 30 0 -30 0 30 60 図 18 PIV と CFD から得られた建物側方の静圧・スカラー風速・風向 90 120 Z' [mm] 150 180 210 240 ものの、平均化時間を比較的長く設定したため、平均風向の算定結果に及ぼした誤差は大きくない。 以 上 よ り、 建 物 モ デ ル 側 方 に 関 し て は X=0mm の PIV 結 果 と 剥 離 域 内 部 の 5 孔 ピ ト ー 管 に よ る 風 向 算定結果を除いて測定結果の信頼性は高いと考えられる。その上で CFD 解析結果を見ると、RSM と LES では実験と良く一致するが、SKE では等方性の渦粘性近似を用いていることにより乱流拡散が大 きくなり速度勾配を小さく評価したと考えられる。また、風向に関しては 90 をまたいで大きく変化 する位置は剥離域の渦中心と考えられるが、SKE ではその位置が壁面に近い位置に見られる。このた め、SKE では風速ベクトルは定性的に良く合っていたものの、剥離域内部の逆流領域のサイズを過小 に評価していることがわかり、それに対して RSM、LES では実験と概ね一致していることがわかる。 5. まとめ 本稿では、比較的簡易な建物モデル内外の通風気流場を対象とした風洞実験と CFD 解析を行うこと で解析精度を検証した。ここでは CFD の乱流モデルに関して、渦粘性モデル、応力方程式モデル(RSM)、 Large Eddy Simulation(LES)と 3 段階に変化させ、それぞれにおいて一般に用いられるモデリング 手法を採用して建物内外の圧力及び風速を比較した。渦粘性モデルの中で今回用いた標準 k-e モデル (SKE)では、気流場は概ね実験と一致するものの、風速・静圧共に詳細な分布は実験結果とは異なる 傾向となる。RSM では概して傾向は実験結果と良く一致するものの、剥離域を過大評価する等、明ら かに実験とは異なる結果が得られる領域があることが示された。LES では実験と比較して差異が目立 つ箇所はなく、最も精度の良い解析結果が示された。本研究で採用したもの以外にも乱流モデルは多 数存在し、適切なモデルは解析対象の性状や着目する位置によっても異なると言える。また、実際の 解析 の際は、求められる解析精度と許容される計算負荷に応じて解析手法を選択する必要がある。 参考文献 1) 石原正雄:「建築換気設計」, 朝倉書店 , 1969. 2) 倉渕隆 , 大場正昭:「風洞実験と LES を併用した通風時の乱流構造解明に関する研究(3)通風の局所相似性に基づく通風量 と流入角の計算法に関する提案」, 日本建築学会大会学術講演梗概集 , D2, pp.645-646, 2001.9. 3) Murakami S., Kato S., Akabayashi S., Mizutani K., Kim Y.D. :「Wind Tunnel Test on Velocity-Pressure Field of CrossVentilation with Open Windows」, ASHRAE Transactions, Volume 97, Part1, pp.525-538, 1991. 4) 小林知広 , 相良和伸 , 山中俊夫 , 甲谷寿史 , 武田尚吾 , 西本真道:「通風時の建物周辺気流に関する風洞実験及び CFD 解析精 度の検証」日本建築学会環境系論文集 , 第 74 巻 , 第 638 号 , pp.481-488, 2009.4. 5) 小林知広 , 甲谷寿史 , 山中俊夫 , 相良和伸 , 桃井良尚 , 浅井香里:「Wind Tunnel Test and Numerical Simulation for Pressure and Velocity Field around Cross-Ventilated Building」日本建築学会環境系論文集 , 第 76 巻 , 662 号 , 2011.4 掲載予定 . 6) Willert C.E. and Gharib M.:「Digital Particle Image Velocimetry」Experiments in Fluids, Vol. 10, No. 4, pp.181-215, 1991. 7) Hart D.P.:「Super-Resolution PIV by Recursive Local-Correlation」, Journal of Visualization, Volume 3, Number 2, pp.187-194, 2000. 8) ポフ:「機械工学における空気力学実験法」, 朝倉書店 , 1969. 9) Wright M.A.:「The evaluation of a simplified form of presentation for five-hole spherical and hemispherical pitometer calibration data」, Journal of Physics E. Scientific Instruments, Volume 3, Number 5, pp. 356-362, 1970. 10) Launder B.E., Reece G.J., and Rodi W.:「Progress in the development of a Reynolds-stress turbulence closure」, Journal of Fluid Mechanics, Volume 68, pp.491-511, 1975. 11) We r n e r H . , We n g l e H . :「 L a rg e - E d d y S i m u l a t i o n o f Tr u b u l e n t F l o w O v e r a n d A r o u n d a C u b e i n a P l a t e C h a n n e l 」 , Proceedings of 8th Symposium on Turbulent Shear Flows, pp.155-168, 1993. 12) 大岡龍三 , 村上周三 , 持田灯:「<u i 'u j '> 輸送方程式中の圧力歪相関項 , wall reflection 項 , 乱流拡散項に関する各種モデルの 評価」, 日本建築学会計画系論文集 , 第 504 号 , pp.55-61, 1998.2. 13) 村上周三:「CFD による建築・都市の環境設計工学」, 東京大学出版会 , 2000. 14) Kato M., Launder B.E.:「The modelling of turbulent flow around stationary and vibrating square cylinder」, Proceedings of 9th Symposium on Turbulent Shear Flows, Paper, 10-4, 1-6, 1993. 15) 村上周三 , 持田灯 , 近藤宏二:「改良 k-e モデルによる 2 次元建物モデル周辺気流の数値計算」, 生産研究(東京大学生産技 術研究所所報), 第 47 巻 , 第 2 号 , pp.29-33, 1995. 謝辞 本研究を推進するにあたり、大変有益なご助言を頂いた大阪大学の山中俊夫教授、相良和伸教授、桃井良尚助教に深謝致しま す。実験及び解析に際して多大なるご尽力を頂いた浅井香里氏(東京大学大学院生、当時大阪大学卒論生)に感謝致します。また、 必要機器の貸し出しに加え、多くのご助言を賜った日本カノマックス株式会社の田中康恵氏、塩崎康弘氏に感謝致します。なお、 本研究の一部は日本学術振興会平成 21 年度科学研究費補助金(特別研究員奨励費 20-912、研究代表者:小林知広)によった。