■VOF 法による津波解析の手法
１ 自由表面流れの数値解法
自由表面を有する津波の 3 次元 CFD（計算流体力学）シミュレーションを実現する数値解析手法について概観する。

波の表面より下は水であり、単純に非圧縮性流体に対する CFD技法（陽解法）が適用できる。 津波のシミュレーションで肝要なのは自由表面の形状を正確に表現する、いわゆる自由表面捕捉(Free Surface Capturing)の技法が重要である。

１－１ 高さ関数法[1]表面位置を固定壁からの 1 価関数（高さ関数：Height Function）で表す。
ｙ＝Ｈ（ｘ、ｔ），∂H/∂ｔ＋∂H/∂ｘ＝０
利点は、計算方法が単純で処理が早く必要メモリーが少なくてすむ事、多次元への拡張が容易な事が挙げられる。

欠点としては表面の勾配が格子対角線Hの傾きを越えられない事、表面が巻き込まれる問題には適用できない事で津波の解析には不向きである。

図 1 高さ関数法の概念

１－２ SMAC(Simplified Marker And Cell )法[2]
流体が占有する領域をマーカー粒子で表現する手法である。 表面はマーカー粒子を含み、かつ隣接するセルが粒子を含まない場所と定義格子点の流速を内挿してマーカー粒子を移動させる。 利点は表面形状に制約がなく複数の流体領域が扱える事、表面の合体／分離に適用できる事である。 欠点は、マーカー粒子の定義に多くのメモリーを必要とする事、個々の粒子の移動・追跡に計算時間を要する事、表面は視覚的にしか把握できない事が挙げられる。

図２ SMAC 法による計算結果果１－３ VOF(Volume Of Fluid)法[3]
基本的な考え方は、セル体積に占める流体体積の比率＝流体占有率（VOF）である事を利用し、流体領域（従って自由表面）を表現する手法である[4]。
SMAC 法と同様に任意の表面形状が扱えるとともに、表面が合体／分離する問題にも適用可能である事や表面形状の変化が追跡でき、表面張力の評価で重要となる自由表面の位置・勾配・曲率などが数値的に求まるなど多くの利点を持つ。

図３ VOF 法による流体セル判定

自由表面の追跡アルゴリズムが単純で計算が速くかつメモリー消費が少ない事や、３次元問題への拡張が容易である事も津波の解析に向いている。
図４では、各セルの VOF の分率を判定して自由表面形状が求まる事が分かる。 図５では津波解析での VOF をボリュームレンダリングして海面形状が補足できる事を示した。

図４ 津波解析での VOF 法による自由表面補足
なお、以上の手法はすべて固定した格子・座標から流体の変形を扱う Euler 的手法である。
流体と共に格子・座標が移動する ALE（Arbitrary Lagrangian Eulerian）法[9]などの Lagrange的手法は、波の巻き込みなどの大変形に適用できないので津波の解析には全く不向きである。

１－４ 密度関数法，Level Set 法[4,5,6]
津波現象のような自由表面を伴う流れ解析の高精度化を実現する上で、自由表面の移動に伴う表面形状の捕捉精度を上げる事が重要である。 これは、移動に伴う形状変化を正確に記述でき、かつシャープな界面を維持できるように数値拡散が抑制されていることを含んでいる。 また、津波の解析では余り関係ないが、自由表面に働く界面張力を精度良く見積もることが重要視される場合がある。
Euler法による固定格子でVOF法で自由表面を記述する場合、数値拡散による”界面のぼやけ”が生じ気液共存域が過大に現れる可能性がある。 また、この”界面のぼやけ”は一旦発生するとなかなか解消しないので、津波解析などでは浸水域を過大に評価して精度を落とす結果につながるので注意が必要である。
VOF法の利用に伴う”界面のぼやけ”を抑制する方法として密度関数法やLevel Set法が考案されている。 密度関数法はVOF関数や実際の密度の変わりに密度関数Φを解く方法であり、密度関数Φ（0<φ<1）の判別は気相：Φ=0，界面：Φ=0.5，液相：Φ=1の原則でセル毎に行われる。 この時、各時間ステップ毎の移流拡散計算で過剰に移流した密度関数（ひいてはVOF値）を各セルの界面からの距離などに基づいて補正する事が可能であり、”界面のぼやけ”を抑制する事が出来る。
Level Set関数値（界面からの距離を表す関数）を用いた場合には、界面曲率κの算出は、界面方向の単位法線方向ベクトルn (nx ,ny , nx) を用いて次式で計算され、正確に求める事が出来るとされる。

f ：VOF 関数値、φ ：レベルセット関数値
VOF 関数値より Level Set 関数値を用いて界面曲率を算出する方がより正確であるとされるが形式的な差はない。

液体Level Set 関数 = -1
Level Set 関数 = 0
Level Set 関数 = +1

気体

図５ 界面における VOF 値と自由表面の向きと Level Set 関数

２ 津波荷重モデル
建築で使用されている津波荷重モデルは、津波波圧（荷重）が浸水深（地上での波の高さ）の 3 倍に比例するというもので、朝倉ら[13]によって検証された。

q( z )   g (3  z )2q(z)：構造設計用の進行方向の津波波圧(kN/m )3

ρ：水の単位体積質量(t/m )津波波圧

2g：重力加速度(m/s )建物



設計用浸水深3ρg
図６ 朝倉による津波荷重モデル

図4.6.3-1 (1)式による津波波圧
この式に従った津波荷重値が、本当に数値シミュレーションから求められるかが最初の検証のポイントであった。 また、図 8 に示す津波特有の荷重の時間推移が数値シミュレーションで再現可能であるかも注目された。 津波が最初に建物に衝突した際の衝撃段波波圧、及びその後、一旦波圧が減衰した後の押し波による最大重複波圧が特徴的であり、これは自由表面として波が建物に当たり一部は乗り越え一部は側方に回り込むがために起こる複雑な現象である。 荷重のオーダーとしては、津波の速度を 10m/sec とした場合 10t/㎡程度になる。

