waterPMLコマンドで用いている計算式
(1) 微分方程式
(Formula used in waterPML command;
(1) Differential equations)
差分法を用いた地震波動場の計算には
変位を用いる定式化と速度・応力を用いる定式化があるが、
waterPMLコマンドではPerfectly Matched Layer(PML)との融合性を考えて
後者を用いている。
本マニュアル中の全ての計算式において縮約は用いないものとする。
Finite difference method (FDM) algorithms for the seismic wavefield
are classified into two categories:
one uses the displacement
and the other does the velocity and stress.
The waterPML command uses the velocity-stress approach,
which better accommodates the Perfectly Matched Layers (PML).
The summation convension for repeated subscripts
is not applied throughout this documentation.
この節ではPMLを含めた微分方程式を導出する。
最終的に使用する式は運動方程式(\ref{eq.motion.Vik.sumT.PML})と
構成則(\ref{eq.constitutive.dTijl.sumV.PML})である。
これらの方程式の解を
(\ref{eq.sum.Vik})(\ref{eq.sum.Tijl})
に代入すれば速度場・応力場が得られる。
弾性定数として(\ref{eq.Cijpq.use})式を利用し、
PML領域における吸収の強さを表す\(\alpha^m\)としては
(\ref{eq.alpha.use})式を利用する。
This section derives the differential equations including the PML.
The equations finally used are
the equation of motion (\ref{eq.motion.Vik.sumT.PML})
and the constitutive low (\ref{eq.constitutive.dTijl.sumV.PML}).
Inserting them into (\ref{eq.sum.Vik}) and (\ref{eq.sum.Tijl})
gives the velocity and stress fields.
Eqs. (\ref{eq.Cijpq.use}) and (\ref{eq.alpha.use})
are used for the elastic constants
and the strength of absorption in the PML volume (\(\alpha^m\)),
respectively.
1-1. 運動方程式
(Equation of motion)
変位の\(x_i\)方向成分を\(U_i\)、
応力の\((x_i,x_k)\)成分を\(\tau_{ik}\)として、
弾性体の運動方程式は
\[\begin{equation}
\rho \ddot{U_i} = \sum_k \tau_{ik,k} + f_i
\label{eq.motion.Ui}
\end{equation}\]
と書ける。
ここで\(\rho\)は密度、
\(f_i\)は震源として与える等価体積力の\(x_i\)方向成分、
「\(^\dot{ }\)」と「\(_{,k}\)」はそれぞれ時間微分と\(x_k\)方向の空間微分である。
(\ref{eq.motion.Ui})式は速度の\(x_i\)方向成分\(V_i\)を用いて
以下のように書き直せる。
\[\begin{equation}
\rho \dot{V}_i = \sum_k \tau_{ik,k} + f_i
\label{eq.motion.Vi}
\end{equation}\]
これを微分する変数ごとに分離して
\[\begin{equation}
\rho \dot{V}_i^k = \tau_{ik,k} + \delta_{k0}f_i
\label{eq.motion.Vik}
\end{equation}\]
\[\begin{equation}
V_i = \sum_k V_i^k
\label{eq.sum.Vik}
\end{equation}\]
のように書くこともできる。
ここでは\(f_i\)は\(x_0\)での微分の式に含めた。
\(\delta_{k0}\)はクロネッカーのデルタである。
The equation of motion of an elastic medium
is given by (\ref{eq.motion.Ui}),
where \(U_i\) is the \(x_i\)-component of the displacement,
\(\tau_{ik}\) is the \((x_i,x_k)\)-component of the stress,
\(\rho\) is the density of the medium,
\(f_i\) is the \(x_i\)-component of the equivalent body force
exerted at the source,
and “\(^\dot{}\)” and “\(_{,k}\)” are
time- and spatial-derivatives in \(x_k\) direction, respectively.
Eq. (\ref{eq.motion.Ui}) can be rewritten as (\ref{eq.motion.Vi}),
where \(V_i\) is the \(x_i\)-component of the velocity.
This equation can be divided into the terms
for each “\(_{,k}\)” as
(\ref{eq.motion.Vik}) and (\ref{eq.sum.Vik}).
Here \(f_i\) was included in the formula for \(x_0\)th derivative,
and \(\delta_{k0}\) is the Kronecker delta.
1-2. 構成方程式
(Constitutive law)
変位と応力の間の構成則は弾性定数\(C_{ijpq}\)を用いて
\[\begin{equation}
\tau_{ij} = \sum_{p,q} C_{ijpq}\frac{U_{p,q}+U_{q,p}}{2}
\label{eq.constitutive.Tij.epsilon}
\end{equation}\]
と書けるが、\(U_{p,q}\)の係数を\(C’_{ijpq}\)と書くことにすれば
\[\begin{equation}
\tau_{ij} = \sum_{p,q} C’_{ijpq} U_{p,q}
\label{eq.constitutive.Tij.U}
\end{equation}\]
と書き直せる。
実際に使用するのは等方弾性体の式
\[\begin{equation}
C_{ijpq}=C’_{ijpq}=\lambda\delta_{ij}\delta_{pq}
+\mu(\delta_{ip}\delta_{jq}+\delta_{iq}\delta_{jp})
\label{eq.Cijpq.use}
\end{equation}\]
(\(\lambda\), \(\mu\): ラメ定数)であるので
(\ref{eq.constitutive.Tij.U})式は
\[\begin{equation}
\tau_{ij} = \sum_{p,q} C_{ijpq} U_{p,q}
\label{eq.constitutive.Tij}
\end{equation}\]
と書ける。
(\ref{eq.constitutive.Tij})式の両辺を時間\(t\)で微分して
\[\begin{equation}
\dot{\tau}_{ij} = \sum_{p,q} C_{ijpq} V_{p,q}
\label{eq.constitutive.dTij}
\end{equation}\]
を得る。これを微分する変数ごとに分離して
\[\begin{equation}
\dot{\tau}_{ij}^l = \sum_p C_{ijpl} V_{p,l}
\label{eq.constitutive.dTijl}
\end{equation}\]
\[\begin{equation}
\tau_{ij} = \sum_l \tau_{ij}^l
\label{eq.sum.Tijl}
\end{equation}\]
のように書くこともできる。
The constitutive law between the displacement and stress
is given by Eq. (\ref{eq.constitutive.Tij.epsilon}),
where \(C_{ijpq}\) is the elastic constant.
This equation can be rewritten as (\ref{eq.constitutive.Tij.U}),
where \(C’_{ijpq}\) be the coefficient for \(U_{p,q}\).
For the elastic constant,
Eq. (\ref{eq.Cijpq.use}) is used,
where \(\lambda\) and \(\mu\) are the Lame constants.
Then Eq. (\ref{eq.constitutive.Tij.U}) can be rewritten as
(\ref{eq.constitutive.Tij}).
Differentiating this equation gives (\ref{eq.constitutive.dTij}),
which can further divided by the terms
for each “\(_{,l}\)” as
(\ref{eq.constitutive.dTijl}) and (\ref{eq.sum.Tijl}).
1-3. 方程式系
(The equation system)
(\ref{eq.motion.Vik})(\ref{eq.sum.Vik})
(\ref{eq.constitutive.dTijl})(\ref{eq.sum.Tijl})式は
1つの閉じた方程式系をなしている。
変数の数を減らすために4本の式から\(V_i\)と\(\tau_{ij}\)を消去して
\[\begin{equation}
\rho \dot{V}_i^k = \sum_l \tau_{ik,k}^l + \delta_{k0}f_i
\label{eq.motion.Vik.sumT}
\end{equation}\]
\[\begin{equation}
\dot{\tau}_{ij}^l = \sum_p \left(C_{ijpl} \sum_k V_{p,l}^k\right)
\label{eq.constitutive.dTijl.sumV}
\end{equation}\]
と書き直しておく。
これらは1つの閉じた方程式系をなしており、
初期条件・境界条件を与えて解くことができる。
Eqs. (\ref{eq.motion.Vik}), (\ref{eq.sum.Vik}),
(\ref{eq.constitutive.dTijl}), and (\ref{eq.sum.Tijl})
constitute a closed equation system.
To reduce the number of variables,
let us remove \(V_i\) and \(\tau_{ij}\) from the equations.
The results are Eqs. (\ref{eq.motion.Vik.sumT})
and (\ref{eq.constitutive.dTijl.sumV}).
They constitute a closed equation system
that can be solved with initail and boundary conditions.
1-4. PMLの導入
(Introducing the PML)
(\ref{eq.motion.Vik.sumT})(\ref{eq.constitutive.dTijl.sumV})式の
\(\PartialDiff{}{x_k}\)を
\(\frac{1}{S_k(x_k)}\PartialDiff{}{x_k}\)で
置き換ることによってPMLを実現する。
ここで
\[\begin{equation}
S_k(x_k)
\equiv 1-\frac{i\alpha^k(|x_k-x_k^b|)}{\omega}
= \frac{i\omega+\alpha^k(|x_k-x_k^b|)}{i\omega}
\label{eq.Sk}
\end{equation}\]
である。
\(\alpha^k\)はPML領域における波動の吸収の強さを表す位置の関数であり、
時刻にはよらないものとする。
\(x_k\)軸に直交する物理領域端の外側のPML領域内でのみ
\(\alpha^k\neq 0\)であるものとし、
その関数形としてはFesta and Nielsen (2003)を参照して
\[\begin{equation}
\alpha^k(\xi)=A\frac{V_p(x_k^b)}{L_{pml}}\left(\frac{\xi}{L_{pml}}\right)^n
\label{eq.alpha.use}
\end{equation}\]
を用いる。
ここで\(x_l^b\)は最寄りの物理領域端の座標、
\(L\)は最寄りの物理領域端からの距離、
\(V_p\)はP波速度、
\(L_{pml}\)はPML領域の\(x_k\)方向の厚さ、
\(A\)と\(n\)は定数である。
The PML is realized by replacing \(\PartialDiff{}{x_k}\)
in Eqs. (\ref{eq.motion.Vik.sumT}) and (\ref{eq.constitutive.dTijl.sumV})
with \(\frac{1}{S_k}\PartialDiff{}{x_k}\),
where \(S_k\) is defined by Eq. (\ref{eq.Sk});
\(\alpha^k\) is the strength of absorption in the PML volume
that is a function of position but is independent of time.
A non-zero \(\alpha^k\) is used only in the PML volumes
outside of the outer planes of the physical volume
perpendicular to the \(x_k\)-axis.
Following Festa and Nielsen (2003),
Eq. (\ref{eq.alpha.use}) is used for \(\alpha^k\),
where \(x_k^b\) is
the coordinate of the nearest outer boundary of the physical volume,
\(L\) is the distance from the boundary,
\(V_p\) is the P-wave velocity,
\(L_{PML}\) is the thickness of the PML volume in the \(x_k\) direction,
and \(A\) and \(n\) are constants.
(\ref{eq.motion.Vik.sumT})式の\(\PartialDiff{}{x_k}\)を
\(\frac{1}{S_k}\PartialDiff{}{x_k}\)で置き換えると
\[\begin{equation}
\rho \dot{V}_i^k = \frac{1}{S_k}\sum_l \tau_{ik,k}^l + \delta_{k0}f_i
= \frac{i\omega}{i\omega+\alpha^k} \sum_l \tau_{ik,k}^l + \delta_{k0}f_i
\label{eq.motion.Vik.sumT.arrange1}
\end{equation}\]
となる。
但しこれ以降、簡単のため引数の\(x_k\)などは省略する。
(\ref{eq.motion.Vik.sumT.arrange1})式を周波数領域で書くと
\[\begin{equation}
\rho i\omega V_i^k =
\frac{i\omega}{i\omega+\alpha^k} \sum_l \tau_{ik,k}^l + \delta_{k0}f_i
\label{eq.motion.Vik.sumT.arrange2}
\end{equation}\]
両辺に\((i\omega+\alpha^k)/i\omega\)を掛けて
\[\begin{equation}
\rho (i\omega+\alpha^k) V_i^k
= \sum_l \tau_{ik,k}^l + \frac{i\omega+\alpha^k}{i\omega}\delta_{k0}f_i
= \sum_l \tau_{ik,k}^l + \delta_{k0}\left(1+\frac{\alpha^k}{i\omega}\right)f_i
\label{eq.motion.Vik.sumT.arrange3}
\end{equation}\]
再び時間領域に直すと
\[\begin{equation}
\rho \left(\PartialDiff{}{t}+\alpha^k\right) V_i^k
= \sum_l \tau_{ik,k}^l + \delta_{k0}\left(1+\alpha^k \int dt \right)f_i
\label{eq.motion.Vik.sumT.PML.noAssumption}
\end{equation}\]
となる。
\(\alpha^k\)は物理領域内でゼロ、
等価体積力\(f_i\)はPML領域内でゼロであるので、
両者の積\(\alpha^kf_i\)は全領域でゼロである。
したがって(\ref{eq.motion.Vik.sumT.PML.noAssumption})式で
\(\alpha^k \int dt f_i=0\)とおくことができて
\[\begin{equation}
\rho \left(\PartialDiff{}{t}+\alpha^k\right) V_i^k
= \sum_l \tau_{ik,k}^l + \delta_{k0} f_i
\label{eq.motion.Vik.sumT.PML}
\end{equation}\]
を得る。
Replacing \(\PartialDiff{}{x_k}\) by \(\frac{1}{S_k}\PartialDiff{}{x_k}\)
in Eq. (\ref{eq.motion.Vik.sumT}) yields
(\ref{eq.motion.Vik.sumT.arrange1});
the arguments (e.g., \(x_k\)) are omitted hereafter for simplicity.
This equation is rewritten in the frequency domain as
(\ref{eq.motion.Vik.sumT.arrange2}).
Multiplying \((i\omega+\alpha^k)/i\omega\)
for the both sides of this equation gives
(\ref{eq.motion.Vik.sumT.arrange3}).
Transforming this equation to the time domain gives
Eq. (\ref{eq.motion.Vik.sumT.PML.noAssumption}).
Because \(\alpha^k\) is zero in the physical volume
and the equivalent body force \(f_i\) is zero in the PML volume,
the product \(\alpha^kf_i\) is zero throughout the entire computational volume.
Therefore, \(\alpha^k \int dt f_i=0\) holds
in Eq. (\ref{eq.motion.Vik.sumT.PML.noAssumption}),
giving Eq. (\ref{eq.motion.Vik.sumT.PML}).
同様に(\ref{eq.constitutive.dTijl.sumV})式の
\(\PartialDiff{}{x_l}\)を\(\frac{1}{S_l}\PartialDiff{}{x_l}\)で置き換えると
\[\begin{equation}
\left(\PartialDiff{}{t}+\alpha^l\right) \tau_{ij}^l
= \sum_p C_{ijpl} \left( \sum_k V_{p,l}^k\right)
\label{eq.constitutive.dTijl.sumV.PML}
\end{equation}\]
が得られる。
In the same way, replacing \(\PartialDiff{}{x_l}\)
by \(\frac{1}{S_l}\PartialDiff{}{x_l}\)
in Eq. (\ref{eq.constitutive.dTijl.sumV}) gives
(\ref{eq.constitutive.dTijl.sumV.PML}).