地震波伝播数値シミュレーション講習会 2025年度 初級編

地震波伝播 数値シミュレーション講習会 2025年度 初級編 前田 拓人 弘前大学大学院理工学研究科

弾性波動論の基礎 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

問題設定：地震波動の数値シミュレーションとは？ p 構成関係式（ここでは弾性体） p 運動方程式 " ! #" "&"" "&"# "&"$ ! = + + + +" "$ ! "' ") "* " ! ## "&"# "&## "&#$ ! = + + + +# "$ ! "' ") "* " ! #$ "&"$ "&#$ "&$$ ! = + + + +% ! "$ "' ") "* !! , !" , !# 弾性変位 #!! , #"" , ### 法線応力 #"# , #!# , #!" 剪断応力 $ 2025-10-30 質量密度 !!! = # + 2& '!! + # '"" + '## !"" = # + 2& '"" + # '!! + '## !## = # + 2& '## + # '!! + '"" !"# = 2&'"# , !!# = 2&'!# , !!" = 2&'!" '!! , '"" , '## 法線ひずみ '"# , '!# , '!" 剪断ひずみ %, & Laméの定数 p これらが地震波動伝播を記述する偏微分方程式 p 解析解は限られた場合にしか存在しない p 離散化された位置と時間 (+$ , ,% , -& , .' ) について 近似的に微分方程式を解くことが目的となる 地震波伝播数値シミュレーション講習会 2025年度初級編

