2次元波動方程式の解の導出と性質

#Abstract

この記事は, 数理物理 Advent Calendar 2018 に空白があったため, 遅ればせながらお邪魔するものです. 放射条件を伴なう2次元全空間の波動方程式

\begin{align}
\partial_t^2 G - \beta^2 \left( \partial_x^2 + \partial_z^2 \right) G & = \delta(x) \delta(z) \delta(t), & (t \ge 0) \\
G \left( x,z,t \right) & = 0, & (\left| x \right|, \left| z \right| \rightarrow \infty )\\
G \left(x,z,t \right) = \partial_t G \left(x,z,t \right) & = 0 & ( t \lt 0 )
\end{align}

の基本解 $G$ を求める Cagniard-De Hoop 法を紹介し, その発展版として2層媒質中での波動場について部分的に述べます.

2次元全空間の波動方程式の基本解は, まず3次元全空間での基本解を求めておいて, その上で一様波源が1直線上に並んでいると仮定し解を1方向に積分すれば得ることができます. しかし別の方法として, 2次元の方程式に対し時間と空間それぞれの変数について Laplace 変換や Fourier 変換を施して常微分方程式に帰着し, これを解いた後に逆変換する Cagniard-De Hoop 法というものがあります. この手法の利点は, 一様全空間ではなく_2層媒質のような構造中における反射波・屈折波をも取り扱い可能である_ という点です. つまり解の対称性を仮定する必要がありません.

本稿では, Cagniard-De Hoop 法の要点を掴むために均質な全空間での解の導出を(pdfファイルにて)紹介した後, デモンストレーションとして半無限2層媒質中における波動場の数値解を紹介します. 解の導出と数値シミュレーションの節は独立なので, どちらか一方だけでも読んで頂ければ幸いです.

#Cagniard-De Hoop 法の紹介

数式ばかりなので以下の文章に纏めました:
Cagniard-De Hoop 法による2次元波動方程式の基本解の導出[pdf]

この内容を理解するための知識は, 主に複素関数論と Fourier 解析学, そして若干の超関数論です. 先の話を少しだけしておくと, $x$ 方向に Fourier 変換しておいて $z$ 方向には変換していないので, 例えば $z = -h$ のような定数で表わされる反射境界がある場合に, そこでの境界条件を Fourier 領域で合わせることができるような構造になっています. この性質を利用した発展版については, 上記 pdf ファイル中でも参照した Aki & Richards の教科書などに解説があります(が, 同書は媒質境界を透過する屈折波を求める手続きや, 次節で紹介する広角反射波の扱いについて記述が不足しているので, いずれ解説したいと思います).

#2層媒質中での波動場

水面や異種材料の接着面, あるいは地質学的構造境界など, 波の伝播速度の不連続面では, ご存知の通り波の反射と屈折が生じます. 反射波・屈折波の性質は, 境界に対し平面波が入射する場合には Snell の法則で簡潔に記述できますが, 点震源から放射された円筒波の場合にはそれほど単純ではありません. 今回は数値解を基に, その性質のひとつである Head wave (先頭波), あるいは Wide-angle reflection (広角反射波) と呼ばれる現象を紹介します.

4000x1000 grids, カラーマップ変更 pic.twitter.com/vanJEHCpqT

— S.Hirano (@Bimaterial) 2018年8月27日

上記 tweet 中の動画では, 領域の上下で波の速度が異なっており, 速度境界の若干上の点震源から波が放射されています. 領域上端は自由端で, ここに到達する波を観察していると, 動画序盤では震源からの直達波が真っ先にやって来ますが, 動画後半では震源から遠く離れた場所に対し, 速度境界からの波が一番最初に到達する様子が見えます. 遠方へは高速な下部領域経由で届いたほうが早いためです. 波は走時を極小化するような経路を見つける(= Fermat の原理, 正確には停留点ですが)というわけです. 数値計算によらなくとも Cagniard-De Hoop 法でこのような波動場が得られることは, またいずれ何かの機会に…

動画からも見て取れるように, 遠方で観測される Head Wave の振幅は直達波に比べるとかなり小さいのですが, それでもこの性質を利用して, 弾性体内部に存在する物質の不連続面を検出することができます. これは地球の内部構造推定や地下資源探査に有用な知識ですので, みなさんもこれを応用して石油王を目指してください.

最後に, 上記シミュレーションのための julia のコードを置いておきます. 領域内部では空間2次精度の Staggered Grid による差分法を, 上端以外の外部境界では Perfectly Matched Layer による吸収端を導入しています, julia を勉強し始めた頃に作成したもので, 試しに上付き/下付き文字を多用してみたのですが, 手間の割に可読性が良いかというと微妙な気もしますね… 描画は julia ではなく, バイナリファイルを出力しておいて gnuplot に任せました.

