waterPMLコマンドで用いている計算式

(4) 地形・水域の表現方法と媒質境界での境界条件

(Formula used in waterPML command; (4) The representation methods for the topography and water-filled regions and the conditions on material boundaries)

waterPMLコマンドでは個々の格子セルに対して 「固体」「水」「真空」のいずれかの属性を割り当てる。 したがって地表面や水面は格子セルの境界面で近似されることになる(図1)。 以下では「固体」「水」「真空」の属性別の取り扱い、 およびそれらの境界の取り扱いについて述べる。
The waterPML command assigns one of the “solid”, “water” or “vacuum” attributes to each grid cell. Therefore the ground and water surfaces are approximated at the boundaries of grid cells (Fig. 1). Below, treatments for each attributes (“solid”, “water” and “vacuum” cells) and their boundaries.


図1. waterPMLコマンドにおける地形と水域の表現方法の模式図。
Fig. 1. Schematic illustration for the representations of the topography and water-filled regions in waterPML command.

4-1. 格子セルの中心点での物性値 (Physical properties at the centers of grid cells)

4-1-1. 固体のセル (Solid cells)

地下構造の設定ファイル (小領域分割 または 層構造) において地下の全ての地点の 密度\(\rho\)、P波速度\(V_p\)、S波速度\(V_s\)が定義される。 これらを元に物理領域内の固体の格子セル中心点での密度とラメ定数 \(\lambda=\rho(V_p^2-2V_s^2)\), \(\mu=\rho V_s^2\) を設定する。
A configuration file for the subsurface structure (either subdomain or layer) defines the density \(\rho\), P-wave velocity \(V_p\), and S-wave velocity \(V_s\) over the entire subsurface. Based on these configurations, the density and Lame constants \(\lambda=\rho(V_p^2-2V_s^2)\) and \(\mu=\rho V_s^2\) at the center of each grid cell of solid type in the physical volume are determined.

4-1-2. 水のセル (Water cells)

水域の設定ファイル において水の密度\(\rho_w\)と音速\(a_w\)が定義される。 音速はP波速度に対応し、S波は水中を伝播しない。 そこで、物理領域内の水の格子セル中心点での密度を\(\rho_w\)、 ラメ定数を\(\lambda=\rho_w a_w^2\), \(\mu=0\) と設定する。
A configuration file for water-filled regions defines the density \(\rho_w\) and sound velocity \(a_w\) of water. The sound velocity corresponds to a P-wave velocity, while the S-wave does not propagate in the water. Therefore, the density \(\rho_w\) and Lame constants \(\lambda=\rho_w a_w^2\) and \(\mu=0\) are set at the center of each grid cell of water type in the physical volume.

4-1-3. 真空のセル (Vacuum cells)

真空での計算は必要ない。 そこで、真空の格子セルにおいては密度やラメ定数を設定せず、 真空であることを表すフラグのみを設定する。
Computation is not needed in the vacuum. Therefore, densities and Lame constants are not defined in vacuum cells. Instead, a flag is defined to indicate that these cells are vacuum.

4-2. 格子セルの境界での物性値 (Physical properties on the boundaries of grid cells)

での密度とラメ定数が計算に必要になる。 周囲の格子セル中心点(真空やPML領域を除く)での密度やラメ定数を平均した値を これらの地点での密度、ラメ定数として用いる。
Densities and Lame constants are required at:
They are computed by averaging the densities and Lame constants at the centers of the surrounding cells (excluding vacuum cells and cells in the PML volume).
なお 固体のセルと水のセルの境界(水底)においては ラメ定数\(\mu=0\)とする(平均ではない)。 この点を除けば水のセルも固体のセルと同様に扱う。 この方法が水域の取り扱いとして適切であることが 岡本・竹中(2005)で示されている。
Water cells are treated in the same manner as solid cells except that the Lame constant \(\mu\) is set to zero (not an average) on the boundary between solid and water cells (water floor). Okamoto and Takenaka (2005) showed the validity of this treatment.

4-3. PML領域内での属性と物性値 (Attributes and physical properties in the PML volume)

PML領域内での格子セルの属性(固体/水/真空の区別)や密度とラメ定数の値は 最寄りの物理領域内の地点に合わせる。 The attribute of each grid cell (the distinction of solid, water, or vacuum cells) and the density and Lame constants in the PML volume are set to be equal to those at the nearest location in the physical volume.

4-4. 地表面や水面の境界条件 (Boundary conditions on the ground and water surfaces)

地表面(固体セルと真空セルの境界)や水面(水セルと真空セルの境界)においては 自由表面境界条件(トラクション0)が成り立つようにする必要がある。 waterPMLコマンドではこの条件を以下の方法で達成する。
The free surface boundary condition (zero traction) is applied on the ground surface (a boundary between solid and vacuum cells) or on the water surface (a boundary between water and vacuum cells). This boundary condition is achieved in the waterPML command by the following way.

4-4-1. 応力の非対角成分 (Off-diagonal stress components)

