概要
今回の記事では、地震波の伝播シミュレーション（断層滑りによる地震発生の2次元シミュレーション）について説明したいと思います。
土木の分野では、直接使えないシミュレーションだと思います。
基礎方程式
地震波伝播シミュレーションで用いる基礎方程式は、以下の線形弾性体の運動方程式を用いる。
\begin{align}
\frac{\partial u_i}{\partial t} &= v_i \\
\rho \frac{\partial v_i}{\partial t} &= \frac{\partial }{\partial x_j} \sigma_{ij}
\end{align}
$u$ は変位であり、外力（重力）は無視した。$\sigma$ は応力であり、ラメ定数を用いて
\begin{align}
\sigma_{ij} &= \lambda \varepsilon_{ll} \delta_{ij} + 2 \mu
\varepsilon_{ij}\\
\varepsilon_{ij} &= \frac{1}{2} \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_j} \right)
\end{align}
と書ける。応力は、対称テンソルである。
今回の記事では、2次元シミュレーションとする。基礎方程式の $x$ 成分を分解すると （$x,y$ と表記された添字については、アインシュタインの縮約表記を使わない）
\begin{align}
\rho \frac{\partial v_x}{\partial t} &= \frac{\partial }{\partial x} \sigma_{xx} + \frac{\partial }{\partial y} \sigma_{xy}\\
&= \frac{\partial }{\partial x}\left\{\lambda \varepsilon_{xx} + \lambda \varepsilon_{yy} + 2 \mu \varepsilon_{xx}\right\} + 2 \mu \frac{\partial }{\partial y} \varepsilon_{xy} \\
&= \lambda
\left(\frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_y}{\partial x \partial y} \right)+ 2 \mu \frac{\partial^2 u_x}{\partial x^2} + \mu\left(\frac{\partial^2 u_y}{\partial x \partial y} +\frac{\partial^2 u_x}{\partial y^2} \right)\\
&= (\lambda+\mu)\left\{ \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_y}{\partial x \partial y}\right\} + \mu \left(\frac{\partial^2 }{\partial x^2} +\frac{\partial^2 }{\partial y^2} \right) u_x
\end{align}
$y$ 成分も同様に
\begin{align}
\rho \frac{\partial v_y}{\partial t}
&= (\lambda+\mu)\left\{ \frac{\partial^2 u_x}{\partial x \partial y} + \frac{\partial^2 u_y}{\partial y^2}\right\} + \mu \left(\frac{\partial^2 }{\partial x^2} +\frac{\partial^2 }{\partial y^2} \right) u_y
\end{align}
と分解できる。数値計算においては、変位の２階微分を離散化しなくても、変位から応力を計算し、応力の1階微分を離散化して計算できる。今回の記事では、上記の式を離散化する。
空間方向の離散化
空間方向の離散化は、基本的に、中心差分を用いる。２変数 $x,y$ の関数 $f(x,y)$ を用意する。微小変化を $\Delta x , \Delta y $ として、関数 $f(x+\Delta x,y+\Delta y)$ を展開し
\begin{align}
f(x+\Delta x,y+\Delta y) &= f(x,y)+\left(\Delta x \frac{\partial}{\partial
x} +\Delta y \frac{\partial}{\partial y} \right)f(x,y) \\
&+ \frac{1}{2!} \left(\Delta x \frac{\partial}{\partial
x} +\Delta y \frac{\partial}{\partial y} \right)^2f(x,y) + \dots \\
&= f(x,y)+\frac{\partial f}{\partial x}\Delta x +\frac{\partial f}{\partial y} \Delta y \\
&+\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 +\frac{\partial^2 f}{\partial x \partial y}\Delta x \Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots
\end{align}
$２$次近似において、$x$ の$1$階微分は、
\begin{align}
f(x+\Delta x,y)-f(x-\Delta x,y) &= f(x,y)+\frac{\partial f}{\partial x}\Delta x
+\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 + \dots \\
&-f(x,y)+\frac{\partial f}{\partial x}\Delta x
-\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial f}{\partial x}\approx \frac{f(x+\Delta x,y)-f(x-\Delta x,y)}{2\Delta x }
\end{align}
$y$ の$1$階微分は、
\begin{align}
f(x,y+\Delta y)-f(x,y-\Delta y) &= f(x,y)+\frac{\partial f}{\partial y}\Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
&-f(x,y)+\frac{\partial f}{\partial y}\Delta y
-\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial f}{\partial y}\approx \frac{f(x,y+\Delta y)-f(x,y-\Delta y)}{2\Delta y }
\end{align}
$x$ の$２$階微分は、
\begin{align}
f(x+\Delta x,y)+f(x-\Delta x,y) &= f(x,y)+\frac{\partial f}{\partial x}\Delta x
+\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 + \dots \\
&+f(x,y)-\frac{\partial f}{\partial x}\Delta x
+\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial^2 f}{\partial x^2}\approx \frac{f(x+\Delta x,y)-2f(x,y)+f(x-\Delta x,y)}{\Delta x^2 }
\end{align}
$y$ の$２$階微分は、
\begin{align}
f(x,y+\Delta y)+f(x,y-\Delta y) &= f(x,y)+\frac{\partial f}{\partial y}\Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
&+f(x,y)-\frac{\partial f}{\partial y}\Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial^2 f}{\partial y^2}\approx \frac{f(x,y+\Delta y)-2f(x,y)+f(x,y-\Delta y)}{\Delta y^2 }
\end{align}
$x,y$ 微分は、
\begin{align}
&f(x+\Delta x,y+\Delta y)+f(x-\Delta x,y-\Delta y)-f(x+\Delta x,y-\Delta y)-f(x-\Delta x,y+\Delta y) \\
&= f(x,y)+\frac{\partial f}{\partial x}\Delta x +\frac{\partial f}{\partial y} \Delta y +\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 +\frac{\partial^2 f}{\partial x \partial y}\Delta x \Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
&+f(x,y)-\frac{\partial f}{\partial x}\Delta x -\frac{\partial f}{\partial y} \Delta y +\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 +\frac{\partial^2 f}{\partial x \partial y}\Delta x \Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
&-f(x,y)-\frac{\partial f}{\partial x}\Delta x +\frac{\partial f}{\partial y} \Delta y -\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 +\frac{\partial^2 f}{\partial x \partial y}\Delta x \Delta y
-\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
&-f(x,y)+\frac{\partial f}{\partial x}\Delta x -\frac{\partial f}{\partial y} \Delta y -\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}\Delta x^2 +\frac{\partial^2 f}{\partial x \partial y}\Delta x \Delta y
-\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial^2 f}{\partial x \partial y}\approx \frac{f(x+\Delta x,y+\Delta y)+f(x-\Delta x,y-\Delta y)-f(x+\Delta x,y-\Delta y)-f(x-\Delta x,y+\Delta y) }{4\Delta x\Delta y }
\end{align}
と求まる。
今回に記事では、境界条件を計算する際に$y$ の$1$階微分の計算が必要である。しかし、両側に格子点がなく中心差分が取れないので、片方に偏るように
\begin{align}
f(x,y-\Delta y) &=f(x,y)-\frac{\partial f}{\partial y}\Delta y
+\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
f(x,y-2\Delta y)&= f(x,y)-2\frac{\partial f}{\partial y}\Delta y
+2\frac{\partial^2 f}{\partial y^2}\Delta y^2 + \dots \\
\end{align}
\begin{align}
\therefore \frac{\partial f}{\partial y}\approx \frac{3f(x,y)-4f(x,y-\Delta y)+f(x,y+2\Delta y) }{2\Delta y }
\end{align}
を使う。
時間方向の離散化
時間方向の離散化は、前進差分を用いる。
\begin{align}
\frac{\partial u_i}{\partial t} &\approx \frac{u_i(t+\Delta t)-u_i(t)}{\Delta t} \\
\frac{\partial v_i}{\partial t} &\approx \frac{v_i(t+\Delta t)-v_i(t)}{\Delta t}
\end{align}
離散化された基礎方程式
上記の方法で、基礎方程式を離散化する。グリッドの $x,y$ 方向の格子点を $i,j$ に対応させ、時間ステップを $n$ とする。すると、速度の時間微分の項は、
\begin{align}
\frac{\partial v_x}{\partial t}
&= \frac{\lambda+\mu}{\rho}\left\{ \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_y}{\partial x \partial y}\right\} + \frac{\mu}{\rho} \left(\frac{\partial^2 }{\partial x^2} +\frac{\partial^2 }{\partial y^2} \right) u_x \\
\frac{\partial v_y}{\partial t}
&= \frac{\lambda+\mu}{\rho}\left\{ \frac{\partial^2 u_x}{\partial x \partial y} + \frac{\partial^2 u_y}{\partial y^2}\right\} + \frac{\mu}{\rho} \left(\frac{\partial^2 }{\partial x^2} +\frac{\partial^2 }{\partial y^2} \right) u_y
\end{align}
\begin{align}
\therefore \ \ & v_x^{n+1}(i,j) = v_x^{n}(i,j) + \frac{(\lambda+2\mu)\Delta t}{\rho \Delta x^2}\left\{u_x^n(i+1,j)-2u_x^n(i,j) +u_x^n(i-1,j) \right\} \\
& + \frac{(\lambda+\mu)\Delta t}{4\rho \Delta x \Delta y}\left\{u_y^n(i+1,j+1)-u_y^n(i+1,j-1)-u_y^n(i-1,j+1) +u_y^n(i-1,j-1) \right\} \\
& + \frac{\mu\Delta t}{\rho \Delta y^2}\left\{u_x^n(i,j+1)-2u_x^n(i,j) +u_x^n(i,j-1) \right\} \\
\therefore \ \ & v_y^{n+1}(i,j) = v_y^{n}(i,j) + \frac{(\lambda+2\mu)\Delta t}{\rho \Delta y^2}\left\{u_y^n(i,j+1)-2u_y^n(i,j) +u_y^n(i,j-1) \right\} \\
& + \frac{(\lambda+\mu)\Delta t}{4\rho \Delta x \Delta y}\left\{u_x^n(i+1,j+1)-u_x^n(i+1,j-1)-u_x^n(i-1,j+1) +u_x^n(i-1,j-1) \right\} \\
& + \frac{\mu\Delta t}{\rho \Delta x^2}\left\{u_y^n(i+1,j)-2u_y^n(i,j) +u_y^n(i-1,j) \right\} \\
\end{align}
変位の時間微分の項は、
\begin{align}
&\frac{\partial u_x}{\partial t} = v_x \ , \ \frac{\partial u_y}{\partial t} = v_y \\
\end{align}
\begin{align}
\therefore & u_x^{n+1}(i,j) = u_x^{n}(i,j) + v_x^{n+1}(i,j)\Delta t \\
\therefore & u_y^{n+1}(i,j) = u_y^{n}(i,j) + v_y^{n+1}(i,j)\Delta t
\end{align}
となる。
時間発展部分のコードは、こんな感じだろう。
!=====================================================================
! 支配方程式
!=====================================================================
subroutine elasticity_eq(n)
USE params
implicit none
integer:: i,k,n
! 境界条件の設定
call DirichletBoundaryCondition(n)
call NeumannBoundaryCondition(n)
! 波動方程式
do i = 1,nx-1
do k = 1,nz-1
! 震源域は更新しない
if (Epicenter_flag(i,k) == 1) then
! 震源域の計算
call Epicenter(n,i,k)
else
! 速度の更新
vx(n+1,i,k) = vx(n,i,k) + (dt*(lame+2.0*mu)/(rho*dx**2))*(ux(n,i+1,k)-2.0*ux(n,i,k)+ux(n,i-1,k)) &
& + (dt*(lame+mu)/(4.0*rho*dx*dz))*(uz(n,i+1,k+1)-uz(n,i+1,k-1)-uz(n,i-1,k+1)+uz(n,i-1,k-1)) &
& + (dt*mu/(rho*dz**2))*(ux(n,i,k+1)-2.0*ux(n,i,k)+ux(n,i,k-1))
vz(n+1,i,k) = vz(n,i,k)+ (dt*(lame+2.0*mu)/(rho*dz**2))*(uz(n,i,k+1)-2.0*uz(n,i,k)+uz(n,i,k-1)) &
& + (dt*(lame+mu)/(4.0*rho*dx*dz))*(ux(n,i+1,k+1)-ux(n,i+1,k-1)-ux(n,i-1,k+1)+ux(n,i-1,k-1)) &
& + (dt*mu/(rho*dx**2))*(uz(n,i+1,k)-2.0*uz(n,i,k)+uz(n,i-1,k))
! 変位の更新
ux(n+1,i,k) = ux(n,i,k) + dt*vx(n+1,i,k)
uz(n+1,i,k) = uz(n,i,k) + dt*vz(n+1,i,k)
endif
enddo
enddo
end subroutine
境界条件
今回の記事では、2次元のシミュレーションなので、４つの境界を定める必要がある。x方向の格子点 $i$ の範囲を $0,1,2,\cdots,n_x$ とし、x方向の格子点 $j$ の範囲を $0,1,2,\cdots,n_y$ とする。
左右の側面と底面は固定境界条件を課す。つまり、左右の側面の境界条件は、
\begin{align}
u_x^{n}(0,j) = 0 \ , \ u_x^{n}(n_x,j) =0 \\
u_y^{n}(0,j) = 0 \ , \ u_y^{n}(n_x,j)=0
\end{align}
底面の境界条件は、
\begin{align}
u_x^{n}(i,0) = 0 \ , \ u_y^{n}(i,0) =0 \\
\end{align}
と書ける。
上面は、地表に対応しており、地表付近の変位は最も知りたい値となる。地表では、地盤が自由に動けるように自由境界条件を課す。自由表面における応力条件は、
\begin{align}
\sigma_{xy} = 0 \ , \ \sigma_{yy} =0 \\
\end{align}
である。変位で表すと
\begin{align}
\sigma_{yy} &= \lambda \varepsilon_{xx} + \lambda \varepsilon_{yy} + 2 \mu \varepsilon_{yy}= \lambda \frac{\partial u_x}{\partial x } + (\lambda +2\mu)\frac{\partial u_y}{\partial y} =0
\\
\sigma_{xy} &= 2 \mu \varepsilon_{xy}=\mu\left(\frac{\partial u_y}{\partial x } +\frac{\partial u_x}{\partial y} \right)=0
\end{align}
となる。自由表面は、$j = n_y$ の部分なので、その部分だけ離散化すると、
\begin{align}
&\lambda\frac{u_x^{n}(i+1,n_y)-u_x^{n}(i-1,n_y) }{2\Delta x} +(\lambda +2\mu)
\frac{3u_y^{n}(i,n_y)-4u_y^{n}(i,n_y-1)+u_y^{n}(i,n_y-2) }{2\Delta y }=0
\\
&\frac{u_y^{n}(i+1,n_y)-u_y^{n}(i-1,n_y) }{2\Delta x} +\frac{3u_x^{n}(i,n_y)-4u_x^{n}(i,n_y-1)+u_x^{n}(i,n_y-2) }{2\Delta y }=0
\end{align}
となる。これらの方程式は、簡単には解けないので、連立方程式を解く必要がある。
上記の方程式を解く上で、簡単な方法として、ガウス＝ザイデル法や逐次加速緩和法がある。それらを適応するために、上記の離散化された方程式を、
\begin{align}
u_y^{n}(i,n_y) &= \frac{1}{3}\left\{ 4u_y^{n}(i,n_y-1)-u_y^{n}(i,n_y-2)
-\frac{\lambda}{\lambda +2\mu}\frac{\Delta y}{\Delta x}\left( u_x^{n}(i+1,n_y)-u_x^{n}(i-1,n_y)\right) \right\} \\
u_x^{n}(i,n_y) &= \frac{1}{3}\left\{ 4u_x^{n}(i,n_y-1)-u_x^{n}(i,n_y-2)
- \frac{\Delta y}{\Delta x}\left( u_y^{n}(i+1,n_y)-u_y^{n}(i-1,n_y)\right) \right\} \\
\end{align}
と変形する。逐次加速緩和法では、反復回数を $m$ として、
\begin{align}
\tilde{u}_y^{(m+1)n}(i,n_y) &= \frac{1}{3}\left\{4u_y^{(l)n}(i,n_y-1)-u_y^{(l)n}(i,n_y-2)
-\frac{\lambda}{\lambda +2\mu}\frac{\Delta y}{\Delta x}\left( u_x^{(l)n}(i+1,n_y)-u_x^{(l)n}(i-1,n_y)\right) \right\} \\
u_y^{(m+1)n}(i,n_y) &= u_y^{(m)n}(i,n_y) +\omega \left(\tilde{u}_y^{(m+1)n}(i,n_y) -u_y^{(m)n}(i,n_y) \right) \ \ , \ \ (0<\omega<2)
\end{align}
そして、
\begin{align}
\tilde{u}_x^{(m+1)n}(i,n_y) &= \frac{1}{3}\left\{ 4u_x^{(l)n}(i,n_y-1)-u_x^{(l)n}(i,n_y-2)
- \frac{\Delta y}{\Delta x}\left( u_y^{(l)n}(i+1,n_y)-u_y^{(l)n}(i-1,n_y)\right) \right\} \\
u_x^{(m+1)n}(i,n_y) &= u_x^{(m)n}(i,n_y) +\omega \left(\tilde{u}_x^{(m+1)n}(i,n_y) -u_x^{(m)n}(i,n_y) \right) \ \ , \ \ (0<\omega<2)
\end{align}
と更新する。$l$ は、更新の順番により、$u_x^{(m+1)n}(i,n_y),u_x^{(m+1)n}(i,n_y)$ より先に更新された場合は $l=m+1$ となり、他の場合は $l=m$ となる。$\omega = 1$ のとき、ガウス＝ザイデル法となる。
境界部分のコードは、こんな感じだろう。
!=====================================================================
! 境界条件
!=====================================================================
subroutine DirichletBoundaryCondition(n)
USE params
implicit none
integer,intent(in) :: n
integer:: i,k
! ディリクレ条件
! 底面の固定
do i = 0,nx
ux(n,i,0) = 0.0
uz(n,i,0) = 0.0
enddo
! 側面の固定
do k = 0,nz
! 左
ux(n,0,k) = 0.0
uz(n,0,k) = 0.0
! 右
ux(n,nx,k) = 0.0
uz(n,nx,k) = 0.0
enddo
end subroutine
!=====================================================================
! 境界条件
!=====================================================================
subroutine NeumannBoundaryCondition(n)
USE params
implicit none
integer,intent(in) :: n
integer:: i,iter
real(dp) :: norm
real(dp) :: uxm,uxb,uzm,uzb
integer,parameter :: itermax = 1000000
real(dp),parameter :: omega = 1.0
real(dp),parameter :: eps_sor = 1.0e-6
! SOR 法の繰り返し計算
iter = 0
do while (iter < itermax )
norm = 0.0
iter = iter + 1
do i = 1,nx-1
! sigma_xy = 0
uxm = (1.0/3.0)*(4.0*ux(n,i,nz-1) - ux(n,i,nz-2) - (dz/dx)*( uz(n,i+1,nz)-uz(n,i-1,nz) ) )
uxb = ux(n,i,nz)
ux(n,i,nz) = (1.0-omega)*uxb + omega*uxm
norm = norm + (ux(n,i,nz)-uxb)**2
! sigma_yy = 0
uzm = (1.0/3.0)*(4.0*uz(n,i,nz-1) - uz(n,i,nz-2) - (lame/(lame+2.0*mu) )*(dz/dx)*( ux(n,i+1,nz)-ux(n,i-1,nz) ) )
uzb = uz(n,i,nz)
uz(n,i,nz) = (1.0-omega)*uzb + omega*uzm
norm = norm + (uz(n,i,nz)-uzb)**2
enddo
! 誤差が小さいなら終了
if ( norm < eps_sor ) exit
enddo
end subroutine
震源域
以上の計算を行えば、地震動のシミュレーションができるが、地震を発生させなければならない。地震の発生は、断層が滑ることで起こるように設定する。
断層が一本の直線で表されたとし、両端の座標が $(x_0,y_0)$ と $(x_1,y_1)$ とする。座標変換を行い、新たな変数 $\xi,\eta$ を
\begin{align}
\xi =& (x-x_0)\cos \theta + (y-y_0)\sin \theta \\
\eta =& (y-y_0)\cos \theta - (x-x_0)\sin \theta \\
\end{align}
とする。また、$\theta ,r $ は、
\begin{align}
x_1-x_0 &= r \cos \theta \ \ , \ \ y_1-y_0 = r \sin \theta \\
r &= \sqrt{(x_1-x_0)^2 + (y_1-y_0)^2}
\end{align}
から計算する。そして、 断層の幅 $2d$ とするならば、$\xi,\eta$ は、
\begin{align}
0 \leq \xi \leq r \ \ , \ \ -d \leq \eta \leq d
\end{align}
を満たすように震源域を設定すれば、一本の直線の断層が設定できる。
震源地に入る格子点において変位を、
\begin{align}
u_x = \frac{U\eta}{2d} \cos\theta \ \ , \ \ u_y = \frac{U\eta}{2d} \sin\theta
\end{align}
と設定する。断層滑り $U$ は、破壊が瞬間的に伝播するとして、
\begin{align}
U =
\begin{cases}
U = 0 & (t< 0) \\
U = U_0\frac{t}{T_0} & (0 < t < T_0 ) \\
U = U_0 & (t > T_0 )
\end{cases}
\end{align}
断層面上で一様に生じるとする。$U_0$ は最終的な滑り量の大きさ、$T_0$ は滑りにかかる時間である。
コードは、こんな感じだろう。
!=====================================================================
! 初期パラメータの設定
!=====================================================================
SUBROUTINE Init_param()
USE params
implicit none
integer:: n
integer:: i,k
real(dp) :: eta,xi
real(dp) :: x,z
! ラメ定数
mu = beta**2*rho
lame = alpha**2*rho-2.0*mu
! 震源域の計算に必要なパラメータ
r = sqrt((x1-x0)**2 +(z1-z0)**2 )
sin_t = (z1-z0)/r
cos_t = (x1-x0)/r
! 初期化
ux(:,:,:) = 0.0
uz(:,:,:) = 0.0
vx(:,:,:) = 0.0
vz(:,:,:) = 0.0
! 震源域のフラグ
Epicenter_flag(:,:) = 0
! 震源域の設定
n=0
do i = 0,nx
do k = 0,nz
x = i*dx
z = k*dz
xi = (x-x0)*cos_t+(z-z0)*sin_t
eta = (z-z0)*cos_t-(x-x0)*sin_t
! 震源域の範囲
if ((0 <= xi).and.(xi <= r) ) then
if ((-d <= eta).and.(eta <= d) ) then
Epicenter_flag(i,k) = 1
call Epicenter(n,i,k)
endif
endif
enddo
enddo
end subroutine
!=====================================================================
! 震源域
!=====================================================================
SUBROUTINE Epicenter(n,i,k)
USE params
implicit none
integer,intent(in) :: n,i,k
real(dp) :: eta, x,z
x = i*dx
z = k*dz
eta =(z-z0)*cos_t-(x-x0)*sin_t
! 時間ごとに滑り量を定義
if (n*dt <= 0.0 ) then
vx(n,i,k) = 0.0
vz(n,i,k) = 0.0
ux(n,i,k) = 0.0
uz(n,i,k) = 0.0
elseif ( (0.0 < n*dt).and.(n*dt <= T0) ) then
vx(n,i,k) = (U0/T0)*eta*cos_t/(2.0*d)
vz(n,i,k) = (U0/T0)*eta*sin_t/(2.0*d)
ux(n,i,k) = (U0*n*dt/T0)*eta*cos_t/(2.0*d)
uz(n,i,k) = (U0*n*dt/T0)*eta*sin_t/(2.0*d)
elseif (T0 < n*dt ) then
vx(n,i,k) = 0.0
vz(n,i,k) = 0.0
ux(n,i,k) = U0*eta*cos_t/(2.0*d)
uz(n,i,k) = U0*eta*sin_t/(2.0*d)
endif
end subroutine
計算結果
シミュレーションの計算領域は、横幅 $100\mbox{km}$ 、深さ $50\mbox{km}$ とし、グリッドのセルサイズは $\Delta x, \Delta y = 1\mbox{km} $ とする。
地震波の速度は、
\begin{align}
\alpha &= \sqrt{\frac{\lambda+2\mu}{\rho}} = 6\mbox{km} \\
\beta &= \sqrt{\frac{\mu}{\rho}} = 3.5\mbox{km}
\end{align}
とかけ、$\alpha$ は縦波（P波）に対応し、$\beta$ は横波（S波）に対応する。計算するときは、
\begin{align}
\mu &= \beta^2\rho \\
\lambda &= \alpha^2 \rho-2\mu
\end{align}
とする。
波動方程式における安定性条件は、
\begin{align}
\frac{\alpha \Delta t}{\Delta x} \leq 1
\end{align}
なので、これを満たすように、
\begin{align}
\Delta t = 0.01 < \frac{\Delta x}{\alpha} = 0.1666\ldots
\end{align}
と設定する。
震源域は、両端の座標を $(x_0,y_0)=(47\mbox{km},22\mbox{km})$ と $(x_1,y_1)=(53\mbox{km},28\mbox{km})$ 、最終的な滑り量の大きさは $U_0 = 0.0004\mbox{km}$ 、滑りにかかる時間 $T_0=0.5\mbox{s}$ と設定する。そして、断層の幅は $d=1\mbox{km}$ とした。
シミュレーション結果は、以下である。
Density と書いてあるが、ヒートマップは密度ではなく、x方向の変位である（3dplot 形式で保存してしまったため）。震源域から、断層滑りが生じ、波が伝播していく様子が確認できる。
最後に
今回の記事では、地震波の伝播シミュレーションについて説明をした。 今回の記事では、一般格子を用いているが、精度を向上させるには、スタガード格子を用いる方が良いと考えられる。
今後、機会があれば、今回のシミュレーションとデータ同化（4次元変分法）を使い、ラメ定数$\lambda,\mu$ などの物性値を推定できるか試してみたい。（今回の離散化の方法を使うとアジョイント方程式は簡単に求まる。）
参考文献
記事は、以下の書籍を参考にしています。
自然災害のシミュレーション入門
井田 喜明(著)