非線形長波モデルと流体粒子法による津波シミュレータの開発 I_ m ρ v p h g a b a 2h b r ab a b Fang W r ab h 5 Wendland 1995 q= r ab /h a d W r ab h

(1)

非線形長波モデルと流体粒子法による津波シミュレータの開発

Development of a Tsunami Simulator Integrating the Smoothed-Particle Hydrodynamics

Method and the Nonlinear Shallow Water Wave Model

諏訪多聞

１

・今村文彦

２

・菅原大助

３

Tamon SUWA, Fumihiko IMAMURA, and Daisuke SUGAWARA

We develop a tsunami simulator integrating a 3-D fluid simulation technology that runs on large-scale parallel computers using smoothed-particle hydrodynamics (SPH) method, together with a 2-D tsunami propagation simulation technique using a nonlinear shallow water wave model. We use the 2-D simulation to calculate tsunami propagation of scale of about 1000km from epicenter to near shore. The 3-D SPH method can be used to calculate the force that a tsunami can exert on a building, and to simulate flooding patterns in urban area of at most km scale. By applying the processing power of computers to the technologies resulting from this research, we seek to contribute to improved disaster preparedness and disaster mitigation through a better understanding of tsunami's mechanisms.

1. はじめに 2011年3月11日に発生した東日本大震災による津波が東北地方を中心に甚大な被害をもたらした経験を踏まえ，将来の巨大地震・津波による被害を最小限にとどめる対策を効果的に実施するために，大規模シミュレーション技術を被災メカニズムの解明や高精度な被害予測に活用することが，強く求められている．そこで，我々は波源域から沿岸までをシームレスに扱うための津波シミュレータの開発を行ったので，その成果を報告する．粒子ベースの津波解析において，境界からの一定流量を与えることで津波を表現した場合，引き波の影響など，取り入れることが難しい要素がある．そのため，現実の津波による建造物の被災メカニズム解明や，越流・浸水過程の解析を行うためには，波源モデルなどの津波に関する想定を考慮する必要があると我々は考えた．そこで，波源から沿岸までを非線形長波モデル（TUNAMI-N2）によって解析し，その結果得られる波高・流速を初期条件・境界条件として沿岸近くでの三次元流体計算を行う手法を開発した．三次元流体計算のためには，Smoothed Particle Hydrodynamics（SPH）法を用いた流体粒子法を採用した．流体粒子法による計算は，水塊分裂や飛沫が発生する状況に対して適用性が高く，砕波や越流などの三次元挙動の再現に適すると考えられる．平面二次元計算と流体粒子法を連結して解析する手法としては，例えば五十里ら（2006）によるBoussinesqモデルとMoving Particle Semi-implicit（MPS）法とを連結する手法がある．これは三次元計算領域の境界に仮想造波１理博富士通（株）２正会員工博東北大学教授災害科学国際研究所３工博東北大学助教災害科学国際研究所板を設定し，Boussinesqモデルによって得られる波高と流速を再現するようにその仮想造波板を制御するというものであり，任意の向きの入射波を扱えるという特徴がある．我々は，津波のように計算領域外から大量の流体が継続的に流入し続ける状況を扱う状況では，境界からの流入を直接扱える接続手法が適しているものと考え，以下に述べる手法を提案する． 2. 手法 波源から沿岸域までの津波伝播を解析するために，我々は以下に述べる非線形長波モデルを採用した． ………（1） ………（2） ………（3） ここで，hは水位，MとNはそれぞれxとy軸方向の線流量， Dは全水深，nはManningの粗度係数である．我々は有限 差分法を用いて上記の式を数値的に解くソルバーの一つ TUNAMI-N2を採用した．TUNAMI-N2では水位と線流量についてスタガード格子上で値を評価した中心差分法を使い，連続方程式を差分化する．運動方程式に対しても同様に差分化を行うが，移流項については，風上差分式を，摩擦項については陰解法を採用する（Gotoら， 1997）．我々はSPH法の運動方程式と連続方程式として，以下の形式を用いた（Suwaら，2013）.

(2)