応力の非対角成分は4つの格子セルが接する辺の中心点で定義されており、地表面や水面が含まれる。 これらの地点において応力が初期値(0)から更新されないようにすることで 自由表面境界条件を実現する(図2)。
Off-diagonal stress components are difined at the center of each edge where four grid cells contact, which includes the ground and water surfaces. The free surface boundary condition is achieved by skipping updates for the off-diagonal stress components from their initial value of zero at these locations (Fig. 2).
プログラム内での具体的な実装方法としては、 応力の定義点のみの通し番号igt(方法3-1) と 半格子点の通し番号ig(方法2) の対応表(igt2ig)を定義する際に 固体内、水中、水底(固体セルと水セルの境界)のみに igtを割り当てるようにしている。
The program implements this approach within the definition of igt2ig, which relates a consecutive index igt for stress definition points (method 3-1) and a consecutive index ig for half-grid nodes (method 2); the values of igt are assigned only to the locations in the solid, in the water, or at the water floor (a boundary between solid and water cells).


図2. 応力の非対角成分を計算する地点 (▪)、 計算をスキップして初期値(0)のままにする地点 (▫)。
Fig. 2. Locations where the off-diagonal stress components are computed (▪) and where the computation is skipped to keep them at their initial values of zero (▫).

4-4-2. 応力の対角成分 (Diagonal stress components)

応力の対角成分は格子セルの中心点で定義されるため、 地表面や水面において応力の対角成分を直接0にすることができない。 そこで、地表面や水面を挟んだ反対側(真空内)の格子セル中心点での法線応力が 固体・水側の格子セル中心点での法線応力の\(-1\)倍になるようにする(図3)。
Because diagonal stress components are defined at the centers of grid cells, they cannot be directly set as zero on the ground and water surfaces. Therefore, the waterPML command assumes that the diagonal stress component at the center of the grid cell that is in the opposite side of the ground or water surface (i.e., the cell in the vacuum) is \(-1\) times that at the center of the solid or water cell (Fig. 3).
実際には真空セル中心点での法線応力を直接設定するのではなく、 地表面や水面における速度の法線方向成分の計算時に間接的にこの条件を与える。 2節で導出した(9)式 \[\begin{eqnarray} V_i^k(\posx_{i-1/2},t) &=& c_0^k(\posx_{i-1/2})V_i^k(\posx_{i-1/2},t_{-1}) \nonumber \\ & & +c_1^k(\posx_{i-1/2})\left[ \sum_l \left\{\tau_{ik}^l\left(\posx_{i-1/2}^{k+1/2},t_{-1/2}\right) -\tau_{ik}^l\left(\posx_{i-1/2}^{k-1/2},t_{-1/2}\right)\right\} \right. \nonumber \\ & & \left. +\delta_{k0}|\Delta\posx_k|f_i(\posx_{i-1/2},t_{-1/2}) \right] \label{eq.Vik} \end{eqnarray}\] において\(\posx_{i-1/2}\)が地表面または水面に位置する場合、 \[\begin{equation} \tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right)= -\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) \label{eq.boundary.normal} \end{equation}\] とおくことで \[\begin{equation} \tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right) -\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) =2\tau_{ii}^l\left(\posx_{i-1/2}^{i+1/2},t_{-1/2}\right) =-2\tau_{ii}^l\left(\posx_{i-1/2}^{i-1/2},t_{-1/2}\right) \label{eq.boundary.normal.apply} \end{equation}\] となり、法線応力としては固体・水側のセルの値のみを用いて (\ref{eq.Vik})式の計算が可能になる。 そして(\ref{eq.boundary.normal})式は図3の状況を表しているので 自由表面境界条件が間接的に満たされる。 このアプローチはOhminato and Chouet (1997)を踏襲したものである。
In practice, the diagonal stress components at the centers of vacuum cells are not directly given; instead, this condition is indirectly used in the computation of the normal velocity component as follows. Suppose that a position \(\posx_{i-1/2}\) is on the ground or water surface in Eq. (9) of section 2 (reproduced as Eq. \ref{eq.Vik}). Assuming that normal stress components around this location satisfy Eq. (\ref{eq.boundary.normal}), we obtain Eq. (\ref{eq.boundary.normal.apply}). Using this relation, the diagonal stress components are needed only in the solid or water side to compute Eq. (\ref{eq.Vik}), and the free surface boundary is indirectly satisfied owing to Eq. (\ref{eq.boundary.normal}). This approach was introduced by Ohminato and Chouet (1997).


図3. 地表面(や水面)で法線応力を0にするための考え方。 ▪ が計算に登場する固体セル内の法線応力、 ▫ が仮想的に考える真空中の法線応力であり、 後者が前者の\(-1\)倍であるものとして 境界面での速度の法線成分(青)を計算する。
Fig. 3. A schematic illustration for the free surface boundary condition for a normal stress component on a ground (or water) surface. An imaginary normal stress in the vacuum (▫) is assumed to be \(-1\) times that in the solid cell (▪) in the computation of the normal velocity component on the boundary (blue).