2.2 テーマ【II-2】 市街地の風環境解析法 2.2.1 はじめに 市街地上空および建物周辺の風況に関する数値流体解析法を構築する．解析対象空間は多 数の建物のスカイラインを境界とする複雑な形状である．本研究では，有限要素法に基づく 数値流体解析を行うことを方針とし，(1) 解析モデルである有限要素メッシュの作成法，(2) 乱流変数の自然対数をとる変数を用いることで数値的安定性を高める logarithmic form によ る乱流解析法，の２つの課題を進めている． 2.2.2 市街地上空の有限要素メッシュ作成法 a) メッシュ作成法の基本的方針 都市上空の風況解析のための解析メッシュはさまざまな方法で作成されている．構造格子 を基本とする差分法による解析では建造物周りの基本メッシュを重合する方法が多く用いら れている[2-2-1,2-2-2]．有限体積法や有限要素法は非構造格子による解析が行われるため， 複雑形状まわりの解析対象領域をメッシュ分割する汎用的な手法が適用される傾向にある [2-2-3, 2-2-4]．しかし汎用的な手法は万能の手法ではないので，解析対象が特定の工学的 問題に限定されている場合に，その対象に応じてチューニングされたメッシュ生成法が開発 されて利用されている[2-2-5]． 本研究の解析対象は，市街地上空の建造物まわりの空間である．また本研究プロジェクト II のテーマ【II-1】で研究される都市の立体構造に関する高度データを積極的に活用できる ようにしたい．これらの状況から，本研究の解析メッシュ作成の基本方針として，まず，図 2-2-1(a)に示すように，地図情報に基づいて地表面における建物周辺空間を２次元的にメッ シュ分割する（これを底面メッシュと呼ぶ） ．次に図 2-2-1(b)に示すように底面メッシュイ メージを上方に平行移動して有限要素を必要な高度まで積み上げる，という２段階の手順の メッシュ作成法を開発することとした． 本年度は，任意形状断面の複数の建物まわりの空間に底面メッシュを作成する方法の作成 を試みた．底面メッシュ作成のアルゴリズムは，清水ら[2-2-6]が開発した，実空間の幾何情 報をいったん写像空間に配置する方法を用いた．写像空間から解析メッシュを作成する段階 では，有限要素のひずみが少なくなるように支配方程式の選択を検討し，さらに要素のひず みレベルに応じたメッシュ分割の補正法を工夫した． (a) (b) 図 2-2-1 風況解析用解析メッシュ作成のコンセプト (a) 底面メッシュ，(b) 上方への要素の積み上げ 1 / 12 b) メッシュ作成アルゴリズム 底面のメッシュ分割アルゴリズムの手順は以下のとおりである． 図 2-2-2 アルゴリズムの手順 1) ２次元の実空間 (x, y ) における解析領域内の建築物の外形を線分のみで構成する．これ を幾何モデルとよぶ [図 2-2-2(a)]． 2) 各線分を直交座標軸 (ξ ,η ) に平行な線分に割り当て，長さを決定する．このモデルを認 識モデルとよぶ [図 2-2-2(b)]． 3) 認識モデルを直交格子状に要素分割する．このモデルを写像モデルとよぶ [図 2-2-2(c)]． 4) 写像モデルの節点と実空間の節点との間の写像関係に関する方程式を解くことで実空間 での節点座標値を求め，解析メッシュを作成する [図 2-2-2(d)]． c) 認識モデルにおける辺の方向の決定 幾何モデル[図 2-2-2(a)]の各辺を写像空間での直交座 標軸 (ξ ,η ) のどちらかの方向に割り当てる．その可能性 ( ) を示す適応度 Pξ , Pη をメンバシップ関数を用いて表す． 図 2-2-3 に示すように，幾何モデルの各辺と x 軸との成 す角度を θ x とし，図 2-2-4(a) に示すメンバシップ関数 より Pξ の初期値を求める．また，対象とする辺と y 軸と なす角度 θ y により決まる適応度を Pη と表す．なお Pη = 1 − Pξ である． 2 / 12 図 2-2-3 辺の角度と挟角 図 2-2-4 角度の決定に用いるメンバシップ関数 続いて図 2-2-3 に示す隣接する２辺 AB の挟角 θα から図 2-2-4(b) に示すメンバシップ関数 より Pα を求める． Pα は２辺 A, B が同一方向に写像される可能性が高いか，２辺 A, B が異 なる方向に写像される可能性が高いかを表す値である． 各辺の適応度の初期値を決定したら， 次式により適応度を修正する． (i) θα ≥ 0.5 の場合 ΔPξ′B = β (Pα − 0.5)PξA PξB ΔPξ′′B = − β (Pα − 0.5)PηA PηB (ii) θα < 0.5 の場合 ΔPξ′B = β (0.5 − Pα )PηA PξB ΔPξ′′B = − β (0.5 − Pα )PξA PηB (1a) (1b) (2a) (2b) β は修正速度を決定する係数で 0.05 としている．これより，辺 A の影響による辺 B の適 応度の修正を次のように行う． PξB (i +1) = PξB + ΔPξ′B + ΔPξ′′B (3) すべての辺の Pξ の値が 0 か 1 に収束するまで式(1)∼(3)の計算を繰り返す． Pξ = 1 の場合 は ξ 軸に平行な辺に写像され， Pξ = 0 ならばη 軸に平行な辺に写像される． d) 認識モデルの辺の長さの決定 幾何モデルの各辺を認識モデルの ξ 軸あるいはη 軸のいずれかに平行になるように方向を 決めた段階では，図 2-2-5(a) のように図形が閉じていない．これを認識モデルにおいて図 2-2-5(b)のように閉じた図形にするために，認識モデルにおける各辺の辺長を次の手順によっ て決定する． (a) (b) 図 2-2-5 (a) 直交方向に辺の方向を決定した段階の幾何モデル，(b) 認識モデル． 3 / 12 幾何モデルの各方向別の辺の長さの和 S i を次式(4)により求める． Ni S i = ∑ sij (i = 1 − 4) (4) j =1 ここで，sij は幾何モデルでの線分の長さを示す．添え字 i は ξ 軸の正負の方向およびη 軸の 正負の方向の４つの方向を表している．また N i は i 方向に属する辺の数である． 認識モデルにおける４方向の辺の和 Ti (i = 1 − 4 ) を，次式(5a, b)のように幾何モデルにおけ る ξ 軸の正負方向の辺長和 S1 , S 3 の平均値，およびη 軸の正負方向の辺長和 S 2 , S 4 の平均値 になるようにする． T1 = T3 = (S1 + S 3 ) 2 , T2 = T4 = (S 2 + S 4 ) 2 (5a, b) これより認識モデルにおける辺の長さ t ij を次のように決定する． t ij = (Ti S i ) sij (6) なお，内角が全て 180°以内の図形（凹部がない図形）の場合は，以上に述べた手順では なく，もっと単純な方法をで認識モデルを作成すればよい．すなわち各辺のもともとの適応 度 Pξ が 0 か 1 のいずれに近いかを調べ，近い方の 0 あるいは 1 を選択して認識モデルを作成 する． 図 2-2-6 は円を近似した正 20 角形の幾何モデル(a)と認識モデル(b)である． 図 2-2-6 (a) 円の幾何モデル，(b) 認識モデル e) 解析メッシュの生成 認識モデルを取り除いた領域を直交格子状に要素分割して写像モデルを作成したのち，構 造物の辺上に位置する節点座標を定め，これを境界条件として写像モデルと実空間の写像関 係 x = x(ξ ,η ), y = y (ξ ,η ) を支配する方程式の境界値問題を有限要素法によって解いて解析 メッシュを作成する（図 2-2-7） ．支配方程式には Laplace 方程式と弾性体の方程式を用い， 4 / 12 得られるメッシュを比較検討した． ２種類の支配方程式は以下のように表される． (i) Laplace 方程式 ∂2x ∂ξ 2 + ∂2x ∂η 2 = 0, ∂2 y ∂ξ 2 + ∂2 y ∂η 2 =0 (7) 図 2-2-7 解析メッシュの境界条件 (ii) 弾性体の支配方程式 ∂ ⎧ E ⎛ ∂x ∂y ⎞⎫ ∂ ⎧ E ⎛ ∂x ∂y ⎞⎫ ⎜⎜ ⎟⎬ + ⎜ ⎟⎬ = 0 +ν + ⎨ ⎨ 2 ∂ξ ⎩1 − ν ⎝ ∂ξ ∂η ⎟⎠⎭ ∂η ⎩ 2(1 + ν ) ⎜⎝ ∂η ∂ξ ⎟⎠⎭ ∂ ⎧ E ⎛ ∂x ∂y ⎞⎫ ∂ ⎧ E ⎛ ∂x ∂y ⎞⎫ ⎜ν ⎟⎬ = 0 ⎜ ⎟⎬ + + + ⎨ ⎨ ∂ξ ⎩ 2(1 + ν ) ⎜⎝ ∂η ∂ξ ⎟⎠⎭ ∂η ⎩1 − ν 2 ⎜⎝ ∂ξ ∂η ⎟⎠⎭ (8a, b) ここで， E はヤング率，ν はポアソン比である． 図 2-2-8(a)は Laplace 方程式，図 2-2-8(b)は弾性体の方程式( E = 1.0, ν = 0.45 )を解いて作成 した解析メッシュを示す．Laplace 方程式より作成したメッシュはメッシュが物体内部に食い 込んでしまっているのに対し．弾性体の方程式の結果はメッシュが物体内部に食い込むこと がなかった．また，ポアソン比が 0.5 に近いほど非圧縮性の変形となってより良好なメッシ ュを得ることが予想される．図 2-2-8 (c)はポアソン比ν = 0.499 としたときの解析メッシュで ある．図 2-2-8(b)と比べてゆがみが低減されている結果になっている． (a) Laplace 方程式 (b) 弾性体方程式（ポアソン比 0.45） (c) 弾性体方程式（ポアソン比 0.499） (d) 弾性体方程式（重み操作後） 図 2-2-8 解析メッシュと拡大図 5 / 12 f) ひずみの平滑化 物体の角付近の要素は，写像前の矩形からのひずみが大きくなる．要素のひずみが領域全 体で平滑化するように，以下の操作を行った[2-2-7]．弾性体の方程式(8a, b) を解くことで作 成される解析メッシュの各要素のひずみ εx = ∂y ∂x ∂x ∂y ,εy = , γ xy = + ∂ξ ∂η ∂ξ ∂η (9a-c) を算出し，式(10)で表されるひずみのノルム ε e を要素ごとに求める． 2 ε e = ε x2 + ε y2 + γ xy (10) ひずみのノルム ε e が最大となる要素を基準として各要素の重み C e を式(11)より求める． Ce = 1 + α εe (11) ε e,max ここで α は重みの値が過大にならないように設定した係数で，本研究では α = 20 とした． 各要素のヤング係数を C e E として弾性体の方程式(8a, b)を再度解く．その結果が図 2-2-8(d) である．この方法によってメッシュのゆがみをさらに改善することができる． g) 駿河台キャンパス地区の底面メッシュ 以上の方法を用いて，実在市街地の底 面メッシュを作成した．対象領域は，図 2-2-9 に示す駿河台キャンパス地区とし た． 図 2-2-10(a) が解析対象領域の幾何モ デル，図 2-2-10(b)が認識モデルである． 解析対象領域内の物体数は 24 である． 物 体が複数あること，および物体が解析対 図 2-2-9 駿河台キャンパスの周辺 （長方形内が解析メッシュの作成対象領域） 象の境界に接していること，のいずれも問題なく本手法が適用できることが示されている． 図 2-2-11 が認識モデルの物体外部空間をメッシュ分割した写像モデルである．総節点数 9411，総要素数 8452 である．この写像モデルを，弾性体方程式でポアソン比 0.499 として有 限要素法によって解いて作成した解析メッシュが図 2-2-12(a)である．また，図 2-2-12(a)の各 要素のゆがみから重みを評価して，再度弾性体方程式を解いて得られた解析メッシュが図 2-2-12(b)である．あまり大きな変化はないが，鋭角部付近のメッシュが改善されている． 6 / 12 (a) (b) 図 2-2-10 駿河台キャンパス地区の (a)幾何モデル，(b) 認識モデル 図 2-2-11 駿河台キャンパスの写像モデル (a) (b) 図 2-2-12 駿河台キャンパスの解析メッシュ (a) 重み操作前，(b) 重み操作後 h) 解析メッシュ作成法のまとめ 市街地の風況解析に用いる解析メッシュのうち，底面メッシュを作成する方法を，清水ら の方法に独自の工夫を加えて作成した．特に，写像モデルから解析メッシュを作成するにあ たって，ポアソン比が 0.5 に近い弾性体方程式を使用すること，およびいったん作成した解 析メッシュの各要素のひずみから要素ごとの重みを評価して重み付ヤング率を用いて計 s ん しなおすことにより，ゆがみが少ないより良好なメッシュを得ることができる．本方法が， さまざまな底面形状を有する複数の建造物を含む解析対象領域に適用できることを駿河台キ ャンパスを対象とする解析メッシュを作成することによって確認した． 7 / 12 2.2.3 有限要素法による風況解析の構築 a) 風況解析法 市街地上空および建物周辺の風の解析は，３次元の解析対象領域における乱流解析を必要 とする．本研究では，工学的な問題に多用されている乱流モデルによる解析法を構築する． 乱流モデルを用いる解析は，多くの変数が相互に関連する非線形な過程であり，数値的に不 安定にならないような解析過程を必要とする．本研究では，プロジェクト代表者が非流線形 物体まわりの乱流解析のために具体化してきた Logarithmic form を用いた k-ε 系モデルによ る解析法 [2-2-8] を具体化する方針である．本年度は，これまでの２次元乱流解析法を３次 元化する課題に着手し，推進してきた．以下に，２次元解析の範囲であるが，本解析法の要 点を記す． b) Logarithmic form k-ε 系モデルは，乱流エネルギー k ，エネルギー散逸率 ε ，時間スケールτ などを変数と して構成される．これらは物理的に正の値のみをとる変数であるが，数値計算の過程で負の 値をとる場合がある．このとき，負の値となった変数を，適当に設定した小さな正の値へ強 制的に修正する，clipping という物理的に意味を持たない操作が一般的に行われる．この clipping を回避するために Ilinca ら[2-2-9]によって Logarithmic form が提案されている． Logarithmic form は乱流量を表す変数の自然対数をとった新たな変数を定義し，その変数に関 する方程式を解く手法である． 乱流エネルギー k およびエネルギー散逸率 ε の自然対数をとり，式(12a, b) に示す新たな変 数 K , E を定義する．この定義から K , E が負になっても k , ε は正のままであることが保証 される． K = ln (k ), E = ln (ε ) (12a, b) 本来の変数である k , ε の微係数は，新しい変数 K , E の微係数との間に次の関係がある． ∂K 1 ∂k , = ∂t k ∂t ∂K 1 ∂k , = ∂x j k ∂x j ∂E 1 ∂ε , = ∂t ε ∂t ∂E 1 ∂ε = ∂x j ε ∂x j (13a-c) この関係を標準 k-ε モデルの k , ε の輸送方程式に適用すると以下の方程式が導かれる． ∂K ∂K ∂ +U j = ∂t ∂x j ∂x j ⎛ ν ⎞⎛ ∂K + ⎜⎜ν + t ⎟⎟⎜ σ k ⎠⎜⎝ ∂x j ⎝ ⎡⎛ ν t ⎞ ∂K ⎤ ⎟⎟ ⎢⎜⎜ν + ⎥ σ ⎥ k ⎠ ∂x j ⎦ ⎣⎢⎝ 2 ⎞ ⎛ ∂U j ⎟ + e − Kν t ∂ ⎜ ∂U i + ⎟ ∂x j ⎜⎝ ∂x j ∂xi ⎠ 8 / 12 K ⎞ ⎟ − Cμ e ⎟ νt ⎠ (14) ∂E ∂E ∂ +U j = ∂t ∂x j ∂x j ⎛ ν ⎞⎛ ∂E + ⎜⎜ν + t ⎟⎟⎜ σ k ⎠⎜⎝ ∂x j ⎝ ⎡⎛ ν t ⎞ ∂E ⎤ ⎟⎟ ⎢⎜⎜ν + ⎥ σ ⎥ k ⎠ ∂x j ⎦ ⎣⎢⎝ 2 ⎞ ⎛ ∂U j ⎟ + Cε 1e − Kν t ∂ ⎜ ∂U i + ⎟ ∂x j ⎜⎝ ∂x j ∂xi ⎠ ⎞ ⎟ − Cε 2 e E − K ⎟ ⎠ (15) 流速，圧力は対数をとらず，運動方程式は変更しない． c) 解析スキーム 運動方程式および式(14), (15)を有限要素法により離散化した．空間は双線形四角形要素で 分割し，流速 U i ，対数をとった乱流エネルギー K ，エネルギー散逸率 E は双線形分布，圧 力 P ，渦粘性係数ν t は要素内一定分布の補間関数を用いた．また，安定化手法として SUPG 法を適用し，時間方向の離散化には Predictor-Corrector 法を用いた． d) バックステップ流れの解析 バックステップ流れの解析条件を 図 2-2-13 に示す．笠木ら[2-2-10]の実 験を比較対象とした．レイノルズ数 は 5,500 である．固体境界である上 下境界の境界条件には壁関数を用い た．図 2-2-14(a)は乱流エネルギー k およびエネルギー散逸率 ε の自然対 数をらない通常の解析法（以下， Γ1 : U max = 80.37 cm/s , I u = 2 % Normal form と称する） の定常状態の Γ2 : σ ij ⋅ n j = ∂k ∂ni = ∂ε ∂n j = 0 流線，図 2-2-14(b) には Logarithmic Γ3 : 壁関数 form で会席して得られた定常状態 の流線を示す．ステップ高さで無次 図 2-2-13 バックステップ流れの解析条件 元化した渦の再付着距離は， 実験値が 6.5 であるのに対し， Normal form が 5.7， Logarithmic form が 6.1 であり，標準 k-ε モデルを用いた結果としてはどちらも良好であるといえる． 図 2-2-14 定常状態の流線：(a) Normal form, (b) Logarithmic form 9 / 12 (a) (b) 図 2-2-15 過渡的状況における乱流量の分布：(a) Normal form の渦粘性係数分布 ［上］および矢印の位置の乱流量の水平方向分布［下］ ，(b) Logarithmic form の 渦粘性係数分布［上］および矢印の位置の乱流量の水平方向分布［下］ ． e) Logarithmic form の効果 ( 図 2-2-15 の上段には，流れが定常状態となる途中の渦粘性係数ν t = C μ k 2 ε , C μ = 0.09 ) の分布図を示した．下段には，上段の図に矢印で示したステップ高さにおけるステップ後方 の乱流量を示した．なお k , ε は対数目盛りでプロットしている．Normal form では図の中ほ どに過大なν t が生じているが，Logarithmic form ではそれが生じていない．このν t はステッ プの先端から生じ， k , ε の発達に伴い移流していく．Logarithmic form は全体で k , ε のバラ ンスがとれているのに対し，Normal form はステップからの水平距離が 120∼150cm の付近で k , ε のバランスが崩れており，これによりν t が過大に評価された．このように Normal form では k , ε のバランスが崩れ，過大なν t が生じたが，Logarithmic form ではこれが解消される ことがわかった．この要因としては，clipping が回避されたことや， k , ε の自然対数である K , E の方程式を解いているため，変数の勾配が全体で低下したことが考えられる． Logarithmic form が数値的な安定性を高め，特に非定常な流れのときに効果が顕著であること が示された． 2.2.4 まとめ 以上， 平成 16 年度は市街地の風況解析に必要となる解析メッシュ作成法および数値流体解 析法の２つの課題に取り組んだ．それぞれの課題について独自の観点を含む数値的手法を検 討し， その基本的な有効性を確認した． 平成 17 年度は， ２つの手法を 3 次元解析用に展開し， 引き続き改良・工夫を加えていく予定である． 10 / 12 参考文献 [2-2-1] 中山浩成，田村哲郎，奥田泰雄：実在都市での危険性物質の大気拡散に関する LES， 第 18 回風工学シンポジウム論文集，pp.121-126, 2004. [2-2-2] 佐々木澄，持田灯，吉野博，岩田達明：仙台市中心部の街路樹がストリートキャニオ ン内の風環境と空気環境に及ぼす影響の解析，第 18 回風工学シンポジウム論文集， pp.99 -104, 2004. [2-2-3] 浜田秀敬，桜井紘己，樫山和男，谷口健男：CAD/GIS による都市環境シミュレーシ ョンのための形状モデリング手法の構築，計算工学講演会論文集，Vol.9, pp.765-768, 2004. [2-2-4] 樫山和男：国際フォーラム「土木工学における計算力学手法の新展開」の企画調査， 平成 16 年度科学研究費補助金（基板研究(C)（企画調査） ）研究成果報告書，第 8 章， 2005. [2-2-5] 真下啓治，森脇清明，大林宏司，谷口健男：回路基板の幾何学的特性を活用した３次 元有限要素生成法の提案，計算工学講演会論文集，Vol.9, pp.781-784, 2004. [2-2-6] 清水ひろみ，高橋宏明，千葉矩正，山下禎文：3 次元シェル構造物の有限要素自動生 成，日本機械学会論文集（Ａ編） ，Vol.59, No.565, 1993. [2-2-7] 野村卓史，西村拓：ALE 有限要素流れ解析のためのメッシュ変形パターン生成法， 土木学会論文集 No.455/I-21, pp.65-74, 1992. [2-2-8] 長谷部寛，野村 卓史：k-εモデルにおける Logarithmic form の有効性の検討と非定 常流れへの適用，構造工学論文集 Vol.51A, pp.921-932, 2005. [2-2-9] F.Ilinca, D.Pelletier: Positivity preserving formulations for the k‒εmodel of turbulence，AIAA Journal , Vol.36, No.1, pp.44 – 50, 1998. [2-2-10] N. Kasagi, S.Kawara and A.Matsunaga: Turbulence measurement in a separated and reattaching flow over a backward - facing step with the aid of three-dimensional particle tracking velocimetry, The Symposium of the Society of Instrument and Control Engineers, 1991. [2-2-11] 加藤真志：修正２方程式乱流モデルによる角柱の基本空力特性に関する研究，名古 屋大学学位論文, 1997. 発表論文リスト [2-2-1] 長谷部寛，野村卓史：k-εモデルにおける Logarithmic form の有効性の検討と非定常 流れへの適用，構造工学論文集 Vol.51A, pp.921-932, 2005. [2-2-2] 小野大輔，長谷部寛，野村卓史：都市の風況解析に向けた三次元メッシュ生成アルゴ リズムの構築，第 32 回土木学会関東支部技術研究発表会講演概要集，I-21, 2005. [2-2-3] 小野大輔，長谷部寛，野村卓史：都市の風況解析に向けた任意形状の建物を含む領域 の解析メッシュ作成法， 第 60 回土木学会年次学術講演会講演概要集 （発表予定） ， 2005. 11 / 12 テーマ【II-2】 市街地の風環境解析法 発表論文抜粋（別紙参照） 12 / 12
© Copyright 2024 ExpyDoc