弾性体の運動方程式 in a nutshell for 1D p 連続体の仮想的な断面には単位面積あたり応力 (stress) ! " p 面の外側向き（膨張する向き）を正符号にとる p 長さ $" 面積 $% の棒の運動方程式は，質量密度を & として 3(4 0121+ = ! + + 1+ 12 − ! + 12 ( 3. （質量 x 加速度 = 力） *( p 両辺を $%$" で割って，$" → 0 の極限 3 ( 4 ! + + 1+ − ! + 3! 0 ( = → 3. 1+ 3+ p 単位体積あたりの運動方程式には，応力の 空間勾配が内力として働く 2025-10-30 # ( + *( # ( *+ ( ( + *( 地震波伝播数値シミュレーション講習会 2025年度初級編

ひずみと構成関係式 in a nutshell for 1D p 棒に応力がかかると，それに応じて変形：長さ $" が $" + $* に膨張 p ひずみ (strain) ：単位距離あたりの変形率 1+ + 17 − 1+ 17 '= = 1+ 1+ p 座標位置 " における変位で書くと 4 + + 1+ − 4(+) 34 '= = 1+ 3+ p 実験によると変形が微小なうちは，ひずみと 応力が比例関係にある p ! = ,- *( → *( + *! ( # ( p この前提のモデルが線形弾性体 # ( + *( *+ ( 2025-10-30 ! ( + *( ( + *( 地震波伝播数値シミュレーション講習会 2025年度初級編

運動方程式から波動方程式 34 3 ( 4 3! p 1次元の運動方程式 0 ( = , 構成関係式 ! = 8' , ひずみの定義 ' = 3+ 3. 3+ p 全部まとめて整理すると… 3 ( 4 1 3! 1 3 1 3 34 = = 8' = 8 ( 3. 0 3+ 0 3+ 0 3+ 3+ p 3(4 1 3 34 = 8 3. ( 0 3+ 3+ : （不均質媒質中の）波動方程式 p さまざまな物理学の分野で，時空間にまたがって情報を運ぶ波動現象 を記述する偏微分方程式（の一形態） p 弾性体の運動方程式と構成関係式は波動現象を表現する！ 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

1次元波動方程式
p

1 .2(", 5)
15 .

1 1
12 ", 5
=
,
& 1"
1"

,が場所によらない場合は

平面波解 2 ", 5 = ; exp ?(@" − A5)
" 5
= ; exp 2C? −
D E

1 .2(", 5) , 1 .2 ", 5
=
.
15
& 1" .
p 位相速度 . =

,/& の波動伝播

p dʼAlembert の解
2 ", 5 = 8 " − .5 + :(" + .5)
" 軸の ± 方向への平行移動
p （1次元では幾何減衰はない）
2025-10-30

! = 3 km/s

地震波伝播数値シミュレーション講習会 2025年度初級編

速度・応力型方程式への変換 p 1次元の方程式 1 .2(", 5) 1! ", 5 12 ", 5 & = , ! ", 5 = , . 15 1" 1" p 変位速度 F ", 5 = 12(", 5)/15 を使って方程式を記述（速度・応力型） 1F 1 1! 1! 1F = , =, 15 & 1" 15 1" 左辺に時間微分・右辺に空間微分の共通の構造 → 計算機で解きやすい p 3次元でも同様（抜粋） 1F/ 1!// 1!/0 1!/1 & = + + 15 1" 1G 1H 1F0 1F1 1!// 1F/ = D + 2I +D + , 15 1" 1G 1H 2025-10-30 1!/1 1F1 1F/ =I + 15 1" 1H 地震波伝播数値シミュレーション講習会 2025年度初級編

連続→離散へ p 1次元の運動方程式と構成関係式 1F 1 1! 1! 1F = , =, 15 & 1" 15 1" p 質量密度 & と弾性係数 , は，固体地球部の物性の地球物理学的指標 p 一般には場所に依存（不均質性） & = &("), , = , " p そうすると，微分方程式に解析解がない！（既知の数式で解けない） p 運動方程式の極限をとる前のモデルでは？ 3 ( 4 ! + + 1+ − ! + 0 ( = 3. 1+ p " と " + $" における値をばらばらに 解けば，微分しなくても良いのでは…？ 2025-10-30 *( # ( + *( # ( *+ ( ( + *( 地震波伝播数値シミュレーション講習会 2025年度初級編

差分近似と離散化 p 空間の微分を2点の値の「差」で表現する , Δ" Δ" 12 ", 5 2 "+ ,5 − 2 " − ,5 ! ", 5 = , ⟹ ! ", 5 = Δ" 2 2 1" p 連続で微分可能な関数 8(") の導関数の定義 $8 8 " + Δ" − 8 " 8 " + Δ"/2 − 8 " − Δ"/2 = lim = lim 2/→3 2/→3 $" Δ" Δ" p 微分の数学的定義から Δ" ⟶ 0 の極限を外し，有限の Δ" で代替 p 差分近似 Finite Difference Approximation p この近似を用いた近似解法を 差分法 Finite Difference Method という p 方程式の求解を Δ" 刻みの点群に限定する，ということでもある（離散化） 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

差分法の基本概念 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

差分近似とその精度 p 微分の数学的定義から Δ" ⟶ 0 の極限を外し，十分に小さい有限の Δ" で代替 p どれくらい小さくすれば十分なのか？…Taylor級数展開から検討 1 $.8 1 $48 7 58 $8 1 $ .+ 4+⋯= S 8 " + Δ" = 8 " + Δ" + Δ" Δ" Δ" . 4 5 $" 2! $" 3! $" T! $" 563 8 " + Δ" − 8 " $8 1 $.8 $8 ⟹ = + Δ" + ⋯ = + U(Δ") . Δ" $" 2! $" $" ↑「オーダー記号」 Δ(の1乗に比例する項の意 $8 8 " + Δ" − 8 " p は導関数 の近似評価で，Δ" に比例した誤差を持つ $" Δ" p Δ" を小さくしていけば，もとの微分に帰着することが保証された！ 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

3種類の差分 p 連続で微分可能な関数 8 に対しては，以下は全部等価 $8 8 " + Δ" − 8 " 8 " − 8 " − Δ" = lim = lim 2/→3 $" 2/→3 Δ" Δ" 8 " + Δ" − 8 " − Δ" 8 " + Δ"/2 − 8 " − Δ"/2 = lim = lim 2/→3 2/→3 2Δ" Δ" p しかし，有限な Δ" に対してはこれらの 値は全部違う！ p それぞれ前進差分・後退差分・中心差分 と呼ぶ 分野によっては「風上差分」「風下差分」と呼ぶことも 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

精度評価 p 前進・後退・中央差分で結果にどれだけ違いがあるか？ $8 1 $.8 . 1 $48 4 9 8 " ± Δ" = 8 " ± Δ" + Δ" ± Δ" + U Δ" $" 2! $" . 3! $" 4 前進 後退 中心 8 " + Δ" − 8 " $8 1 $.8 $8 = + Δ" + ⋯ = + U(Δ") . Δ" $" 2 $" $" 8 " − 8 " − Δ" $8 1 $.8 $8 = − Δ" … = + U(Δ") . Δ" $" 2 $" $" 8 " + Δ" − 8 " − Δ" $8 1 $48 . $8 .) = − Δ" … = + U(Δ" 2Δ" $" 6 $" 4 $" p 中心差分は Δ" . に比例する誤差・・・系統的に精度が高い 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

スタガードグリッド差分 p 数値シミュレーションでは，あらゆる " について微分方程式を解くかわりに， 離散化された ": = XΔ" (X = 1,2, ⋯ ) における値を求める： 8(:<=) − 8(:?=) $8 8 XΔ" + Δ" − 8 XΔ" − Δ" Y ≃ = $" /6:2/ 2Δ" 2Δ" p スタガードグリッド（食い違い格子） $8 8 " + Δ"/2 − 8 " − Δ"/2 Y ≃ $" /6:2/ Δ" 8(:<=/.) − 8(:?=/.) = Δ" 使われていない！ @+1 @−1 ( 2Δ( @ − 1/2 @ + 1/2 ( Δ( p 中心差分より精度は上がるが，8 と 1/ 8 の位置がズレる 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

高次精度差分 p Taylor展開をより高次まで実施し，高次項も0になるように公式を作る p 1階の導関数だけを残すために，より多くの情報が必要．以下は4次公式 3Δ+ 3 1: 9 1( : ( 9 1) : ) 27 1 * : * + : +− =: + − Δ+ + Δ+ − Δ+ + Δ+ + A Δ+ 2 2 1+ 8 1+ ( 16 1+ ) 128 1+ * Δ+ 1 1: 1 1( : ( 1 1) : ) 1 1* : * + : +− =: + − Δ+ + Δ+ − Δ+ + Δ+ + A Δ+ 2 2 1+ 8 1+ ( 48 1+ ) 384 1+ * Δ+ 1 1: 1 1( : ( 1 1) : ) 1 1* : * + : ++ =: + + Δ+ + Δ+ + Δ+ + Δ+ + A Δ+ 2 2 1+ 8 3+ ( 48 3+ ) 384 1+ * 3Δ+ 3 1: 9 1( : ( 9 1) : ) 27 1 * : * + : ++ =: + + Δ+ + Δ+ + Δ+ + Δ+ + A Δ+ 2 2 1+ 8 1+ ( 16 1+ ) 128 1+ * 1: 1 9 Δ+ Δ+ = : ++ −: +− 1+ Δ+ 8 2 2 1 3Δ+ 3Δ+ − : ++ −: +− 24 2 2 + A Δ+ * ※ E Δ( $ /Δ( = E Δ( % 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

高次精度差分 p 一般には，[次精度のスタガードグリッド差分（ただし[ は2の倍数）は //( 1: 1 1 1 / = C D, : + + E − Δ+ − : + − E − Δ+ 1+ Δ+ 2 2 ,-. と表される． p ただし， \GH は差分計算のための係数．たとえば 2次公式 D. = 1 4次公式 * D. ( 9 * 1 = , D( = − 8 24 （参）16次公式の一部 41409225 = 33554432 143 &' ,( = − 167772160 &' ,& 7桁にわたる変動 …計算機上の精度に限界 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

差分方程式の解法 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

空間微分の離散化 p 1次元の運動方程式と構成関係式 1F 1 1! 1! 1F = , =, をそのまま差分化 15 & 1" 15 1" 3!(+, .) F + + Δ+/2, . − F(+ − Δ+/2, .) = 8(+) 3. Δ+ 3F(+, .) 1 ! + + Δ+/2, . − !(+ − Δ+/2, .) = 3. 0(+) Δ+ p これでは，2本の式の F と ! の位置が整合的ではない！ p 構成関係式全体を Δ"/2 平行移動 3!(+ − Δ+/2, .) F + + Δ+, . − F(+, .) = 8(+ − Δ+/2) 3. Δ+ 3F(+, .) 1 ! + + Δ+/2, . − !(+ − Δ+/2, .) = 3. 0 Δ+ 2025-10-30 変位速度と応力が 互い違いに配置される 地震波伝播数値シミュレーション講習会 2025年度初級編

時間差分 p 時間についての微分も，同じように中心差分（間隔 Δ5）で近似 p 左辺の時間差分の中心時刻と右辺の時刻が一致 1 Δ+ Δ. Δ+ Δ. ! +− ,. + −! +− ,. − Δ. 2 2 2 2 時間についても互い違い Δ+ F + + Δ+, . − F(+, .) =8 +− 2 Δ+ F +, . + Δ. − F +, . 1 1 H+ Δ. H+ Δ. = ! ++ ,. + −! +− ,. + Δ. 0 + Δ+ 2 2 2 2 p もっとも未来の時間の情報について解くと，時間発展方程式を得る． Δ+ Δ. Δ+ Δ. 8(+ − Δ+/2) ! +− ,. + =! +− ,. + + F +, . − F + − Δ+, . Δ. 2 2 2 2 Δ+ 1 1 H+ Δ. H+ Δ. F +, . + Δ. = F +, . + ! ++ ,. + −! +− ,. + Δ. 0(+) Δ+ 2 2 2 2 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

グリッド番号によるコンパクト表現 p 8 " = XΔ", 5 = ]Δ5 = 8(:) と略記 (I) Δ+ Δ. Δ+ Δ. 8(+ − Δ+/2) ! +− ,. + =! +− ,. + + F +, . − F + − Δ+, . Δ. 2 2 2 2 Δ+ 1 1 H+ Δ. H+ Δ. F +, . + Δ. = F +, . + ! ++ ,. + −! +− ,. + 0 Δ+ 2 2 2 2 (5) (56./() (53./() !(23./() = !(23./() + 8(23./() (56./() 5 F(2) − F 23. Δ+ (56./() 1 !(26./() − !(23./() (56.) (5) F(2) = F(2) + Δ. 0(2) Δ+ 2025-10-30 Δ. Δ. 見やすくなるだけでなく コンピュータ・プログラムとの 親和性が Up 地震波伝播数値シミュレーション講習会 2025年度初級編

時間発展のようす (56./() (53./() !(23./() = !(23./() (5) +8(23./() (56.) F(2) 5 F(2) − F 23. Δ+ (5) = F(2) (56./() + Δ. (56./() 1 !(26./() − !(23./() 0(2) Δ+ Δ. p 周辺情報を用いながら 過去から未来に発展 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

時間発展の陽と陰 時間中心差分 時間後退差分 J K + ΔK − J(K) ΔK =L K+ ΔK 2 J K + ΔK − J(K) = L(K + ΔK) ΔK p 時間微分を中心差分で近似した結果は未来について陽に解ける（陽的） Δ+ Δ. Δ+ Δ. 8(+ − Δ+/2) ! +− ,. + =! +− ,. + + F +, . − F + − Δ+, . Δ. 2 2 2 2 Δ+ 1 1 H+ Δ. H+ Δ. F +, . + Δ. = F +, . + ! ++ ,. + −! +− ,. + Δ. 0(+) Δ+ 2 2 2 2 p もし後退差分なら？…全空間グリッドの連立方程式になってしまう（陰的） Δ+ Δ+ 8(+ − Δ+/2) ! +− , . + Δ. = ! + − ,. + F +, . + Δ. − F + − Δ+, . + Δ. Δ. 2 2 Δ+ 1 1 H+ H+ F +, . + Δ. = F +, . + ! ++ , . + Δ. − ! + − , . + Δ. Δ. 0(+) Δ+ 2 2 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

コンパクト表現：グリッド vs セル (.) (./&/!) (.+&/!) &(*+&/!) = &(*+&/!) + 7(*+&/!) . 8(*) − 8 *+& Δ' (./&/!) Δ$ (./&/!) 1 &(*/&/!) − &(*+&/!) (./&) (.) 8(*) = 8(*) + Δ$ !(*) Δ' p ½ グリッドはコンピュータ上のリストや配列では表現不能 p X − 1 Δ" < " ≤ XΔ" の範囲 を ? 番目のセル（Cell）と定義．時間も同様 (1) (1/&) &(0) (1) = &(0) + 7(0) 1 8(0) − 8 0+& Δ' (1/&) (1/&) 8(0) (1) = 8(0) + ( Δ' Δ$ p 番号がすべて整数になり，計算機 プログラムにしやすい形になった 2025-10-30 I− 3 1 Δx I− Δx ( 2) 2) (I − 1) Δx IΔx I− 3 2 (1/&) 1 &(0/&) − &(0) !(0) Δ$ i−1 I−1 I− 1 2 i I I+ 1 3 Δx I+ Δx ( 2) 2) (I + 1) Δx (I + 2) Δx ( I+ 1 2 i+1 I+1 I+ 3 2 i+2 I+2 x Grid Cell 地震波伝播数値シミュレーション講習会 2025年度初級編

メモリ（記憶）領域とメモリ更新 (1) (1/&) (1) (1) Δ' (1/&) (1/&) & − & 1 (0/&) (0) &(0) = &(0) + 7(0) (1/&) = 8(0) + 8(0) 1 8(0) − 8 0+& !(0) Δ' Δ$ Δ$ p 応力 !(M) と 変位速度 F(M) (? = 1 ⋯ ]) をセットとして扱う (5) (5) p 時間ステップ T から T + 1 にすべて更新したら，T 番目の値は不要 p 過去の状態をすべて覚えておく必要はない！ p 実際の計算は同じ変数を更新するスタイル &(0) ← &(0) + 7(0) 2025-10-30 8(0) − 8 0+& Δ' Δ$ (= = 1 ⋯ ?) ⇒ 1 & 0/& − & 0 8(0) ← 8(0) + Δ$ (= = 1 ⋯ ?) !0 Δ' 地震波伝播数値シミュレーション講習会 2025年度初級編

模擬コード p 応力 ! ，変位速度 F ，剛性率 , ，質量密度 & をリスト（配列）で表現 p プログラミング言語の A=B は等式ではなく「AをBという値に設定する」の意 do i=1, N for i in range (1,N+1): &(0) = &(0) + 7(0) 8(0) − 8 0+& Δ' Δ$ for i in range (1,N+1): 1 & 0/& − & 0 8(0) = 8(0) + Δ$ !0 Δ' &(0) = &(0) + 7(0) 8(0) − 8 0+& Δ' Δ$ end do do i=1, N 1 & 0/& − & 0 8(0) = 8(0) + Δ$ !0 Δ' end do ※ Pythonのrangeによるループは通常 0 から始まるが，ここではFortranにあわせて1始まりとしている． 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

差分法の限界と分散関係式 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

差分法の限界（１） 150 x [km] p 右図は両方とも1次元波動のスタガード グリッド差分法シミュレーション 200 p 上側：ベル型の形状が一定の速度でx軸 の両側に伝わっていく 100 50 0 p これは波動方程式の期待通り p 波形の形が変わってしまった p どうしてこうなった？ 0 10 15 20 15 20 t [s] 200 150 100 50 0 0 2025-10-30 5 Δ( = 2.0 km distance [km] p 下側：上と同じシミュレーションを， 空間グリッド間隔を10倍にしたもの Δ( = 0.2 km 5 10 time [s] 地震波伝播数値シミュレーション講習会 2025年度初級編

差分法の限界（２） 150 x [km] p 右図は両方とも1次元波動のスタガード グリッド差分法シミュレーション 200 p 上側：ベル型の形状が一定の速度でx軸 の両側に伝わっていく 100 50 0 p これは波動方程式の期待通り p 計算結果が爆発（発散） p 地震波速度を大きくしても同様の現象 p どうしてこうなった！！？？ 0 10 15 20 15 20 t [s] 200 150 100 50 0 0 2025-10-30 5 ΔK = 0.05001 秒 distance [km] p 下側：上と同じシミュレーションで， 時間ステップを 0.02% 増やしたもの ΔK = 0.05 秒 5 10 time [s] 地震波伝播数値シミュレーション講習会 2025年度初級編

方程式の近似と解の近似 基礎方程式 求解 方程式を近似 （困難） 差分方程式 求解 真の解 差分解 解を近似？？ p 『運動方程式の近似解』 『運動方程式を近似後の数値解』とは限らない！ p 解析的に解ける均質媒質で解の比較から検討する 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

分散関係式 p シミュレーションが正しく「波動らしい」振る舞いをするなら，平面波 F = F7 exp L M+ − N. , ! = !7 exp L M+ − N. が解になるはず p 離散化された方程式に上記の解を代入して整理すると，以下の関係を得る AΔ5 . @Δ" sin 2 = b. 2 (b ≡ ,/&) (A) . . Δ5 Δ" 2 2 p 角周波数 A と波数 @ の関係式： 一般に分散関係式という sin. p もともとの波動の分散関係式は A. = b.@ . p (A) 2025-10-30 (B) (B) のとき，数値シミュレーションが波動現象の良い近似になっている 地震波伝播数値シミュレーション講習会 2025年度初級編

分散関係式のざっくり解析 p 一般に " ≪ 1 に対して sin " ≃ " NΔ. MΔ+ sin( 2 = R( 2 Δ. ( Δ+ ( 2 2 sin( p 差分法の分散関係式において AΔ5/2 ≪ 1 かつ @Δ"/2 ≪ 1 であれば， ( NΔ. NΔ. sin( 2 ≃ 2 (, = N Δ. ( Δ. ( 2 2 ( MΔ+ MΔ+ sin( 2 ≃ 2 ( = M Δ+ ( Δ+ ( 2 2 p 時間・空間グリッド間隔が十分に小さければ元の分散関係 A. = b.@ . に一致 p 差分方程式の解は波動方程式の解の近似になっている，と言えそう 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

安定条件 p 差分方程式の分散関係式を変形すると . AΔ5 Δ5/2 @Δ" . . . sin = b sin . 2 Δ"/2 2 AΔ5 sin 2 Δ5 @Δ" =b sin Δ" 2 ≤1 p bΔ5 ≤ Δ" が必須：この条件が満たされないと『波』としての解を持てない p 意味：e" は時間 e5 で波が伝播する距離よりも大きくなければならない p （参考）3次元M次精度空間微分公式の場合 H/. Δ5 ≤ iQRS S \GH G6= ?= 1 1 1 + .+ . . Δ" ΔG ΔH B234 媒質中最大速度 6 ,5 差分公式の係数 p 等方4次公式なら \=9 = 9/8, \.9 = −1/24 で Δ5 ≤ 0.495 Δℎ/ iQRS Δ' = Δ) = Δ* = Δℎ 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

波長条件 A 2 Δ5 @Δ" arcsin b sin p 分散関係式を変形すると， . = = @ @Δ5 Δ" 2 p 地震波速度が波数（波長）に依存する：周波数分散性（数値分散） p @Δ" ≪ 1すなわち D ≫ Δ" のとき 2 Δ5 @Δ" .= arcsin b sin @Δ5 Δ" 2 2 Δ5 @Δ" ≃ arcsin b ≃b @Δ5 Δ" 2 (sin " ≃ arcsin " ≃ " for " ≪ 1) p 意味：e" は再現したい波の波長よりも十分細かくなくてはならない p 1次元M次精度公式の場合 H/. 2 Δ5 1 [H] .= arcsin b S \G sin 2C s − @Δ" @Δ5 Δ" 2 G6= 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

数値分散の影響 p 波長と空間グリッド 間隔の比 Λ に依存 p 安定条件からの最大時 間ステップと時間ス テップの比にも依存 p 2次元以上では方向に も依存してしまう （数値異方性） p 波長は場所による p あらゆる場所で 条件を満たせるよう に設定するしかない 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

2次元・3次元問題 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

3次元問題と2次元断面：基礎方程式 p 運動方程式と構成関係式（速度・応力型）…3次元ではこれをそのまま解く VW& V#!# V#"# V### VW" V#!" V#"" V#"# VW! V#!! V#!" V#!# $ = + + , $ VK = V( + VX + VY , $ VK = V( + VX + VY VK V( VX VY V### VW# VW! VW" VW" VW# V#"" VW" V#!! VW! VW! VW# = % + 2& +% + = % + 2& +% + , = % + 2& +% + , VK VY V( VX VK V( VX VY VK VX V( VY V#"# VW" VW# =& + , VK VY VX V#!# VW! VW# =& + , VK VY V( V#!" VW! VW" =& + VK VX V( ( p 伝播について興味のある方向に " 軸を取る p 例えば震源・観測点を含むラディアル方向 X p このとき， G = 0 断面における G 方向の 物理量変化は小さい，とみなせる Y 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

2次元断面とP-SV/SH 分解 p "H 断面において 1/1G = 0 を要請すると，方程式系が2つの独立な部分に分解 VW! V#!! V#!# $ = + , VK V( VY $ VW& V#!# V### = + VK V( VY V#!! VW! VW# = % + 2& +% VK V( VY V### VW# VW! = % + 2& +% VK VY V( V#!# VW! VW# =& + VK VY V( VW" V#!" V#"# $ = + VK V( VY VW" V#!" VW" V#"# =& =& , VK VY VK V( 2025-10-30 P-SV SH 弾性変位ポテンシャル表現 0 0 [ = ∇] + ∇× 0 + ∇×∇× 0 ` _ ] (, Y , `((, Y), _ (, Y がそれぞれ P, SV, SHに対応 地震波伝播数値シミュレーション講習会 2025年度初級編

SH波とスタガードグリッド p すべてが中心差分になるように 変数がレイアウトされている p 確認してみよう 3F" 3!!" 3!"# 0 = + 3. 3+ 33!!" 3F" =& 3. 3+ 3!"# 3F" =& 3. 3- 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

P-SV波とスタガードグリッド p すべてが中心差分になるように 変数がレイアウトされている p 確認してみよう 3F! 3!!! 3!!# 0 = + 3. 3+ 33F8 3!!# 3!## 0 = + 3. 3+ 33!!! 3F! 3F# = # + 2& +# 3. 3+ 33!## 3F# 3F! = # + 2& +# 3. 33+ 3!!# 3F! 3F# =& + 3. 33+ 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

x y x 3次元計算におけるスタガードグリッドレイアウト y J∆y I∆x 図 3.1 3 次元問題のスタガードグリッドの配置 p 煩雑なためセル (?, u, @) の変数のみ表示 p P-SV と SH の断面を y方向に並べた形 (a) 2D P-SV ıii vz P-SV (b) 2D SH SH vx vy ıxy ıxz Voxel (I,J,K) vz K∆z σyz σxz z σii vx vy σxy ıyz (K-1)∆z (J-1)∆y p 軸対称異方性までは表現できるが 21弾性係数の異方性ではうまくいかない 2025-10-30 x y J∆y I∆x (I-1)∆x x z y 地震波伝播数値シミュレーション講習会 2025年度初級編

SH波方程式系の離散化：グリッド表記 (56./() (5) (53./() 3!!" 3F" =& 3. 3+ !!"(2,:3./() − !!"(2,:3./() 3!"# 3F" =& 3. 3- !"#(23./(,:) − !"#(23./(,:) Δ. (56./() = &(2,:3./() Δ+ (5) (53./() Δ. (5) F" 26./(,:3./( − F" 23./(,:3./( = &(23./(,:) (5) F" 23./(,:6./( − F" 23./(,:3./( Δ- 3F" 3!!" 3!"# 0 = + 3. 3+ 3(56.) 0(23./(,:3./() Δ. (56./() = 2025-10-30 (5) F"(23./(,:3./() − F"(23./(,:3./() (56./() !!"(2,:3./() − !!"(23.,:3./() Δ+ (56./() + (56./() !"#(23./(,:) − !"#(23./(,:3.) Δ- 地震波伝播数値シミュレーション講習会 2025年度初級編

SH波方程式系の離散化：セル表記 ('6.) (') 3!!" 3F" =& 3. 3+ !!"($,&) − !!"($,&) 3!"# 3F" =& 3. 3- !"#($,&) − !"#($,&) Δ. ('6.) (') = &($,&) (') Δ. (') F" $6.,& − F" $,& Δ+ (') = &($,&) (') F" $,&6. − F" $,& Δ- (", $ − 1) 3F" 3!!" 3!"# 0 = + 3. 3+ 3('6.) 0($,&) (') F"($,&) − F"($,&) 2025-10-30 (", $) (" + 1, $) Δ. ('6.) = (" − 1, $) ('6.) !!"($,&) − !!"($,&3.) Δ+ ('6.) + ('6.) !"#($,&) − !"#($,&3.) (", $ + 1) Δ地震波伝播数値シミュレーション講習会 2025年度初級編

SH波方程式系の離散化：時間発展（更新）表現 ('6.) (') !!"($,&) − !!"($,&) Δ. ('6.) 0($,&) (') = &($,&) (') F"($,&) − F"($,&) Δ. Δ. ('6.) Δ+ !"#($,&) ← !"#($,&) + &($,&) F"($,&) ← F"($,&) + Δ+ !!"($,&) − !!"($,&3.) !!"($,&) ← !!"($,&) + &($,&) (') !"#($,&) − !"#($,&) ('6.) = ('6.) (') F" $6.,& − F" $,& F" $6.,& − F" $,& Δ+ F" $,&6. − F" $,& ('6.) + (') = &($,&) (') F" $,&6. − F" $,& Δ- ('6.) !"#($,&) − !"#($,&3.) Δ- Δ. Δ. 1 Δ!!"($,&) − !!"($3.,&) 0 $,& Δ+ + !"#($,&) − !"#($,&3.) Δ- Δ. これでそのままプログラミングできる形式になった！ 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

地震波数値シミュレーションの要素技術 地表＋海底の境界条件 震源における 地震波の励起 波動伝播モデル 特に内部減衰 モデル外周部：吸収境界条件 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編

まとめ：より現実的な数値シミュレーションに向けて p 差分法による地震波数値シミュレーションの核心部分を解説しました： p スタガードグリッドによる離散化の方法を学びました p また，離散化にともなう限界と制約があることも学びました p 一方，まだまだ他にもたくさんの要素技術が必要となります： p 震源における地震波励起の方法 p 地表・海底における物理的境界条件とその実現方法 p 計算領域外周における吸収境界条件とその実現方法 p 内部減衰の実現方法と実装 p 大規模な数値シミュレーションを実現する並列計算 p （これらについては中級編で概略を取り扱います） 2025-10-30 地震波伝播数値シミュレーション講習会 2025年度初級編