図７ 有川らによる津波荷重曲線[14]３ 避難ビル津波荷重シミュレーション
次に定量的に津波荷重の評価が可能か検証するために、3 階建て避難ビル(軒高 11m)について、海岸際にビルが立地している仮定での解析を行った（図 8, 9）。

（（

図８ 避難ビル立面図・平面図

図９ 避難ビル解析モデル

シミュレーション結果の一部を図 10, 11 に示す。 窓開口が無い場合には、建物前面で流れが減衰し低層部で大きな波圧が発生する事が確認できる。 この波圧の鉛直分布は朝倉の式で求められる浸水深 3 倍の値に一致する（図 12）。 窓開口がある場合には、建物前面の波圧は小さくなるものの内壁で水を受け止めてているために内壁に大きな荷重が掛かる。 また水平荷重の時間推移では（図 11）、衝撃段波波圧や最大重複波圧が図 7 のとおり再現出来ていることが分かる。

窓開口無し

窓開口有り

衝撃段波波圧
最大重複波圧

窓開口無し

図 10

窓開口有り

避難ビル表面流速分布・圧力分布

図 11

水平荷重時間推移

図 12 は波高 8m の場合の波圧の断面分布である。 建物
前面の波圧は窓無しの場合では、衝突直後に低層部に荷重が集中、浸水深に徐々に比例窓有りでは建窓無し，5 秒

物の中に水が投げれこみ、その傾向が低層部で顕著である
ため 1 階部分が 2 階部分より圧力が高まり 2 階床スラブを下から突き上げるのでこの部窓無し，10 秒位が破壊する可能性がある。

図 12

津波荷重力の時間推移 −異なる高さの部位の応答

図 12 には、浸水深 2m，波の速度 4.7m/sec（Fr 数 1.5）の場合の異なる高さで求めた建物前面水平荷重の時間推移について示した。 衝撃段波波圧発生時（左上グラフ）には、建物低層部に集中的に高い荷重が働いている事が分かる。
津波の建物衝突からやや時間の経過した最大重複波圧発生時（右上グラフ及び右下グラフ）には、建物前面には鉛直方向に線型の分布を持つ荷重が働くことが見てとれる。 この傾向は建物の浸水深 3 倍に比例するライン（グラフ中、赤点線）に沿っており、これは得られた荷重が朝倉の式に一致する事を示している。
また衝撃段波波圧発生直後の荷重減衰期には、荷重の鉛直分布は朝倉の式とは異なって浸水深の 3 倍以下の荷重となりそれも高さに比例しない。 これは、建物に津波が当たった後、建物を乗り越える波や側方に回りこむ波が発生し極めて複雑な流動が建物周囲に起こっている事に起因すると考えられる。
建築物全面の窓開口の有無によって津波荷重や最大波高を軽減させる効果があることは確認したが、開口部の大きさによって建築物内に流れ込む津波の流れが変化することが予想される。 そのためどのように津波荷重が低減されるのか調べる必要があり、建築物前面の開口部の面積を変えた場合の津波のシミュレーションを実施し、建築物全体に作用する津波力の変化を求め開口部の影響を調べた。
図 13 上は建築物の前面の開口面積比と最大津波力との関係で示した。 点線は開口部がない場合（開口面積比０％）を最大として線形的な低減を示したが、シミュレーション結果もほぼ同じ点線上にあることが分かった。 つまり開口面積比が 40％くらいまでであれば、建築物前面の開口面積比に応じて津波力が線形的に低減されるものと考えられる。 しかしながら、

建築物前面の開口面積比が大きくなっても、内壁や背面壁がある場合にはその部材に津波力が作用するので、必ずしも 0 に漸近しないと推定される 。 ピロティを設けた場合は水流がそこを優先的に抜けるために、津波力は大きく減衰する。
図 13 下は建築物表面の津波圧力を積分し求めた津波力の時間変化である。3 秒前後に衝撃的な津波圧力が発生しているが建築物全体ではそれほど大きくはないが、その後 5 秒前後で津波力は最大となっている。また、建築物前面の開口面積比によって、津波力が低減している。

最大津波力︵q(z)   g (3  z)

4.00E+07

Ｃｉ
2.2

3.00E+07

1.5
2.00E+07

1.00E+07

0.6N)

ピロティ構造0.0

0.00E+00
0.0

20.0

40.0

60.0

80.0

100.0

開口面積比μ(%)

津波力︵

4.00E+070.0
26.2
3.00E+07

35.0
43.5N)

51.6

2.00E+07

59.762.91.00E+070.00E+00
05101520時間(s)図 13 水平荷重に対する開口面積比の影響

３ 避難ビル津波荷重シミュレーション[1] C.W.Hirt， B.D.Nichols and N.C.Romero，LA-5852，1975 年[2] A.A.Amsden and F.H.Harlow， LA-4376, 1970 年[3] B.D.Nichols， C.W.Hirt and R.S.Hotchkiss， LA-8355，1980 年[4] C.W.Hirt，“Numerical Fluid Dynamics：A Short Course，”Lecture Notes in Tokyo，Flow Science，Inc.（1982）[5] 高谷 材料とプロセス，22 巻 1 号, ROMBUNNO.TO16, ISSN1882-8922[6] 姫野ら 日本機械学会論文集 B 編，76 巻 765 号．P778，2010 年[7] 朝位ら 土木学会論文集，810 号，P127,2006 年