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を含めた微分方程式を導出する。 最終的に使用する式は運動方程式( $20$ )と構成則( $21$ )である。 これらの方程式の解を (

4

)(

11

) に代入すれば速度場・応力場が得られる。弾性定数として(

7

)式を利用し、 PML領域における吸収の強さを表す

α^{m}

としては (

15

)式を利用する。
This section derives the differential equations including the PML. The equations finally used are the equation of motion ( $20$ ) and the constitutive low ( $21$ ). Inserting them into (

4

) and (

11

) gives the velocity and stress fields. Eqs. (

7

) and (

15

) are used for the elastic constants and the strength of absorption in the PML volume (

α^{m}

), respectively.

1-1. 運動方程式 (Equation of motion)

変位の

x_{i}

方向成分を

U_{i}

、応力の

(x_{i}, x_{k})

成分を

τ_{i k}

として、弾性体の運動方程式は

\begin{matrix} (1) & ρ \ddot{U_{i}} = \sum_{k} τ_{i k, k} + f_{i} \end{matrix}

と書ける。ここで

ρ

は密度、

f_{i}

は震源として与える等価体積力の

x_{i}

方向成分、「

^{\dot{}}

」と「

_{, k}

」はそれぞれ時間微分と

x_{k}

方向の空間微分である。 (

1

)式は速度の

x_{i}

方向成分

V_{i}

を用いて以下のように書き直せる。

\begin{matrix} (2) & ρ {\dot{V}}_{i} = \sum_{k} τ_{i k, k} + f_{i} \end{matrix}

これを微分する変数ごとに分離して

\begin{matrix} (3) & ρ {\dot{V}}_{i}^{k} = τ_{i k, k} + δ_{k 0} f_{i} \end{matrix}

\begin{matrix} (4) & V_{i} = \sum_{k} V_{i}^{k} \end{matrix}

のように書くこともできる。ここでは

f_{i}

は

x_{0}

での微分の式に含めた。

δ_{k 0}

はクロネッカーのデルタである。
The equation of motion of an elastic medium is given by (

1

), where

U_{i}

is the

x_{i}

-component of the displacement,

τ_{i k}

is the

(x_{i}, x_{k})

-component of the stress,

ρ

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. (

1

) can be rewritten as (

2

), where

V_{i}

is the

x_{i}

-component of the velocity. This equation can be divided into the terms for each “

_{, k}

” as (

3

) and (

4

). Here

f_{i}

was included in the formula for

x_{0}

th derivative, and

δ_{k 0}

is the Kronecker delta.

1-2. 構成方程式 (Constitutive law)

変位と応力の間の構成則は弾性定数

C_{i j p q}

を用いて

\begin{matrix} (5) & τ_{i j} = \sum_{p, q} C_{i j p q} \frac{U_{p, q} + U_{q, p}}{2} \end{matrix}

と書けるが、

U_{p, q}

の係数を