…………（4） ………（5） ここで，m，ρ，v，p，hはそれぞれ粒子の質量，密度， 速度，圧力，定数ベクトルg→は重力加速度である．添え 字a，bは粒子のインデックスで，和は粒子aから半径2h の範囲にある粒子bにわたって取る．また，rabは粒子a， b間の相対位置ベクトルである．式（4）が質量の連続の式，式（5）が運動方程式である．本研究において式（5）を用いた理由は，これがエネルギー保存則から導かれたものであり（Fangら，2009），数値的エネルギー減衰の防止に適していると考えたからである．なお，越波や波高の予測に関する中川ら（2011）の結果によれば，粘性の影響が些少であったため，式（5）では粘性項を省略 した．関数W（rab，h）は和の重みを表すカーネル関数であり，次の5次スプラインを採用した（Wendland， 1995）． …………（6） ここで，q=|rab|/hと表記した．また，adはW（rab，h）を空間全体で積分した値が1となるように決まる規格化係 数であり，3次元では21/（16π h3_{）である．粒子の圧力p} は密度ρから状態方程式 ………（7）を用いて算出する微圧縮近似を用いた．ここで，ρ0は基準密度で水の物性値である1000kgm-3_{を用い，p} 0は基準圧力で，2.25×107_{Paとした．p} 0は数値計算上のパラメータであり，典型的な運動エネルギーを持ったSPH粒子の密度が基準密度から大きく外れない程度に大きな値を選んだ．非線形長波モデルとSPH法とを連結させるため，非線形長波モデルの2次元計算結果をSPH法の流体粒子計算の初期条件と境界条件として与えることとした．SPH法を用いて解析する領域内に，初期条件として非線形長波モデルから得られた流速を初速度として持ち，粒子配置として水底から2次元計算で得られた水位までに等間隔で 図-2 三次元SPH計算による津波の浸水の様子．縦軸，横軸は解析領域のx，y空間座標．上部の網掛けは陸域，点がSPH粒子，下部の斜 めの網掛けが二次元非線形長波モデルによる浸水域を示す． 図-1 東松島地区における検証計算で用いた地図データ．（a）-（d）は二次元非線形長波モデルの解析で用いた水深・標高データ．左から順に，1350m，450m，150m，50m の格子データ．（e）は三次元SPH 計算を実行した領域（地図（d）の中の長方形で示した領域）の標高データ

(3)

置くこととした． 境界近く（初期粒子間隔dxの範囲）にある粒子につい ては2次元計算で得られた流速をSPH法計算の全時間において与えるものとした．今回用いた2次元計算では流体速度は鉛直方向に一様，かつ，水平速度成分のみを持つものと仮定したため，境界条件としてSPH粒子に与える速度は以下のとおりとする．速度はSPH粒子の鉛直位置座標によらず，水平位置座標のみによるものとし，また，与えられる速度は水平速度成分のみを持ち，鉛直速度成分は0とする．密度については基準密度ρ0を与えるものとする．境界における，流入・流出を表現するために境界面を SPH粒子の断面積程度の面要素に区切り，その各面要素の位置において，SPH粒子の生成・消滅を行うモデルを導入した．流出のためには単に境界位置を越えた粒子を消滅させればよい．流入のためには，この面要素の位置 における流入・流出体積Vを毎タイムステップ更新し， 必要に応じてSPH粒子を生成することとした．具体的な手順としては，各面要素の位置において，2次元伝播計 算結果に基づいて設定した速度ベクトルvと，面要素の 面積S，法線ベクトルnを用い，以下の式に基づいてdt後 の流入・流出体積Vを更新する． ………（8） この時，面要素に設定される速度ベクトルvはその水平 位置座標にのみ依存し，鉛直座標には寄らないものとする．ただし，2次元計算の結果得られた水位よりも高い 位置に存在する面要素に対してはv=0とする．更新した 結果，流入・流出体積がSPH粒子一つに対応する体積 （基準体積V0と呼ぶ）よりも大きくなったら，SPH粒子を 一つ生成し，流入・流出体積からV0を減ずる．一方，更 新した結果，流入・流出体積が-V0を下回ったなら，V0を 加える．こうして，流入・流出体積は常に-V0から+V0の範囲に留まる．これは流入については粒子生成によって表現している一方，流出については，境界を越えた粒子を単に消滅させているため，流入・流出体積が負の大きな値を取ると，粒子を消滅させていることとの整合性が取れなくなるためである．新たに生成されたSPH粒子の位置座標は境界の面要素と同じ位置に置き，速度は2次 元伝播計算結果に基づき速度ベクトルvを，密度は基準 密度ρ0を与えるものとする． 3. 結果 （1）東松島地区 本手法の妥当性を検証するため，東松島地区における検証計算を行った．非線形長波モデルの解析の計算条件として，以下を設定した：地図データとして1350m， 450m，150m，50mの格子点上の水深・標高データを用い た．図-1に示す通り，最大の地図は東日本の広域を覆い， 最小の地図では東松島地区に注目している．初期水位形状は東日本大震災の津波初期波形に相当するように選ん だ．図-1（d）に示す0.5km×1.7kmの長方形の領域を三次 元SPH法の対象領域として設定した．この領域を図-1（e） 図-4 広域計算結果と連結した境界流入条件で三次元SPH計算を行った結果（左）と一様流入を仮定した結果（右）の比較 図-3 三次元SPH計算で得られた波形の推移．図-1（e）に示した線で切った断面の図．点がSPH粒子，緑線（下）が地図データによる標 高・水深，赤線（上）が非線形長波モデルによって得られた水位．