#=
   Finite difference solver for anti-plane wave propagation in a semi-infinite bimaterial system.

   Governing equation for x₁∈R, x₂<0, and t > 0:
   ρ ∂ₜ v = ∂x₁ σ₁₃ + ∂x₂ σ₂₃ + f
   ∂ₜ σ₁₃ = μ ∂x₁ v
   ∂ₜ σ₂₃ = μ ∂x₂ v
   where
   (ρ, μ): depth-dependent density and rigidity
   v: velocity in x₃-direction (= ∂ₜ u, where u is displacement in x₃-direction)
   (σ₁₃, σ₂₃) = μ (∂x₁ u, ∂x₂ u): 1,3- and 2,3-component of the shear stress
   f = F(t/T,r/R): body force in x₃-direction, where F(t,r) = t exp(-t) exp(-r²), and r² = x₁² + x₂²
   T and R are the scaling parameter in time and space, respectively

   Boundary conditions:
   σ₂₃ = 0 at x₂ = 0 (i.e., x₂ = 0 is a free surface)
   H: depth of a material interface (i.e., μ = μ⁺ for x₂ > -H and μ = μ⁻ for x₂ < -H)
   v, σ₁₃, σ₂₃→0 at far-field, which is achieved using Perfectly Matched Layer(PML)

   Also, this solver uses the following numerical parameters:
   L₁, L₂: length of the calculation domain along x₁- and x₂-axis, respectively
   N₁: number of grids along x₁-axis
   dx = L₁/N₁: grid interval (calculated)
   N₂ = L₂/dx: number of grids along x₂-axis (calculated)
   Nₜ: number of time steps
   CFL: Courant number
   η: maximum value of the numerical viscocity for PML
   Nₚₘₗ: number of grids for PML
=#

using SpecialFunctions
using Printf

const (L₁,L₂)=(4.0,1.0)   #   length of the domain
const (μ⁺,μ⁻)=(1.0,2.0)   #   rigidity
const (ρ⁺,ρ⁻)=(1.0,1.0)   #   density
const (c⁺,c⁻)=(√(μ⁺/ρ⁺),√(μ⁻/ρ⁻))   #   wave speed
const (x₁⁰,x₂⁰) = (0.5, -0.4) #   hypocenter
const H = -0.5*L₂  #   depth of the material interface
const T = 0.01  #   scaling parameter of the source in time
const R = 0.01  #   scaling parameter of the source in space

const N₁=1001  #   number of grids in x₁ direction
const Nₜ=2*(N₁-1)   #   number of time steps
const CFL=0.5   #   Courant number
const dx=L₁/(N₁-1.0)    #   grid interval
const dt=CFL*dx/max(c⁺,c⁻)  #   time interval
const x₁=0.0:dx:+L₁ #   x₁ coordinate
const x₂=-L₂:dx:0.0 #   x₂ coordinate
const N₂=length(x₂)  #   number of grids in x₂ direction
const t=0.0:dt:Nₜ*dt    #   time
const η = 0.2   #   maximum value of the numerical viscocity for PML
const Nₚₘₗ = 40 #   number of grids for PML

function FDTD()
    x₂ˢᵍ = Staggered(x₂)
    μ = Bimaterial(μ⁺,μ⁻,H,x₂)
    μˢᵍ = Bimaterial(μ⁺,μ⁻,H,x₂ˢᵍ)
    ρ = Bimaterial(ρ⁺,ρ⁻,H,x₂)
    (η₁,η₂,η₁₃,η₂₃) = Viscosity(x₁,x₂)
    (ρ⁻¹ᵣ,μᵣ,μᵣˢᵍ,s) = NumericalCoefficients(ρ,μ,μˢᵍ,x₁,x₂,x₁⁰,x₂⁰,R,dx,dt)

    (u,v,v₁,v₂,σ₁₃,σ₂₃) = (zeros(N₁,N₂),zeros(N₁,N₂),zeros(N₁,N₂),zeros(N₁,N₂),zeros(N₁-1,N₂-1),zeros(N₁-1,N₂-1))

    @time for i = 2:Nₜ
        STF = f(t[i],T)   #   source time function
        @inbounds for j = 2:N₁-1, k = 2:N₂-1
            u[j,k] = u[j,k] + 0.5*v[j,k]*dt #   time-integration of v w/ trapezoidal rule (+half of the previous v)
            BodyForce = STF*s[j,k]
            v₁[j,k] = η₁[j]*v₁[j,k] + ρ⁻¹ᵣ[k]*( σ₁₃[j,k] - σ₁₃[j-1,k] ) + BodyForce
            v₂[j,k] = η₂[k]*v₂[j,k] + ρ⁻¹ᵣ[k]*( σ₂₃[j,k] - σ₂₃[j,k-1] ) + BodyForce
            v[j,k] = v₁[j,k] + v₂[j,k]  #   artificial splitting of v for PML
            u[j,k] = u[j,k] + 0.5*v[j,k]*dt #   time-integration of v w/ trapezoidal rule (+half of the current v)
            σ₁₃[j-1,k] = η₁₃[j-1]*σ₁₃[j-1,k] + μᵣ[k]*( v[j,k] - v[j-1,k] )
            σ₂₃[j,k-1] = η₂₃[k-1]*σ₂₃[j,k-1] + μᵣˢᵍ[k-1]*( v[j,k] - v[j,k-1] )
        end
        if i % 10 == 0
            open("dat/U$(i).dat","w") do uo
                write(uo,u)
            end
            open("dat/V$(i).dat","w") do vo
                write(vo,v)
            end
        end
        @printf("%5d / %5d\r", i, Nₜ);
    end