C_{i j p q}^{'}

と書くことにすれば

\begin{matrix} (6) & τ_{i j} = \sum_{p, q} C_{i j p q}^{'} U_{p, q} \end{matrix}

と書き直せる。実際に使用するのは等方弾性体の式

\begin{matrix} (7) & C_{i j p q} = C_{i j p q}^{'} = λ δ_{i j} δ_{p q} + μ (δ_{i p} δ_{j q} + δ_{i q} δ_{j p}) \end{matrix}

(

λ

μ

: ラメ定数)であるので (

6

)式は

\begin{matrix} (8) & τ_{i j} = \sum_{p, q} C_{i j p q} U_{p, q} \end{matrix}

と書ける。 (

8

)式の両辺を時間

t

で微分して

\begin{matrix} (9) & {\dot{τ}}_{i j} = \sum_{p, q} C_{i j p q} V_{p, q} \end{matrix}

を得る。これを微分する変数ごとに分離して

\begin{matrix} (10) & {\dot{τ}}_{i j}^{l} = \sum_{p} C_{i j p l} V_{p, l} \end{matrix}

\begin{matrix} (11) & τ_{i j} = \sum_{l} τ_{i j}^{l} \end{matrix}

のように書くこともできる。
The constitutive law between the displacement and stress is given by Eq. (

5

), where

C_{i j p q}

is the elastic constant. This equation can be rewritten as (

6

), where

C_{i j p q}^{'}

be the coefficient for

U_{p, q}

. For the elastic constant, Eq. (

7

) is used, where

λ

and

μ

are the Lame constants. Then Eq. (

6

) can be rewritten as (

8

). Differentiating this equation gives (

9

), which can further divided by the terms for each “

_{, l}

” as (

10

) and (

11

1-3. 方程式系 (The equation system)

(

3

)(

4

) (

10

)(

11

)式は 1つの閉じた方程式系をなしている。変数の数を減らすために4本の式から

V_{i}

と

τ_{i j}

を消去して

\begin{matrix} (12) & ρ {\dot{V}}_{i}^{k} = \sum_{l} τ_{i k, k}^{l} + δ_{k 0} f_{i} \end{matrix}

\begin{matrix} (13) & {\dot{τ}}_{i j}^{l} = \sum_{p} (C_{i j p l} \sum_{k} V_{p, l}^{k}) \end{matrix}

と書き直しておく。これらは1つの閉じた方程式系をなしており、初期条件・境界条件を与えて解くことができる。
Eqs. (

3

), (

4

), (

10

), and (

11

) constitute a closed equation system. To reduce the number of variables, let us remove

V_{i}

and

τ_{i j}

from the equations. The results are Eqs. (

12

) and (

13

). They constitute a closed equation system that can be solved with initail and boundary conditions.

1-4. PMLの導入 (Introducing the PML)

(

12

)(

13

)式の

\frac{\partial}{\partial x_{k}}

を

\frac{1}{S_{k} (x_{k})} \frac{\partial}{\partial x_{k}}

で置き換ることによってPMLを実現する。ここで

\begin{matrix} (14) & S_{k} (x_{k}) \equiv 1 - \frac{i α^{k} (| x_{k} - x_{k}^{b} |)}{ω} = \frac{i ω + α^{k} (| x_{k} - x_{k}^{b} |)}{i ω} \end{matrix}

である。

α^{k}

はPML領域における波動の吸収の強さを表す位置の関数であり、時刻にはよらないものとする。

x_{k}

軸に直交する物理領域端の外側のPML領域内でのみ

α^{k} \neq 0

であるものとし、その関数形としてはFesta and Nielsen (2003)を参照して

\begin{matrix} (15) & α^{k} (ξ) = A \frac{V_{p} (x_{k}^{b})}{L_{p m l}} {(\frac{ξ}{L_{p m l}})}^{n} \end{matrix}

を用いる。ここで

x_{l}^{b}

は最寄りの物理領域端の座標、

L

は最寄りの物理領域端からの距離、

V_{p}

はP波速度、

L_{p m l}

はPML領域の

x_{k}

方向の厚さ、

A

と

n

は定数である。
The PML is realized by replacing

\frac{\partial}{\partial x_{k}}

in Eqs. (

12

) and (

13

) with

\frac{1}{S_{k}} \frac{\partial}{\partial x_{k}}

, where

S_{k}

is defined by Eq. (

14

);

α^{k}

is the strength of absorption in the PML volume that is a function of position but is independent of time. A non-zero

α^{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. (

15

) is used for

α^{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_{P M L}

is the thickness of the PML volume in the

x_{k}

direction, and

A

and

n

are constants.

(

12

)式の

\frac{\partial}{\partial x_{k}}

を

\frac{1}{S_{k}} \frac{\partial}{\partial x_{k}}

で置き換えると

\begin{matrix} (16) & ρ {\dot{V}}_{i}^{k} = \frac{1}{S_{k}} \sum_{l} τ_{i k, k}^{l} + δ_{k 0} f_{i} = \frac{i ω}{i ω + α^{k}} \sum_{l} τ_{i k, k}^{l} + δ_{k 0} f_{i} \end{matrix}

となる。但しこれ以降、簡単のため引数の

x_{k}

などは省略する。 (

16

)式を周波数領域で書くと

\begin{matrix} (17) & ρ i ω V_{i}^{k} = \frac{i ω}{i ω + α^{k}} \sum_{l} τ_{i k, k}^{l} + δ_{k 0} f_{i} \end{matrix}

両辺に

(i ω + α^{k}) / i ω

を掛けて

\begin{matrix} (18) & ρ (i ω + α^{k}) V_{i}^{k} = \sum_{l} τ_{i k, k}^{l} + \frac{i ω + α^{k}}{i ω} δ_{k 0} f_{i} = \sum_{l} τ_{i k, k}^{l} + δ_{k 0} (1 + \frac{α^{k}}{i ω}) f_{i} \end{matrix}

再び時間領域に直すと

\begin{matrix} (19) & ρ (\frac{\partial}{\partial t} + α^{k}) V_{i}^{k} = \sum_{l} τ_{i k, k}^{l} + δ_{k 0} (1 + α^{k} \int d t) f_{i} \end{matrix}

となる。

α^{k}

は物理領域内でゼロ、等価体積力

f_{i}

はPML領域内でゼロであるので、両者の積

α^{k} f_{i}

は全領域でゼロである。したがって(

19

)式で

α^{k} \int d t f_{i} = 0

とおくことができて

\begin{matrix} (20) & ρ (\frac{\partial}{\partial t} + α^{k}) V_{i}^{k} = \sum_{l} τ_{i k, k}^{l} + δ_{k 0} f_{i} \end{matrix}

を得る。
Replacing

\frac{\partial}{\partial x_{k}}

\frac{1}{S_{k}} \frac{\partial}{\partial x_{k}}

in Eq. (

12

) yields (

16

); the arguments (e.g.,

x_{k}

) are omitted hereafter for simplicity. This equation is rewritten in the frequency domain as (

17

). Multiplying

(i ω + α^{k}) / i ω

for the both sides of this equation gives (

18

). Transforming this equation to the time domain gives Eq. (

19

). Because

α^{k}

is zero in the physical volume and the equivalent body force

f_{i}

is zero in the PML volume, the product

α^{k} f_{i}

is zero throughout the entire computational volume. Therefore,

α^{k} \int d t f_{i} = 0

holds in Eq. (

19

), giving Eq. (

20

).

同様に(

13

)式の

\frac{\partial}{\partial x_{l}}

を

\frac{1}{S_{l}} \frac{\partial}{\partial x_{l}}

で置き換えると

\begin{matrix} (21) & (\frac{\partial}{\partial t} + α^{l}) τ_{i j}^{l} = \sum_{p} C_{i j p l} (\sum_{k} V_{p, l}^{k}) \end{matrix}

が得られる。
In the same way, replacing

\frac{\partial}{\partial x_{l}}

\frac{1}{S_{l}} \frac{\partial}{\partial x_{l}}

in Eq. (

13

) gives (

21