計算工学講演会論文集 Vol. 21 (2016 年 5 月)

計算工学会

津波避難ビルの損傷評価のための津波伝播解析
Tsunami Propagation Analysis for Damage Estimation of Evacuation Building
田中聖三 1) 荻野弘明 2) 磯部大吾郎 3)
Seizo Tanaka, Hiroaki Ogino and Daigoro Isobe
1)2)

博 (工) 筑波大学 システム情報系 助教（〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: stanaka@kz.tsukuba.ac.jp）

筑波大学大学院システム情報工学研究科（〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: s1420873@u.tsukuba.ac.jp）
3)

博 (工) 筑波大学 システム情報系 教授（〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: isobe@kz.tsukuba.ac.jp）

This paper presents a tsunami propagation analysis system for damage estimation of evacuation building.
For the problem of the wave propagation and estimating the wave force to the building, we have developed
the three dimensional free surface flow analysis code based on the volume of fluid (VOF) method. The
stabilized finite element method and OpenMP/MPI hybrid parallelization were implemented. As for thenumerical example, the estimation of tsunami wave force on the evacuation buildings were carried out. The
drag coeﬃcient obtained by the present method were compared on the several inflow condtions and building
shapes.
Key Words : Tsunami propagation analysis, evacuation building, VOF method

1. はじめに
地震発生時に励起される可能性のある津波災害は，沿岸域に多大な人的・物的被害が引き起こされることは記憶に新しい．津波発生時に沿岸域でとるべき行動は，原則として高台へ避難することであるが，高台から遠く離れた平地部においては，津波避難ビルが建設されている．この津波避難ビルは，津波波力に対して十分な堅牢性を有している必要があり，想定しうる様々な津波波力に対して，設計時に十分検討される必要がある．
これらの多数の津波シナリオを取り扱う検討方法として，数値解析による検討方法は，有効な手法であるといえる．十分な堅牢性を評価するためには，建物の部材の終局状態における損傷を検討しなければならず，そのためには，部材細部に作用する局所的な波力も考慮しなければならない．合わせて，建物の詳細な形状を表現する必要があり，これらの要求を満足する津波伝播解析に必要なモデルの自由度は多大なものとなる．
これらの大規模問題に対して，大規模並列計算機を用いて解析を行うことが可能となってきているが，設計時には部材寸法の変更などによる構造部材の形状変化があり，その度に津波入力の全てのケースを再計算して波力を推定することは，やはり計算コストは莫大なものとなる．
そこで本研究では，津波避難ビルの設計時における損傷評価のための津波伝播解析手法を構築するものである．まず，構造物の形状を詳細に表現し，これらに作用する波力を高精度に評価することが可能な解析手法の構築を行った．大規模問題に対して，OpenMP/MPIによる Hybrid 並列計算法を導入した．また，時々刻々の自由表面流れを解析するために，安定化有限要素法に基づく Volume of Fluid (VOF) 法を実装した．数値解析例として，津波避難ビルに対する津波伝播解析を取り上げ，津波避難ビルの形状による波力の変化について検討を行った．また，津波避難ビルの形状による抗力係数の変化についても検討を行い，設計時に用いる津波波力の推定式の有効性の検討を行った．

2. 自由表面流れ解析手法(1) 流体の基礎方程式
自由表面流れのモデルとして密度と粘性係数が異なる 2 種類の流体を考える．下添字 1 と 2 を使って密度と粘性係数を ρ1,2 と µ1,2 とする．この 2 種類の流体が混ざった解析領域内において，密度と粘性係数を場所の関数とする．すなわち，

ρ = cρ1 + (1 − c) ρ2 ,(1)µ = cµ1 + (1 − c) µ2 ,(2)ここで，c は場所の関数である界面関数であり，1 の流体内では c = 1，2 の流体内では c = 0 となる．流れに伴い界面関数は時間変化するが，この時間変化は次のALE 標記による移流方程式に従うことを仮定する．

∂c+ ū · ∇c = 0∂ton Ωt(3)

ここで u は流速，メッシュ移動速度を û とすると ū = u−ûである．流速 u は次の非圧縮性粘性流体の基礎方程式の解となる．

()∂uρ+ ū · ∇u − f − ∇ · σ(u, p) = 0∂t∇·u=0on Ωt(4)on Ωt(5)ここで，p は圧力，f は外力，そして σ は流体の応力である．

