技術報告 物理探査 第 64 巻第 3 号 (2011) 187-195 頁 BUTSURI-TANSA, Vol. 64, No. 3 (2011) pp. 187-195 地震波動伝播シミュレーションにおける粒子法の適用性に関する研究 武川順一*・三ケ田均*・後藤忠徳* 要 旨 地震国である我が国では，地震発生時の構造物の安全性や液状化対策などを考えることは，安全・ 安心な社会の実現にとって重要な課題である。地震防災を考える時，地震発生により生じる強震動を 正確に予測することは，被害予測や社会基盤施設の復旧シミュレーションなどをおこなう上で最も基 本的な情報となる。 従来から，地震波動伝播シミュレーションにはスタガード格子を用いた差分法が広く用いられてき た。また，自由境界条件の取り扱いの容易さなどから有限要素法の適用も試みられている。本研究で は，粒子法の一種である Hamiltonian Particle Method (HPM) に着目し，地震波動伝播シミュレーショ ンへの適用性を検討する。HPM は，有限要素法と同様に自由境界条件の取り扱いが容易であり，任意 の地表面形状を有する解析対象を扱いやすい。また，離散化には粒子点座標を必要とするのみで，有 限要素法のように節点と要素のコネクティビティを必要としないといった利点も有する。一方で， HPM はシンプレクティックスキームを適用することでエネルギー保存精度に優れた解析が可能であ るが，粒子が局所的な振動を生じてしまい大域的な振動が減衰してしまうという問題がある。この問 題を解決するために，振動抑制のための人工的な力の導入がおこなわれており，梁の自由振動解析な どでその効果が検討されている。本研究では，地盤を伝播する地震波動のシミュレーションをする際 に，この人工的な力がどのように結果に影響を及ぼすかに注目し，詳細な検討をおこなった。また， 地震による強震動を考慮する際には，地表面を伝播する表面波を正確に再現できることは重要である ので，表面波の再現性を他の数値解析手法による結果と比較することで詳細に検討した。その結果， HPM を地震波動場のモデリングに適用する場合には，人工的な力の導入は不可欠であり，その大きさ を規定する係数には媒質のヤング率の 0.5 倍程度の値を使用することが望ましいことがわかった。ま た，水平でない地表面を有するモデルにおける表面波伝播についてもシミュレーションをおこなった ところ，表面波伝播に関しては実体波部分よりも数値的な振動が生じやすく，振源周波数や粒子間隔 を適切に設定する必要があることがわかった。 キーワード：粒子法・数値シミュレーション・弾性波・表面波 1. はじめに （Graves, 1996）。差分法は理論も簡便であり，要求され る計算機性能も比較的低くて済むという利点を持つ。し 地震国である我が国では，安心・安全な社会の構築の かし，凹凸を有する地表面の境界条件を導入することが ために，地震発生時の構造物や社会基盤施設の耐震性を 難しく，任意の地表面形状を解析に取り込むことは容易 考えることは重要な課題である。構造物などの耐震性を でない。一方で，自由境界条件の取り扱いの容易さなど 考慮する際の入力となる強震動予測には，従来から，ス から，有限要素法による地震動シミュレーションも試み タガ ード格子を用 いた差分法が 広く用いられ てきた られている。通常の有限要素法では，節点と要素のコネ 2011 年 5 月 9 日原稿受付；2011年 8 月 10 日受理 * 京都大学大学院工学研究科 〒615-8540 京都府京都市西京区京都大学桂 © 2011 SEGJ 187 188 物理探査 第 64 巻第 3 号 (2011) クティビティの生成に労力を要することや，剛性マトリ 議論する。また，振源の周波数を変化させた計算もおこ クス成分を記憶するために膨大なメモリーを必要とする ない，周波数の変化による影響を検討する。最後に，平 などの問題があるため，ボクセル要素の適用が提案され 坦でない自由表面を有するモデルにおける表面波伝播の ている（Koketsu et al., 2004） 。ボクセル要素を用いるこ シミュレーションをおこない，有限要素法による結果と とでこれらの問題をある程度解決でき，さらに任意のリ 比較することでその適用性を検討する。 ファインメントが可能となるが，重複要素の設定などプ リプロセスにおいてデータ構造は多少複雑となってしま 2. 計算手法 う。また，Cundall and Strack (1979) により提案された， 質点の力学に基づく Discrete Element Method (DEM) も 地震波伝播シミュレーションに適用されており（Toomey 2.1 変形勾配テンソル 変形勾配テンソル F とは，Fig.1 に示すように，物質内 and Bean, 2000），流体流動と固体の連成解析にも応用さ のある微小なベクトル dX が変形を受けた後に dx となっ れている（O’Brien and Bean, 2004a）。しかし，DEM では たとき，dx=F∙dX という線形変換作用素として定義され バネ定数などの微視的パラメータを設定するためにキャ るもので，ある物質点近傍の相対的な変形を表すもので リブレーションを必要とするなどの問題点がある。 ある（久田・野口，1995） 。HPM では解析対象を粒子に 一方で，粒子法と呼ばれる数値解析手法がこれまでに よって離散化し，重み付き最小二乗法により粒子点上で 数多く提案され，様々な物理シミュレーションに適用さ の 変 形 勾 配 テ ン ソ ル を 計 算 し て い る （ Suzuki and れてきている。Smoothed particle hydrodynamics (SPH) Koshizuka, 2008）。粒子 i における変形勾配テンソルを Fi （ Lucy, 1977 ） や Moving particle semi-implicit (MPS) とし，誤差関数 ei を以下のように定義する。 （Koshizuka and Oka, 1996）などが提案され，流体解析や 固体解析，もしくはそれらの連成解析に適用されてきた ଶ 用いて粒子位置で応力・ひずみを計算する Hamiltonian e୧ = ∑୨ห۴୧ ‫ܚ‬୧୨଴ − ‫ܚ‬୧୨ ห w୧୨ particle method (HPM) が提案され，シンプレクティック この誤差関数 ei を最小にするように粒子 i における変 スキームを適用することでエネルギー保存に優れた手法 形勾配テンソル Fi を決定する。ここで，‫ܚ‬୧୨଴ と‫ܚ‬୧୨は粒子 ij （例えば，Chikazawa et al., 2001） 。また，最小二乗法を が弾性体の解析に適用されている（Suzuki and Koshizuka, (1) 間の初期相対位置ベクトルと現在の相対位置ベクトルで 2008）。HPM は粒子の局所的な振動によって大域的な振 ある。wij は粒子 ij 間の重みであり，適当な影響半径内で 動が減衰してしまうという問題点があるため，これを人 のみ値をもつ関数として定義される。式(1)中の総和は， 工的な力を導入することにより抑制する工夫がなされて 粒子 i の影響半径内にある全粒子 j に対して定義される。 いる（Kondo et al., 2010） 。HPM は，有限要素法と同様に 式(1)を最小にする変形勾配テンソルは以下のようにな 自由境界条件の導入が容易であるため，任意の地表面形 る。 状を有する解析対象を容易にモデル化できる（Suzuki et al., 2007）。また，離散化には粒子の位置情報のみを必要 とするため，有限要素法のように節点と要素のコネクテ ィビティを生成する作業をプリプロセスの段階で要求す ることはなく，入力データが簡便になるという利点も有 する。しかし，局所的な振動を抑制するための人工的な ۴୧ = ∑୨ ‫ܚ‬୧୨ ⨂‫ܚ‬୧୨଴ w୧୨ ‫ିۯ‬ଵ ୧ ‫ۯ‬୧ = ∑୨ ‫ܚ‬୧୨଴ ⨂‫ܚ‬୧୨଴ w୧୨ (3) 以上のように，変形勾配テンソルは粒子座標のみによ ってあらわされる。 力の決定には物理的な根拠はなく，解析対象に応じて適 切な値を選ぶ必要がある。そこで，本研究では，地震波 動伝播シミュレーションに HPM を適用する際に，人工 的な力がどのように計算結果に影響を及ぼすかを検討し， HPM の地震波動伝播シミュレーションへの適用性を考 察する。 最初に，等方均質媒質内を伝播する地震波のシミュレ ーションを，人工的な力を考慮した場合と考慮しない場 合でおこない，人工的な力がどのように地震波伝播の再 現性に影響しているかを調べる。続いて，人工的な力の 大きさを変化させた計算結果をスタガード格子差分法に よる結果と比較することで，適切な力の大きさについて (2) Fig.1. Conceptual diagram of deformation gradient tensor. 武川・三ケ田・後藤：地震波動伝播シミュレーションにおける粒子法の適用性に関する研究 189 のヤング率 E と同程度の適切な値を選ぶ必要がある（近 2.2 粒子の運動方程式 式(2), (3)で計算される変形勾配テンソルを用いて，以 藤ほか，2007） 。地震波動伝播シミュレーションをおこな ෡ の値に関しては次章で詳しく検討する。 う上での適切なE 下のように Green-Lagrange ひずみテンソル E および第 人工的な力の大きさは，式(8)からもわかるように，式 2Piola-Kirchhoff 応力テンソル S が計算される。 (2),(3)における変形勾配テンソルの近似が適切である場 合にはゼロとなり，計算結果に影響をおよぼさない。 ۳ = (۴ ୘ ۴ − ۷)/ 2 ‫ = ܁‬2μ۳ + λtr(۳)۷ (4) 3. 数値計算 (5) 式(4), (5)を用いて，変形により弾性体に蓄えられるポ テンシャルエネルギーV は， 本章では，前章で説明した HPM の理論を用いて，実 体波・表面波伝播シミュレーションをおこない，その結 果を差分法や有限要素法による結果と比較することで， ଵ V = ∑୧ ۳୧ : ‫܁‬୧ ΔA୧ ଶ (6) と表現できる。ここで，ΔAi は粒子 i が占める体積をあら わす。式(6)を用いることで，粒子 i の運動方程式が以下 のように導かれる。 HPM の地震波動伝播シミュレーションに対する適用性 の検討をおこなう。 3.1 均質媒質中を伝播する実体波 2.で述べた手法を用いて，Fig.2 に示す均質媒質中を伝 播する地震波のシミュレーションをおこない，スタガー ド格 子を用いた差 分法による結 果と比較する ことで ିଵ ଴ ଴ ρ ∂‫ܞ‬୧ / ∂t = − ∂V/ ∂‫ܠ‬୧ = ∑୨൫۴୧ ‫܁‬୧ ‫ିۯ‬ଵ ୧ ‫ܚ‬୧୨ + ۴୨ ‫܁‬୨ ‫ ۯ‬୨ ‫ܚ‬୧୨ ൯w୧୨ HPM による地震波動場の再現性を検討する。最初に，人 上式にシンプレクティックスキームを適用すること その後，人工的な力を決める係数の大きさ Cart を様々に (7) 工的な力を考慮した場合と考慮しない場合の結果を比較 し， 弾性波動伝播における人工的な力の影響を検討する。 で，エネルギー保存に優れた計算が可能となる。本手法 変化させ， 適切な人工的な力の大きさについて検討する。 の詳細は，文献（Suzuki and Koshizuka, 2008 など）を参 媒質の P 波速度は 4000m/s，S 波速度は 2310m/s，密度は 照されたい。 2700kg/m3 とした。また，振源関数にはリッカーウェブ レット（中心周波数 50Hz）を使用し，振源位置において 2.3 人工的な力 x 方向変位速度場に与えた。初期配置における粒子間距 式(7)をそのまま用いたのでは，粒子が局所的な振動を 起こしてしまうことが知られている。この局所的な振動 離は 1m，時間刻みは 0.1ms とした。影響半径としては初 期粒子間距離の 1.9 倍を用いた。 に物理的な意味はなく，精度のよい計算のためにはこの ここで，影響半径について少し詳しく説明する。影響 振動は抑制されなければならない。本節では，粒子の局 半径を大きくしすぎると，式(2),(3)で計算をおこなう近 所的な振動を抑制する方法について説明する。Kondo et 接粒子数が増え，結果として計算時間の増大を招く。ま al. (2010)は，以下の人工的なポテンシャル力を導入する た，影響半径内の変形が線形に近ければ，式(2),(3)によ ことで HPM における局所的な振動を抑制した。 る変形勾配テンソルの近似精度がよくなる（すなわち， ܎୧ୟ୰୲ = Cୟ୰୲ ∑୨൛൫‫ܚ‬୧୨ − ۴୧ ‫ܚ‬୧୨଴ ൯ + ൫‫ܚ‬୧୨ − ۴୨ ‫ܚ‬୧୨଴ ൯ൟw୧୨ (8) ここで，Cart は人工的な力の大きさを決める係数であ る。この係数は， Cୟ୰୲ = ෡ ୫୉ ୢ మ ஡ ∑ ቚ‫ ܚ‬బ ቚ ୵ ౠ ౟ౠ ౟ౠ (9) となる。ここで，m は粒子の質量，ρ は密度，d は空間 ෡ は人工的な力の大きさを決定するパラ の次元である。E メータでありヤング率 E と同じ次元を持つが，この人工 的な力の大きさの決定に物理的な根拠はなく，解析対象 Fig.2. Geometry of the homogeneous model used in the numerical experiment. 物理探査 190 第 64 巻第 3 号 (2011) ‫ܚ‬୧୨ − ۴୧ ‫ܚ‬୧୨଴ ≈ 0）ことからも，影響半径はできるだけコン 所的な振動に起因すると思われる不自然な振動が確認で パクトなサポートを持つように設定されることが望まし きる。この振動に物理的な意味はなく，応力評価位置が い。 影響半径を初期粒子間距離の 1.9 倍に設定した場合， 粒子位置と一致するために生じる数値的な振動であり， 影響半径内の粒子は 2 次元の場合，上下・左右・斜め方 抑制されるべきものである。一方で，人工的な力を考慮 向の計 8 個となり，対象粒子を正方形状に取り囲む粒子 した場合（Fig.3b）では，P 波と S 波がモデル中央の振源 との相対変位から変形勾配テンソルが計算されることに から 同心円状に伝 播している様 子のみが確認 でき， なる。Kondo et al. (2010) では影響半径として初期粒子間 Fig.3a に見られた局所的な振動は見られない。これによ 距離の 2.5 倍を用いているが，この場合だと隣接粒子数 り，人工的な力の導入が数値的な振動を抑制しているこ は 20 個となり計算時間は増大してしまう。しかし，他の とが確認できる。このように，地震波動場のモデリング 粒子法（elastic lattice model など）においては，対象粒子 に HPM を適用する場合は，人工的な力による局所的な を取り囲むような周囲粒子との相互作用を計算すること 振動抑制は必須であると考えられる。 で十分な精度が得られており（Toomey and Bean, 2000; 次に，人工的な力の大きさを変化させた場合に，それ Valle-Garcia and Sanchez-Sesma, 2003; O’Brien and Bean, がどのようにシミュレーション結果に影響をおよぼすか ෡ を変化さ 検討する。人工的な力の大きさは，式(9)中のE 2004b; O’Brien et al., 2009），HPM においても上述のよう な対象粒子の選択で十分な計算精度が得られるものと推 察される。なお，以下に検討する内容に対しては，影響 半径を 2.5 倍としたケースでも同様の計算をおこなって ෡ をヤング率 E の 0.25 せることで調整できる。本節では，E 倍，0.5 倍，1.0 倍と変化させて検討をおこなった。Fig.4 ෡ を変化させた場合の，受振器位置における y 方向変 はE おり，本研究で得られた結論と同様の結果が得られてい 位速度の受振記録を示している。各図中には同じモデル ることを記しておく。 に対するスタガード格子を用いた差分法による結果を点 Fig.3 に 100ms 後における y 方向変位速度場のスナッ 線により併せて示している。図より，各ケースにおいて プショットを示す。Fig.3a が人工的な力を考慮しない場 0.08s 付近から見られる P 波部分についてはほとんど差 合の結果，Fig.3b が人工的な力を考慮した場合（ただし， ෡ =E）の結果である。図より，人工的な力を考慮しない E が見られないが，0.13s 付近から見られる S 波部分につい 場合（Fig.3a）ではモデル中央部から上下左右方向に局 ෡ = 0.25 × E (a) E (a) without artificial force, ෡ = 0.5 × E (b) E (m/s) (b) with artificial force, ෡ = 1.0 × E (c) E Fig.4. The waveforms of displacement velocity field in the y-direction at the receiver point. The source frequency is Fig.3. The snapshots of displacement velocity field in y-direction with and without artificial force after 100ms. set to 50Hz. Solid and dotted lines represent the results from HPM and FDM, respectively. 武川・三ケ田・後藤：地震波動伝播シミュレーションにおける粒子法の適用性に関する研究 ては各ケースで違いが確認できる。局所的な振動を抑制 191 が大きすぎると結果に悪影響をおよぼすことがわかる。 する人工的な力は変形勾配テンソルの近似誤差 次に，振源関数の周波数を 100Hz として同じモデルに （‫ܚ‬୧୨ − ۴୧ ‫ܚ‬୧୨଴ ）の大きさに依存して作用する。P 波部分に 対して計算をおこなった。Fig.5 にシミュレーションによ ついて各ケースでほとんど差が見られなかったのは，縦 り受振器で得られた結果を示す。先程と同様に，差分法 波による変形は式(2),(3)でうまく近似できているためで による結果も併せて示している。図より，周波数を 50Hz ෡を とした前ケースと結果に大きな違いが確認できる。E あると考えられる（すなわち，‫ܚ‬୧୨ − ۴୧ ‫ܚ‬୧୨଴ ≈ 0）。Kondo et al. (2010) は HPM で梁を伝播する弾性波について計算した 結果を検討している。それによると，P 波伝播について ヤング率 E の 0.25 倍，1.0 倍としたケースでは，S 波部 ෡ をヤング 分に非現実的な振動が確認できる。一方で，E は HPM における変形勾配テンソルの近似は良好であり， 率 E の 0.5 倍としたケースでは他のケースで見られるよ 人工的な力の有無によらず良い結果を示すが，S 波伝播 うな不自然な振動はなく，参考として示した差分法によ については人工的な力を考慮しないと不自然な振動モー る結果とも良い一致を示している。このことから，人工 ドが現れることがわかっている。本研究による結果も， 的な力は小さすぎると局所的な振動を抑えきれず，結果 P 波部分については人工的な力の大きさによる差はほと に悪影響をおよぼすことがわかった。また，P 波部分に んど見られず，S 波部分についてのみ大きな違いが見ら 関しては，各ケースとも差分法による結果と比較的良い れる結果となった。これは，Kondo et al. (2010)における ෡ の大きさによる結果の違 結果と整合するものである。E 一致を示しているが，周波数を 50Hz とした時に比べる とわずかな違いが見られる。これは，振源の周波数を増 ෡ をヤング率 E の 0.25 倍，0.5 倍とした場 いについては，E 加させることで伝播する地震波の波長が短くなるため， 合の結果にはそれほどの差は見られない。しかし，1.0 影響半径内の変形が線形な状態から遠ざかり，式(2),(3) 倍の値を用いた場合では，S 波部分に他のケースとの明 による変形勾配テンソルの近似精度が低下するためであ 確な違いが確認できる。参照として示したスタガード格 ると考えられる。この問題は粒子間隔を小さくすること 子を用いた差分法による結果との誤差から，人工的な力 により解決できると考えられ（Kondo et al., 2010） ，伝播 する地震波の波長を考慮した上で粒子間隔を適切に設定 する必要が示唆される。 このことを確認するため，媒質定数や振源周波数など の条件を先程のモデルと同様にし，粒子間隔だけを半分 に変更したモデルを作成し，シミュレーションをおこな った。Fig.6 に受振器で得られた結果を示す。各ケースに ෡ = 0.25 × E (a) E おいて Fig.5 で見られた非現実的な波形の振動が無くな ෡ の大きさによる結果の違いに っていることがわかる。E ついては，基本的に Fig.4 で見られるものと同じ傾向で ෡ をヤング率 E の 0.25 倍， あることがわかる。すわわち，E 0.5 倍としたケースでは P 波，S 波共に差分法の結果とよ ෡ をヤング率 E の 1.0 倍としたケース く一致しているが，E では S 波部分に大きな違いが見られる。Fig.4 と Fig.6 に おける波長と粒子間隔の相対的な関係は同じであるため， ෡ = 0.5 × E (b) E このような結果が得られたものと推察される。また，本 研究では媒質のポアソン比などの媒質定数も変化させて 同様のシミュレーションをおこなった結果，同様の結果 が得られたことを付記しておく。 以上の結果より，HPM による地震波動場モデリングを おこなう際は，人工的な力の導入は必須であるが，人工 ෡ の大きさにより結果に不自然 的な力の大きさを決めるE ෡ = 1.0 × E (c) E Fig.5. The waveforms of displacement velocity field in the y-direction at the receiver point. The source frequency is set to 100Hz. Solid and dotted lines represent the result from HPM and FDM, respectively. な振動が見られることが確認された。本研究における結 ෡ の大きさはヤング 果より，実体波の伝播については，E 率 E の 0.5 倍程度の値が適切であることがわかった。 3.2 表面波伝播 次に，Fig.7 に示すモデルを用いて，HPM の表面波伝 物理探査 192 第 64 巻第 3 号 (2011) 播シミュレーションに対する適用性を検討する。HPM で あることを意味する。このように，Fig.7 に示すような地 は，解析対象外部の連続体の密度が解析対象の連続体の 表面形状を有するモデルでも，HPM では容易に計算に取 密度に比べて十分に小さい場合，解析対象外部の粒子か り込むことが可能である。さらに，離散化に必要な情報 らの寄与を無視することで自由境界条件を設定すること は粒子の座標のみであり，有限要素法のように節点と要 ができる（鈴木ほか，2005）。地震動伝播シミュレーショ 素のコネクティビティを必要とせず，プリプロセスにお ンにおける地表面境界の場合，一般的に解析対象外部の けるデータは単純になるという利点も有する。 媒質は空気であるため，解析対象である地盤と比較して 本節では，比較に用いる数値解析手法として差分法で 密度は十分に小さいと仮定できる。そのため，解析対象 はなく有限要素法を採用する。要素には４節点四角形要 外部の粒子からの寄与を無視することで自由境界条件の 素を用いた。前節と同様に，人工的な力を決める係数の 導入が可能である。これは，有限要素法のように特別な 大きさ Cart を様々に変化させ，人工的な力の大きさがど 処理を必要とすることなく自由境界条件の導入が可能で のように表面波伝播に影響するかを検討する。媒質定数 は前節と同様とした。また，振源関数としては前節と同 様に 50Hz のリッカーウェブレットを用い，x 方向変位速 度場に与えた。 ෡ をヤング率 E の 0.5 倍とした時の各時刻に Fig.8 に，E おけるスナップショットを示す。図より，HPM により複 雑な地表面形状を有するモデルに対しても，その境界を 走る実体波・表面波の伝播が再現できていることが確認 ෡ = 0.25 × E (a) E ෡ = 0.5 × E (b) E ෡ = 1.0 × E (c) E Fig.6. The waveforms of displacement velocity field in the できる。また，Fig.9 に受振器位置における粒子の運動を (a) after 100ms, (b) after 150ms, y-direction at the receiver point. The particle spacing is 0.5m. (m/s) (c) after 200ms, Fig.8. The snapshots of the displacement velocity field in the Fig.7. Geometry of the numerical model with an arbitrary free surface. y-direction. After (a)7.5s, (b) 10s, (c) 15s, (d) 20s, respectively. 武川・三ケ田・後藤：地震波動伝播シミュレーションにおける粒子法の適用性に関する研究 193 示す。実体波の部分と思われる直線的な運動と，レイリ 速度の記録を示す。参考として，同様のモデルに対して ー波の特徴である楕円運動が確認できる。有限要素法と 有限要素法により計算した波形を点線で並べて示してい 同様に，特別な処理を必要とすることなく自由境界条件 る。図より，各ケースにおいて P 波，S 波，表面波の到 ෡ をヤング率 E の 1.0 倍と 達が順に確認できる。しかし，E を取り入れることができるのが HPM を用いる利点の一 つである。 ෡ をヤング率 E の 0.25 倍，0.5 倍， 続いて Fig.10 に，E 1.0 倍と変化させた時の受振器位置における y 方向変位 したケース（Fig.10c）では，0.3s 付近に見られる表面波 部分に不自然な振動が見られる。有限要素法による結果 との誤差から，前節と同様に大きすぎる人工的な力が結 果に悪影響をおよぼしたものと考えられる。前節と同様 に，S 波部分にも各ケースで多少の違いが見られるが， より大きな違いが表面波部分で生じていることがわかる。 これは，前節での議論と同様，表面波の波長と本研究で 用いた粒子間隔の兼ね合いより変形勾配テンソルの近似 精度が低下したためであると考えられる。このことから， 表面波の伝播も含むような地震波動場のモデリングをお こなう際には，人工的な力の影響はより重要になると考 えられる。 4. まとめ 本研究では，HPM の地震波動場モデリングに対する適 Fig.9. The particle motion at the receiver point. Ellipsoidal motion can be observed. 用性を検討するため，実体波・表面波伝播のシミュレー ションを実施し，差分法や有限要素法の結果と比較する ことで考察した。以下に，本研究により得られた知見を 示す。 ・ HPM を地震波動場モデリングに適用する時，局所的 な振動を抑制する人工的な力を導入しない場合だと 結果に不自然な振動が観測されるため，人工的な力 ෡ = 0.25 × E (a) E の導入は必須であることがわかった。 ෡ は，媒質のヤン ・ 人工的な力の大きさを決定する係数E グ率 E の 0.5 倍程度とすることが望ましく，小さすぎ る，もしくは大きすぎる力を用いると結果に悪影響 をおよぼすことが示された。 ・ 振源周波数を高くしてシミュレーションをおこなっ たところ，より大きな数値振動が確認されたことか ら，影響半径内の変形が線形な状態から遠ざかると 変形勾配テンソルの近似精度が落ち，不自然な振動 ෡ = 0.5 × E (b) E が大きくなることが示された。また，この数値振動 の問題は粒子間隔を小さくすることで解決できるこ とを示した。 ・ 水平でない形状を有する地表面形状をモデル化し， そこを伝播する表面波のシミュレーションをおこな った。その結果，HPM により水平でない地表面上を 伝播する表面波が再現できることがわかった。また， ෡ = 1.0 × E (c) E 表面波伝播に関しては，実体波部分よりも数値振動 の影響を受けやすいことがわかった。 Fig.10. The waveform of the displacement velocity at the receiver point. HPM を用いる利点は，自由境界条件の導入が容易であ 物理探査 194 第 64 巻第 3 号 (2011) ることや，プリプロセスにおいて要素や節点のコネクテ Koshizuka, S. and Oka, Y. (1996) : Moving particle semi-implicit ィビティ作成を必要としないことなどが挙げられる。こ method for fragmentation of incompressible fluid, Nucl. Sci. れらの特徴は， 差分法や有限要素法に対する利点であり， Eng., 123, 421-434. 今後，これらの手法に変わる地震波動場モデリング手法 として利用されることが期待される。 Lucy, L. B. (1977) : A numerical approach to the testing of fission hypothesis, Astron. J., 82, 1013-1024. O’Brien, G. S. and Bean, C. J. (2004a) : A discrete numerical method 謝 辞 for modeling volcanic earthquake source mechanisms, J. 匿名の査読者からいただいた助言は，本稿を改善する 上で非常に有益でした。ここに記して感謝の意を表しま す。 Geophys. Res., 109, B09301. O’Brien, G. S. and Bean, C. J. (2004b) : A 3D numerical elastic lattice method for seismic wave propagation in heterogeneous media with topography, Geophys. Res. Lett., 31, L14608. 参 考 文 献 O’Brien, G. S., Bean, C. J. and Tapamo, H. (2009) : Dispersion analysis and computational efficiency of elastic lattice methods Chikazawa, Y., Koshizuka S. and Oka, Y. (2001) : A particle method for elastic and visco-plastic structures and fluid-structure interactions, Compt. Mech., 27, 97-106. Cundall, P. A. and Strack, O. D. L. (1979) : A discrete numerical model for granular assemblies, Geotechnique, 29, 47-65. Graves, R. W. (1996) : Simulating seismic wave propagation in 3D for seismic wave propagation, Computers & Geosciences, 35, 1768-1775. 鈴木幸人，越塚誠一，岡芳明（2005） ：HMPS (Hamiltonian Moving Particle Semi-implicit) 法の開発（第１報，運動方程式の導出） ， Transactions of JSCES，Paper No.20050016. Suzuki, Y., Koshizuka, S. and Oka, Y. (2007) : Hamiltonian elastic media using staggered-grid finite differences, Bull. Seism. moving-particle Soc. Am., 86, 1091-1106. incompressible fluid flows, Comput. Meth. Appl. Mech. Eng., 久田俊明，野口裕久（1995）：非線形有限要素法の基礎と応用， 丸善. 近藤雅裕，鈴木幸人，越塚誠一（2007） ：最小自乗近似による粒 子法弾性解析手法の振動抑制，Transactions of JSCES，Paper No.20070031. Koketsu, K., Fujiwara, H. and Ikegami, Y. (2004) : Finite-element simulation of seismic ground motion with a voxel mesh, Pure Appl. Geophys., 161, 2183-2198. Kondo, M., Suzuki, Y. and Koshizuka, S. (2010) : Suppressing local particle oscillations in the Hamiltonian particle method for elasticity, Int. J. Num. Meth. Eng., 81, 1514-1528. semi-implicit (HMPS) method for 196, 2876-2894. Suzuki, Y. and Koshizuka, S. (2008) : A Hamiltonian particle method for non-linear elastodynamics, Int. J. Num. Meth. Eng., 74, 1344-1373. Toomey, A. and Bean, C. J. (2000) : Numerical simulation of seismic waves using a discrete particle scheme, Geophys. J. Int., 141, 595-604. Valle-Garcia, R. D. and Sanchez-Sesma, F. J. (2003) : Rayleigh waves modeling using an elastic lattice model, Geophys. Res. Lett., 30, 1866. 武川・三ケ田・後藤：地震波動伝播シミュレーションにおける粒子法の適用性に関する研究 195 Applicability of a Hamiltonian particle method for the simulation of the seismic wave propagation Junichi Takekawa*, Hitoshi Mikada* and Tada-nori Goto* ABSTRACT It is important to simulate a strong ground motion induced by the occurrence of earthquakes for disaster mitigation and prediction. Though finite difference method (FDM) has been used to simulate seismic wave propagation, it is difficult to introduce traction-free boundary conditions in a model with an arbitrary ground surface. Since finite element method (FEM) has an advantage of introduction of traction-free boundary conditions, FEM was applied to simulate seismic wave propagation. However, the connectivity between elements and nodes are needed in pre-processing for FEM analyses. This will lead the complex data structure and time consuming process. On the other hand, many particle methods have been developed and applied to simulate solid analyses. A Hamiltonian particle method (HPM) is one of the particle methods and developed for accurate energy conserving method. In HPM, the traction-free boundary condition is automatically introduced like FEM. Furthermore, the data structure becomes very simple because the positions of particles are only needed in HPM. However, artificial forces are needed for suppressing local particle oscillations in HPM. In the present study, we apply HPM to simulate seismic wave propagation and evaluate the effect of artificial force with comparing to the results from FDM and FEM. The results show that HPM can reproduce seismic wave propagation with sufficient accuracy. Keywords: particle method, numerical simulation, seismic wave, surface wave Manuscript received May 9, 2011; Accepted August 10, 2011. * Kyoto University Kyotodaigaku-Katsura, Nishikyo-ku, Kyoto, 615-8540, Japan © 2011 SEGJ