概要

　今回の記事では、地震波の伝播シミュレーション（断層滑りによる地震発生の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$ などの物性値を推定できるか試してみたい。（今回の離散化の方法を使うとアジョイント方程式は簡単に求まる。）

参考文献

記事は、以下の書籍を参考にしています。

自然災害のシミュレーション入門
井田喜明(著)