(2) 有限要素法による定式化
式 (3) に対して SUPG 法に基づいた安定化有限要素法を適用すると以下の式が得られる．
( h)∫∂cch∗+ ūh · ∇ch dΩ∂tΩ)
nel ∫∑() ( ∂chhhhh+
τcS ū · ∇c∗+ ū · ∇c dΩ∂te=1 Ωenel ∫∑+(6)
τcDC ∇ch∗ · ∇ch dΩ = 0e=1Ωe

ここで τcS と τcDC は安定化パラメータであり，以下のように定義される [3]．
( )2
−1/2
 22||ū|| τcS = +,(7)∆th τcDC =h||ū||,2(8)ここで，h は要素長である．
非圧縮性粘性流体の基礎方程式 (4),(5) に対して
SUPG/PSPG 法に基づく安定化有限要素法を適用すると，以下の弱形式を得る：( h)∫∂uhhhw ·ρ+ ū · ∇u − f dΩ∂tΩt∫∫hh h+ε(w ) : σ(u , p ) dΩ +qh ∇ · uh dΩΩt{図–1 解析モデル (上:平面図，下:側面図)Ωt}1τS ūh · ∇wh + τP ∇qheρ
e=1 Ωt{ ( h)}∂u· ρ+ ūh · ∇uh − f − ∇ · σ(uh , ph ) dΩ∂t∫nel ∫∑hh+τDC ∇ · w ρ∇ · u dΩ =wh · hh dΓ, (9)+nel ∫∑e=1Ωet(Γt )h

ここで τS , τP , τDC は安定化パラメータ [4] である．h は運動方程式 (4) の Numann 境界条件である．
界面関数 c，流速 u, 移流速度 ū および圧力 p の時間積分には，Adams-Bashforth 近似による Crank-Nicolson法を用いる．この結果，次式を得る．
)
1( hch =cn+1 + chn ,(10)2)
1( hun+1 + uhn ,(11)uh =(2)13ūh = uhn − uhn−1 − ûhn ,(12)22ph = phn+1 .

(13)以上の空間・時間方向の離散化によって導かれた連立一次方程式の解法には，GP-BiCG 法を適用する．これらの計算を行うプログラムは，計算時間と計算機記憶容量削減のため，領域分割法に基づく並列計算手法 (MPI により通信) と，Multi-Core Processor に対応した OpenMPを用いたスレッド並列を組み合わせた hybrid 並列計算手法を導入する．

図–2 避難ビルの形状 (左:開口部無し，右:開口部有り)

3. 数値解析例(1) 解析モデル
数値解析例として，津波避難ビルに作用する津波波力の評価問題を取り上げる．解析領域としては，高さ:80m,幅:180m, 奥行き:240m を設定した．避難ビルは，10 階建てを想定し，高さ:40.15m, 幅:32m, 奥行き:18.6m の外形を設定した．解析領域の平面図，側面図を図-1 に示す．また，避難ビルの形状は図-2 に示されるような開口部の無いモデルと開口部の有るモデルの２種類を用意し，津波波力の比較を行う．
入射波の条件として，表-1 に示されるように，水深 8mを一定とし，流入流速を 6,9,12m/s と変化させる 3 Caseを設定し比較を行う．なお，各 Case の流入流速と水深から算出されるフルード数 (Fr) は，0.67〜1.35 となっている．流入流体は海水を仮定し，密度 ρ1 = 1027kg/m3 ，粘性係数 µ1 = 1.0 × 10−3 Pa · s とした．また，気体は表–1 入射波条件

Case 1Case 2
Case 3

水深 (m)
888

流入流速 (m/s)
6.0
9.012.0

フルード数 Fr0.67
1.00
1.35

Case 1

Time= 8.0s

Time=10.0s

Case 2

Case 3
図–4 時刻 10s における避難ビル前面の圧力分布 (左:開口部無し，右:開口部有り)

Time=12.0s
図–3 各時刻における水面形状 (Case2, 左:開口部無し，右:開口部有り)

空気を仮定し，密度 ρ2 = 1.293kg/m3 ，粘性係数 µ2 =1.8 × 10−5 Pa · s とした．
解析に用いた有限要素は 3 次元四面体要素であり，最小メッシュサイズは，流入部の水深以下や避難ビル周辺で約 50cm とした．解析モデルの総節点と総要素数
は，開口部無しモデルで約 300 万節点と約 1700 万要素，開口部有りモデルがそれぞれ，約 370 万と約 2100 万である．
境界条件として，底面，側壁面および避難ビル表面
では slip 条件を課している．今回，液体の流出境界は設定しておらず，流入された海水が後部境界に到達するまでを解析時間とした．

(2) 解析結果
解析結果として，図-3 に Case2 の各時刻における水面形状を示す．左側に開口部無し，右側に開口部有りの結果を示す．これらの図より，開口部無しの場合は，建物壁面に沿って津波の波高が大きくなっているのに対し，開口部有りの場合では，波高が低減されていることがわかる．また，時刻 10s でのビル前面の圧力分布の比較を図-4 に示す．この図は，開口部の有無，入射波条件により比較するものである．入射波条件のフルード数が上がるにつれて，ビルに作用する圧力が上がっていることがわかる．また，開口部無しの場合に比べて，開口部有りの場合の方が低減されているのは，入射波条件によらず同様の傾向であった．
ここで，解析により得られたビル前面に作用する波力の時刻歴を図-5 に示す．図中の実線が解析により得られた波力である．今，簡易的に波力を算定するために，推定波力を以下のように静水圧による波力と抗力の和で算定することを考える．
Fwave = FS + F D ,∫FS = γρgB z dz

(14)(15)h

FD =1γρC D U 2 A2(16)904Case1 w/o Aperture (computed)
Case2 w/o Aperture (computed)
Case3 w/o Aperture (computed)
Case1 w/o Aperture (estimated)
Case2 w/o Aperture (estimated)
Case3 w/o Aperture (estimated)

Wave Force [MN]

7060

Case1 w/o Aperture (computed)
Case2 w/o Aperture (computed)
Case3 w/o Aperture (computed)

3.5
3Drag Coe cient

80

50
4030

2.5
21.5
120
0.5

1000051015

2025300510Time [s]902025304
Case1 w/ Aperture (computed)
Case2 w/ Aperture (computed)
Case3 w/ Aperture (computed)
Case1 w/ Aperture (estimated)
Case2 w/ Aperture (estimated)
Case3 w/ Aperture (estimated)

7060

Case1 w/ Aperture (computed)
Case2 w/ Aperture (computed)
Case3 w/ Aperture (computed)

3.5
3Drag Coe cient

80

Wave Force [MN]

15Time [s]5040302.521.5
1200.5

10000510152025300

Time [s]図–5 避難ビル前面に作用する波力の時刻歴 (上:開口部無し，下:開口部有り)

γ は開口低減係数であり，開口部無しの場合が γ = 1.0，開口部有りの場合で γ = 0.52 となる．また， B は受圧面の幅，U は代表流速，A は受圧面の面積である．特に C D は抗力係数であり，建物の形状に依存する係数である．抗力係数以外のパラメータは入射波条件などにより設定可能であるが，任意形状の建物に対する抗力係数を設定することは難しい．ここで，解析により得られた図-5 の波力から，式 (16) で表される静水圧成分を減じ，式 (16) を用いて，抗力係数を算定することを考える．図-6 に解析により得られた波力から算出された抗力係数の時刻歴を示す．図-5 では，得られた波力は入射条件により大きさが変化しているが，図-6 に示される抗力係数の時刻歴では，その大きさは入射条件の Case によらず，ほぼ同程度の値となっている．抗力係数は建物の形状により決定されるため，これは当然のことである．得られた抗力係数により算定された推定波力を図-5 に破線で示しているが，この図より津波波力を推定できていることがわかる．これは，一度数値解析により，波力を求め，その時系列から抗力係数を算出すれば，簡便に津波波力が推定可能であることを示唆している．このことは，入射波条件が複数有る
場合，3 次元流体解析を行う回数を削減することが可能であるといえ，避難ビルの設計時の構造部材の損傷評価におけるコストの削減に寄与すると考えられる．
4. おわりに
本論文では，津波避難ビルの損傷評価のための津波伝播解析手法の構築を行い，避難ビルの開口部の有無による波力の変化について検討を行った．また，計算された波力より抗力係数を算出し，算出された抗力係数を用いて，推定波力を算定した．算定された推定波51015202530Time [s]図–6 算定された抗力係数の時刻歴 (上:開口部無し，下:開口部有り)

力は津波波力を評価することが可能であった．これにより，避難ビルの設計時の損傷評価のコストが削減されると考えられる．
今後の課題として，津波波力と構造物の連成解析手法の開発を行う予定である．
参考文献[1] 津 波 避 難 ビ ル 等 に 関 わ る ガ
イ ド ラ イ ン,
内 閣 府:http://www.bousai.go.jp/oshirase/h17/tsunami hinan.html(2016 年 3 月 20 日閲覧)
[2] 津 波 避 難 ビ ル 等 の 構 造 上 の 要 件 の 解 説,
日 本 建 築 防 災 協 会:
http://www.kenchikubosai.or.jp/files/2013/11/12 tsunami01.pdf(2016年 3 月 20 日閲覧)
[3] Aliabadi, S. and Tezduyar, T. E. :Stabilized-finiteelement/interface-capturing technique for parallel
computation of unsteady flows with interfaces, Compute Methods in Applied Mechanics and Engineering,
Vol.190, pp.243-261, 2000.
[4] Aliabadi, S., Abedi, J. and Zellars, B. :Parallel finite
element simulation of mooring forces on floating objects, Int. J. Numer. Meth. Fluids, Vol.41, pp.809-822,
2003.
[5] 田中聖三, Sun, F., 堀宗朗, 市村強, Wijerathne,
M.L.L.: 動的津波荷重による構造物の破壊解析の基礎的研究, 土木学会論文集 A1(構造・地震工学),
Vol.69, No.4, I 903-I 908, 2013.