end

function f(t,T)
    tⁿ = t/T
    force = tⁿ*exp(-tⁿ)/T
    return force
end

function g(x₁,x₂,x₁⁰,x₂⁰,R)
    x₁ⁿ = (x₁-x₁⁰)/R
    x₂ⁿ = (x₂-x₂⁰)/R
    distribution = exp(-x₁ⁿ^2-x₂ⁿ^2)/(pi*R^2)
    return distribution
end

function Staggered(x)
    N = length(x)
    xˢᵍ = zeros(N-1)
    for i=1:N-1
        xˢᵍ[i] = 0.5*(x[i+1] + x[i])
    end
    return xˢᵍ
end

function Bimaterial(c⁺,c⁻,d,x)
    R = 2*dx
    N = length(x)
    C = zeros(N)
    for i = 1:N
        ξ = (x[i]-d)/R
        C[i] = c⁻ + (c⁺-c⁻)*(1 + erf(ξ))/2
    end
    return C
end

function Viscosity(x₁,x₂)
    (x₁ˢᵍ,x₂ˢᵍ) = (Staggered(x₁),Staggered(x₂))
    (η₁,η₂,η₁₃,η₂₃) = (ones(N₁),ones(N₂),ones(N₁-1),ones(N₂-1))
    R = Nₚₘₗ*dx
    for i = 1:Nₚₘₗ
        p = pi/2.0/R
        η₁[i+1] = 1-η*cos(p*(x₁[i+1]-x₁[1]))^2  #   left for v₁
        η₁[N₁-i] = 1-η*cos(p*(x₁[N₁-i]-x₁[N₁]))^2   #   right for v₁
        η₂[i+1] = 1-η*cos(p*(x₂[i+1]-x₂[1]))^2  #   bottom for v₂
        η₁₃[i] = 1-η*cos(p*(x₁ˢᵍ[i]-x₁[1]))^2   #   left for σ₁₃
        η₁₃[N₁-i] = 1-η*cos(p*(x₁ˢᵍ[N₁-i]-x₁[N₁]))^2    #   right for σ₁₃
        η₂₃[i] = 1-η*cos(p*(x₂ˢᵍ[i]-x₂[1]))^2   #   bottom for σ₂₃
    end
    return η₁,η₂,η₁₃,η₂₃
end

function NumericalCoefficients(ρ,μ,μˢᵍ,x₁,x₂,x₁⁰,x₂⁰,R,dx,dt)
    ρ⁻¹ᵣ = similar(ρ)
    μᵣ = similar(μ)
    μᵣˢᵍ = similar(μˢᵍ)
    r = dt/dx
    for k = 2:N₂-1
        ρ⁻¹ᵣ[k] = r/ρ[k]
        μᵣ[k] = r*μ[k]
        μᵣˢᵍ[k-1] = r*μˢᵍ[k]
    end
    s = zeros(N₁,N₂)
    for j = 2:N₁-1, k = 2:N₂-1
        s[j,k] = (0.5/ρ[k])*g(x₁[j],x₂[k],x₁⁰,x₂⁰,R)*dt   #   spatial distribution of the source
    end
    return ρ⁻¹ᵣ,μᵣ,μᵣˢᵍ,s
end

FDTD()

param = open("param.dat","w")
println(param,"Lx=",L₁)
println(param,"Ly=",L₂)
println(param,"Nx=",N₁)
println(param,"Ny=",N₂)
println(param,"Nt=",Nₜ)
println(param,"dt=",dt)
println(param,"dx=",dx)
println(param,"cp=",c⁺)
println(param,"cm=",c⁻)