(4)

に示す．三次元SPH計算の初期時刻として，非線形長波モデルによる計算で沿岸に津波が到達する時刻（3,360秒時点）を設定した．流体粒子を水底から1m間隔で非線形長波モデルの結果得られた水位に至るまで格子点上に配置して，初期粒子配置とし，流体粒子の初速度は非線形長波モデルの結果から換算された水平方向速度を与え，鉛直成分は持たないものとして設定した．非線形長波モデル，及び，三次元SPH法による東松島 地区の津波シミュレーション結果を図-2から図-4に示す． 図-2では，上空から見下ろす向きで非線形長波モデルと SPH法のそれぞれで計算した浸水域をプロットした．三次元SPH法の計算を開始した時刻を0secとして20秒おきに表示している．左の粒子分布をプロットしたものが SPH計算の結果であり，右の網掛けで浸水域を示したものが非線形長波モデルによる結果である．これを見ると， SPH法による結果は非線形長波モデルによるものと良い 一致を示している．また，図-3では，三次元SPH計算を 行った領域の鉛直断面図（図-1（e）に示した直線に沿っ て切ったもの）を示す．SPH法による結果を粒子のプロットとして示し，非線形長波モデルによる結果を水位を表す実線として示している．ここでも，SPH法による結果は非線形長波モデルによるものと良い一致を示しており，今回採用した連結手法が問題なく動作することを確認できた． 図-4では，三次元SPH計算として，今回採用した連結 モデルによる結果と，境界からの一様流入を仮定したモ デルの結果を比較した．計算領域は図-1（e）に示した東 松島地区であり，水位（流体粒子の鉛直座標）に基づいてカラーマップを示した．一様流入を仮定していたのでは現れない，側面から回り込んで流入してくる成分が表現されているのが分かる．今回実行した三次元SPH計算においては初期粒子数としておよそ230万体を用い，流入によって増加し40秒時点ではおよそ336万体の流体粒子を扱った．計算時間としては，ワークステーション（Fujitsu PRIMERGY RX200S3） 8ノード（64CPUコア）を用いて20時間程度を要した． （2）一次元伝播計算と二次元SPH法による検証計算 今回開発した連結手法の有効性を検証するため，我々 は以下の検証計算を行った．図-5に本計算に用いた水槽 模型を示す．この模型では，左端に1/10勾配のスロープを置き，右端には鉛直の固定壁を設定した．水槽の深さ と長さはそれぞれ，15mと500mである．図-5は表示の都 合上，縦横の縮尺が大きく異なることを注意する．この 2次元水槽における表面波の伝播を1次元非線形長波モデル，2次元SPH法，二手法の連結のそれぞれで解析し結果を比較した．初期波形として， ………（9）をあたえ，これが孤立波として，スロープに打ちあがっ ていく様子を解析する．ここで，x0=500m，λ=75mであ り，初期波高Aは3.0mと15.0mの2通りを行った． 図-6にx=100m（初期水深0m）の地点における水位履歴 を示す．A=15.0mを右に，A=3.0mの結果を左に示す． SPH法による結果は青の実線で示し，非線形長波モデル による結果は赤の破線で示した．A=15mの場合，SPH法 にのみピークにおいて鋭く立ち上がる構造が見られる一 方，A=3mの場合ではSPH法と非線形長波モデルの結果と は同様の結果を示している．水深と波高の比率が1/10程度の範囲では，長波近似の結果が直接運動方程式を解いたSPH法の結果とよく一致すると考えられる．次に我々は流入・流出が起こる個所で接続することの 妥当性を検証するため，A=3mで，流入・流出境界を x=150mに設定したSPH法計算を実行した．つまり， x=150mより左の領域でのみSPH法の計算を行い，領域の右側からは一次元伝播計算の結果に応じて流入・流出が 図-5 一次元伝播計算と二次元SPH法による検証計算で用いた 水槽模型． 図-6 一次元伝播計算と二次元SPH法による検証計算の結果得られた水位履歴．初期波高A=15m（右）とA=3.0m（左）．縦軸と横軸はそ れぞれ時刻（秒）と水位（m）．

(5)

起こることとした．図-6左のグラフで，連結SPH計算に おけるx=100mでの水位履歴を黒の一点鎖線で示した． 結果は全体をSPH法のみで解析した結果，非線形長波モデルのみで解析した結果のいずれともよい一致を示した．これによって，計算初期において，波の全体が計算領域に含まれていない状況から計算を開始しても，流入・流出を境界において適切に設定することで，波の伝播を扱えることが確認できた． 4. 結論 東松島地区において本手法の妥当性を確認したところ，陸域に到達するまでの伝播が非線形長波モデルと3 次元SPH法計算でよく一致することが示された．このことから，境界の流入・流出量を合わせることで計算領域内においても，二次元長波モデルと整合する流れが実現していると考えられる．また，東松島の伝播における断面波形からは，流入量のみならず，波形についても，想定を反映した結果を与えている可能性が示唆される．ま た，図-4に示されるように，計算領域の一方からの一様 流入を仮定していたのでは現れない，側面から回り込んで流入してくる成分を扱うことができていることも確認できた．今回の計算では，実際の地図情報と広域計算結果に基づいて，8コアのワークステーションを160ノード時間積だけ用いることで，0.5km×1.7kmの領域を1m径の粒子を用いて，3次元SPH法の計算ができることを示した．現在，スーパーコンピュータとして数万ノード規模のシステムが稼働しており，仮に160万ノード時間積という計算リソース（4万ノード程度を2日程度）を用いることができれば，以下のような解析が可能になる．同程度の領域・時間に対しては，0.1m径の高解像度での解析．これは遡上するサージフロントをより正確に扱える可能性がある．0.5m径程度の解像度を用いるならば，10km四方の領域を100秒程度解析することが可能となる．これは湾や港を一つ分扱うのに十分な大きさと言える．今後は今回開発された技術を用いて，実際の被災地のデータを用いたさらなる検証を行い，シミュレータの信頼性を高めることを予定している．さらに将来的には，本シミュレータを活用することで，津波被害予測を高精度化し，災害に強い街づくりに貢献することを目指したいと考えている．参　考　文　献五十里洋行・後藤仁志・酒井哲郎・奥田一弘（2006）：三次元数値波動水槽のための粒子法とBoussinesqモデルとのハイブリッド化，海岸工学論文集，第53巻，pp. 11-15. 中川知和・片岡保人・竹鼻直人・諏訪多聞・風間正喜（2011）：粒子法による護岸越波現象の数値計算，土木学会論文集 B3 （海洋開発），Vol. 67, pp. 268-273.

Fang, J., A. Parriaux, M. Rentschler, and C. Ancey (2009) : Improved SPH methods for simulating free surface flows of viscous fluids, Applied numerical mathematics, Vol. 59, pp. 251-271. Goto, C., Y. Ogawa, N. Shuto, and F. Imamura (1997) : Numerical

method of tsunami simulation with the leap-frog scheme (IUGG/IOC Time Project), IOC Manual, UNESCO, No. 35, 126 p. Suwa, T., T. Nakagawa, and K. Murakami (2013) : A Study of the

Wave Transformation Passing over an Artificial Reef using SPH Method, Journal of computational science and technology, Vol. 7, No. 2, pp. 126-133.

Wendland, H. (1995): Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational mathematics, Vol. 4, pp. 389-396.