# Multiblock SBP-SAT Methodology of Symmetric Matrix Form of Elastic Wave Equations on Curvilinear Grids

Authors: Cheng Sun; Zai-Lin Yang; Guan-Xi-Xi Jiang; Yong Yang
Journal: Shock and Vibration
 (2020)

---

## Abstract

A stable and accurate finite-difference discretization of first-order elastic wave equations is derived in this work. To simplify the origin and proof of the formulas, a symmetric matrix form (SMF) for elastic wave equations is presented. The curve domain is discretized using summation-by-parts (SBP) operators, and the boundary conditions are weakly enforced using the simultaneous-approximation-term (SAT) technique, which gave rise to a provably stable high-order SBP-SAT method via the energy method. In addition, SMF can be extended to wave equations of different types (SH wave and P-SV wave) and dimensions, which can simplify the boundary derivation process and improve its applicability. Application of this approximation can divide the domain into a multiblock context for calculation, and the interface boundary conditions of blocks can also be used to simulate cracks and other structures. Several numerical simulation examples, including actual elevation within the area of Lushan, China, are presented, which verifies the viability of the framework present in this paper. The applicability of simulating elastic wave propagation and the application potential in the seismic numerical simulation of this method are also revealed.


---

## Body

## 1. Introduction

Partial differential equations (PDEs) are widely used in mathematics, physics, and other fields. However, in practical engineering, it is difficult to obtain the corresponding analytical results. Hence, a number of numerical methods have been developed for computational modeling of elastic wave propagation. For example, finite-difference methods (FDMs) [1–3], finite element methods (FEM) [4, 5], spectral element methods (SEM) [6], boundary element methods (BEM) [7–9], finite volume methods (FVM) [10], and pseudospectral methods [11, 12]. However, each method has its own merits and drawbacks.FDMs are commonly used and extremely efficient for wave propagation in seismic engineering applications due to their efficiency and ability to simulate complex layered media, among others. However, it is difficult to manage complex geometries, particularly when using a staggered grid method [13]. An additional problem of this method is that it is unsatisfactory when simulating cracked media. Generally, when the traditional FDM is used, the method to simulate the crack is to fill the crack with a low velocity medium such as water. Meanwhile, grids close to a low velocity medium must be refined for stability, which requires a high volume of calculations and difficult automatic modeling. Additionally, selected FDM methods present difficulties in terms of managing long-time stability, both in terms of spatial- and time-related aspects. These drawbacks give rise to significant limitations for the simulation of large-scale, long-time, and undulating media. Therefore, a numerical method with time and spatial stability is proposed in this work.In the context of numerical simulation, accuracy will have a significant effect on results. Compared to first- and second-order methods, higher-order FDM (HOFDM) can effectively reduce the degree of freedom required for simulating wave propagation [14]. However, HOFDMs also give rise to additional limitations and difficulties, such as discrete approximation of boundary conditions, and the establishment and simulation of complex structures such as cracks. Furthermore, a larger amount of computation and longer computation time are generally required, leading to computation limitations.In large-scale simulation cases particularly, the computational load and memory utilization of traditional methods will face significant challenges. Multiblock simulation methods proposed for composite domains can save processing memory by partitioning the region represented in a large-scale model, thereby greatly improving computing efficiency by enabling parallel calculations between blocks.A relatively appropriate HOFDM that combines initial boundary value conditions involves summation-by-parts (SBP) operators [15–19]. This approach is employed for approximating the derivatives in a spatial domain, together with the weak enforcement of boundaries and interface conditions, using either the simultaneous-approximation-term (SAT) method [20, 21] or the projection method [22–24]. The SBP-SAT method can develop strictly stable proof using the energy method, leading to robustness in spatial and time domains. Some review papers of SBP operators can be found in [25, 26] and examples of SBP-SAT can be found in [27–30].The essential characteristics of SBP operators are that they mimic the integration-by-parts (IBP) property so that it is possible to mimic the energy dissipation of a continuous problem. The basic theory of first-order SBP operators was proposed by Kreiss and Scherer [19] in 1974, whereby the energy method is applied for stability proof, leading to the HOFDM known as strict stability [31, 32].The SBP operators introduced above were derived in a uniform grid with central-difference stencils in the interior of the area and using one-sided difference stencils near the boundaries, thereby maintaining the properties of SBP operators. These operators can be referred to as classical SBP operators. Other types of SBP operators are generalized SBP [25, 33–35], multidimensional SBP [36, 37], upwind SBP [38, 39], and staggered and upwind SBP operators [40]. Generalized SBP operators have one or more characteristics of nonrepeating interior operators, nonuniform nodal distribution, and exclude one or both boundary nodes. Multidimensional SBP operators remove the limitations of classical SBP operators for applying the tensor-product formulation to multidimensional conditions because most classical SBP operators are one-dimensional. Upwind SBP operators, combining flux-splitting techniques for hyperbolic systems were proposed via noncentral-difference stencil in the interior, which introduced artificial dissipation on a nonstaggered or staggered grid.In seismology, and in a number of fields such as acoustics, electromagnetics, and oceanography, a commonly encountered difficulty is the treatment of unbounded domains. The local high-order absorbing boundary condition (ABC) [41–43] and perfectly matched layer (PML) [44–47] are the most commonly used techniques for trimming these domains, and “absorbing” outgoing waves as much as possible, in order to simulate elastic wave motion in finite computation domains. ABC is a boundary condition defined on an artificial boundary in order to make small, spurious reflections occur when elastic waves propagate at the boundary. Absorbing layers were proposed by modifying the underlying equations, in order to rapidly dampen the wave in a layer. Due to the difficulties associated with the transformation between time and frequency domains, energy estimation was not applicable to the PML within the time domain [48]. A local high-order ABC, as well as a nonreflecting boundary condition (BC), has been widely used for their feasibility, intuitive nature, and generalizability when combined with the SAT method, enabling well-proven approximates [49–58]. We employed the nonreflecting BC since it provided a simple formulation of the elastic wave problem. As an extension, a method combining PML with the SBP-SAT (summation by parts + simultaneous approximation terms) methodology of the symmetric matrix form (SMF) of elastic wave equations was developed in a general curvilinear geometry and is presented in this work. Nonreflecting BCs derived in the Laplace space can be observed in [59]. Additionally, a well-proposed SBP-SAT methodology was used on grids with nonconforming interfaces [60–64].The aim of this work was to establish a generic representation of the elastic wave equation using the SMF. As a result, the SMF framework was able to formulate SH and P-SV wave equations in a similar manner, thereby simplifying this calculation, while also facilitating the transformation of the extension of dimensions (two-dimensional [2D] and three-dimensional [3D]). Meanwhile, using the SMF, a general continuous stability analysis is presented, combined with BCs in a general hyperbolic system. An equivalent form can be derived in a curve domain, similar to the form in a cuboid domain, by using curvilinear transformation of the elastic wave equations. Thereafter, a generic stable discretization of the SMF on curvilinear coordinates was developed using SBP-SAT, where it was observed that the SMF is able to provide an equivalent type of the cuboid domain. In this instance, the results of similar stability analyses greatly simplified the calculation and proof process. Hence, this framework is applicable to a wide range of elastic wave equations.The remainder of this paper is presented as follows: in Section2, the SMF of elastic wave equations is established. The well-posed BCs of the SMF of elastic wave equations are proposed using the flux-splitting form provided in Section 3. Section 4 presents the stability analysis using the continuous energy method. A semidiscrete and applied SBP-SAT methodology is presented, and stability analysis within a discrete norm is discussed in Section 5. Numerical studies are performed in Section 6. Subsequently, conclusions are drawn and future research is presented in Section 7.


## 2. SMF of Elastic Wave Equations

### 2.1. Formulation of SH Wave

Let the Cartesian coordinates of a two-dimensional spatial domain defined byx,y and t denotes the time variable. Consider the first-order elastic wave equation of the SH form:(1)ρ∂vz∂t=∂σxz∂x+∂σyz∂y,∂σxz∂t=G∂vz∂x,∂σyz∂t=G∂vz∂x,where vz is the velocity in the z direction of the SH wave, defined by vs=G/ρ, and σxzandσyz are the stresses; ρ is the density of the elastic medium; and G is the shear modulus. Equation (1) can then be revised as a symmetric system that has a form(2)∂ω∂t=Ax∂ω∂x+Ay∂ω∂y,where(3)Ax=0vs0vs00000,Ay=00vs000vs00,ω=ρ2vz12Gσxz12GσyzT.


### 2.2. Formulation of P-SV Wave

The P-SV from the elastic wave equation can be written as(4)∂σxx∂t∂σyy∂t∂σxy∂t=C∂vx∂x∂vy∂y∂vx∂y+∂vy∂x,ρ∂vx∂t=∂σxx∂x+∂σxy∂y,ρ∂vy∂t=∂σxy∂x+∂σyy∂y,where vxandvy are the particle velocities, σij are the stresses, ρ is density, and C=CT is the stiffness matrix. In the context of the isotropic case, the stiffness matrix is described by two independent elastic coefficients, the Lame parameters λ and μ. Note that C is a symmetric positive definite:(5)C=λ+2μλ0λλ+2μ000μ.It will be convenient to use wave velocitiesvp and vs as medium properties:(6)vp=λ+2μρ,vs=μρ.Equation (4) can be formulated as a symmetric system as(7)A^0∂ω^∂t=A^x∂ω^∂x+A^y∂ω^∂y,where ω^=σxx/ρσyy/ρσxy/ρvxvyT,(8)A^0=ab000ba00000c000001000001,A^x=0001000000000011000000100,A^y=0000000001000100010001000,(9)a=vp24vp2vp2−vs2,b=−vp2−2vs24vp2vp2−vs2,c=1vs2.The matrices in (8) are symmetric, and A^0 is positive definite. Equation (7) also can be written as a symmetric first-order hyperbolic system, as a generic system (2):(10)∂ω∂t=Ax∂ω∂x+Ay∂ω∂y.In order to transform (7) to (10), we shall ansatz A^0=PTΛP, where(11)P=1/21/20001/2−1/2000001000001000001,Λ=1/2vp2−vs21/2vs21/vs211.P is an orthonormal matrix that satisfies PTP=I, denoting(12)M=PTΛ−1/2P,M−1=PTΛ1/2P.ThusMTA^0=M−1, and premultiplying M from Equation (10) yields(13)MTA^0ω^t=MTA^xω^x+MTA^yω^y=MTA^xMM−1ω^x+MTA^yMM−1ω^y,as(14)M−1ω^t=MTA^xMM−1ω^x+MTA^yMM−1ω^y.Here, we transform (14) into(15)ωt=Axωx+Ayωy,ω=M−1ω^.MatrixM satisfying the above conditions is(16)M=PTΛ−1/2P=vp2−vs2+vs2vp2−vs2−vs2000vp2−vs2−vs2vp2−vs2+vs200000vs000001000001.


### 2.3. Curvilinear Coordinates and Coordinate Transformation

Regular grids have difficulty in managing numerical simulations for geometrically complex models. Consider the elastic wave equation in a domain with a curved boundary shape, where a transform from the curvilinear domainΩ∈x,y to the cuboid domain Ω′∈ξ,η was effected (Figure 1). Here, we introduced mapping defined byFigure 1
Coordinate transformation.(17)
x
ξ
,
η
,
y
ξ
,
η
⟷
ξ
x
,
y
,
η
x
,
y
.The Jacobian determinant is(18)J=∂x∂ξ∂x∂η∂y∂ξ∂y∂η=xξyη−xηyξ.Therefore, the derivatives in each direction can be written as(19)Jωx=yηωξ−yξωη,Jωy=−xηωξ+xξωη.Using the relations, it follows that(20)Aξ=Axyη−AyxηJ,Aη=−Axyξ+AyxξJ.Thus, we obtain the same formula as (2) and (10):(21)ωt=Aξωξ+Aηωη,where AξandAη are symmetric matrices since Ax and Ay are symmetric as well.


## 2.1. Formulation of SH Wave

Let the Cartesian coordinates of a two-dimensional spatial domain defined byx,y and t denotes the time variable. Consider the first-order elastic wave equation of the SH form:(1)ρ∂vz∂t=∂σxz∂x+∂σyz∂y,∂σxz∂t=G∂vz∂x,∂σyz∂t=G∂vz∂x,where vz is the velocity in the z direction of the SH wave, defined by vs=G/ρ, and σxzandσyz are the stresses; ρ is the density of the elastic medium; and G is the shear modulus. Equation (1) can then be revised as a symmetric system that has a form(2)∂ω∂t=Ax∂ω∂x+Ay∂ω∂y,where(3)Ax=0vs0vs00000,Ay=00vs000vs00,ω=ρ2vz12Gσxz12GσyzT.


## 2.2. Formulation of P-SV Wave

The P-SV from the elastic wave equation can be written as(4)∂σxx∂t∂σyy∂t∂σxy∂t=C∂vx∂x∂vy∂y∂vx∂y+∂vy∂x,ρ∂vx∂t=∂σxx∂x+∂σxy∂y,ρ∂vy∂t=∂σxy∂x+∂σyy∂y,where vxandvy are the particle velocities, σij are the stresses, ρ is density, and C=CT is the stiffness matrix. In the context of the isotropic case, the stiffness matrix is described by two independent elastic coefficients, the Lame parameters λ and μ. Note that C is a symmetric positive definite:(5)C=λ+2μλ0λλ+2μ000μ.It will be convenient to use wave velocitiesvp and vs as medium properties:(6)vp=λ+2μρ,vs=μρ.Equation (4) can be formulated as a symmetric system as(7)A^0∂ω^∂t=A^x∂ω^∂x+A^y∂ω^∂y,where ω^=σxx/ρσyy/ρσxy/ρvxvyT,(8)A^0=ab000ba00000c000001000001,A^x=0001000000000011000000100,A^y=0000000001000100010001000,(9)a=vp24vp2vp2−vs2,b=−vp2−2vs24vp2vp2−vs2,c=1vs2.The matrices in (8) are symmetric, and A^0 is positive definite. Equation (7) also can be written as a symmetric first-order hyperbolic system, as a generic system (2):(10)∂ω∂t=Ax∂ω∂x+Ay∂ω∂y.In order to transform (7) to (10), we shall ansatz A^0=PTΛP, where(11)P=1/21/20001/2−1/2000001000001000001,Λ=1/2vp2−vs21/2vs21/vs211.P is an orthonormal matrix that satisfies PTP=I, denoting(12)M=PTΛ−1/2P,M−1=PTΛ1/2P.ThusMTA^0=M−1, and premultiplying M from Equation (10) yields(13)MTA^0ω^t=MTA^xω^x+MTA^yω^y=MTA^xMM−1ω^x+MTA^yMM−1ω^y,as(14)M−1ω^t=MTA^xMM−1ω^x+MTA^yMM−1ω^y.Here, we transform (14) into(15)ωt=Axωx+Ayωy,ω=M−1ω^.MatrixM satisfying the above conditions is(16)M=PTΛ−1/2P=vp2−vs2+vs2vp2−vs2−vs2000vp2−vs2−vs2vp2−vs2+vs200000vs000001000001.


## 2.3. Curvilinear Coordinates and Coordinate Transformation

Regular grids have difficulty in managing numerical simulations for geometrically complex models. Consider the elastic wave equation in a domain with a curved boundary shape, where a transform from the curvilinear domainΩ∈x,y to the cuboid domain Ω′∈ξ,η was effected (Figure 1). Here, we introduced mapping defined byFigure 1
Coordinate transformation.(17)
x
ξ
,
η
,
y
ξ
,
η
⟷
ξ
x
,
y
,
η
x
,
y
.The Jacobian determinant is(18)J=∂x∂ξ∂x∂η∂y∂ξ∂y∂η=xξyη−xηyξ.Therefore, the derivatives in each direction can be written as(19)Jωx=yηωξ−yξωη,Jωy=−xηωξ+xξωη.Using the relations, it follows that(20)Aξ=Axyη−AyxηJ,Aη=−Axyξ+AyxξJ.Thus, we obtain the same formula as (2) and (10):(21)ωt=Aξωξ+Aηωη,where AξandAη are symmetric matrices since Ax and Ay are symmetric as well.


## 3. Well-Posed Boundary Conditions

A generic representation of the elastic wave equation, referred to as the SMF, is presented. Using this framework, a specific class of the elastic wave equation can be analyzed and approximated. Meanwhile, the implementation of BCs is also simplified. Wave propagation is often extremely large, even in unbounded spatial domains. Considerable physical domains must be replaced by suitable-sized computational domains, by introducing artificial BCs to trim these domains. Moreover, a free-surface BC is necessary for simulating elastic wave propagation at the surface or crack structure. The traditional approach for simulating different BCs is to employ the properties of waves at different boundaries, e.g., the traction perpendicular to the surface of the free boundary is zero. Different BCs have different mathematical formulas, which gives rise to difficulties in the derivation and unity of the elastic wave equations framework. The SMF is not simply a generic representation of the elastic wave equation but also creates a generic system of different BCs. To derive a nonreflecting BC and the free-surface BC, a well-proven characteristic type of ABC was considered, where the number of variables propagated into the domain was equal to the number of variables propagated out the domain [65, 66]. Here we denote four boundaries in a domain: west, east, south, and north (Figure 1). The free-surface boundary condition at the north boundary was considered. A nonreflecting BC was proposed at the west, east, and south boundaries. More details of characteristic BCs are provided in [67, 68]. In this paper, we first introduce the splitting of the coefficient matrix:(22)Aξ±=TAξTΔAξ±ΔAξ2TAξT.Then, having(23)Aξ±=TAξ±ΔAξ±TAξ±T,Aη±=TAη±ΔAη±TAη±T.Using (23), the characteristic variables associated with propagating in and out of the domain at the boundary are(24)pwest=TAξwest−Tωwest,qwest=TAξwest+Tωwestpeast=TAξeast+Tωeast,qeast=TAξeast−Tωeastpsouth=TAηsouth−Tωsouth,qsouth=TAηsouth+Tωsouthpnorth=TAηnorth+Tωnorth,qnorth=TAηnorth−Tωnorth.Well-posed BCs posit that only incoming variables of the correct quantity must be prescribed on the boundary. Thus, the well-posed problem can be specified by the incoming characteristic variables composed ofT, which is(25)p=Rq,where R serves as a reflection coefficient, i.e.,R∈−1,1. Aimed at different BCs, R = −1 represents a free-surface BC and a nonreflecting BC, as defined by R = 0.


## 4. Continuous Energy Analysis

In this section, an energy analysis is presented for the SMF of the elastic wave equation (10), alongside the BC in (25). The energy method was based on constructing a norm for the given problem that did not grow over time. To simplify notation in the forthcoming energy analysis, a number of definitions are necessary. The result at the grid point was stacked as a vector, where u1T=u11,u12,…,u1k and u2T=u21,u22,…,u2k. The inner product was given by u1,u2A=∫xlxru1TAxu2dx, and Ax=ATx>0 in the domain u1,u2∈Lxl,xr. Thus, the corresponding norm was expressed as uA2=u,uA. Here, we again included the SMF of elastic wave equations in a bounded domain for the sake of convenience:(26)ωt=Aξωξ+Aηωη.Multiplying (26) by ωT, setting data to zero and adding its transpose, then integrating over the domain yield(27)∫ωtTω+ωTωtdS=∫ωTAξωξ+ωTAηωη+ωξTAξω+ωηTAηωdS.Equation (27) is equivalent to(28)∂∂tω2=BTwest+BTeast+BTsouth+BTnorth,where ∂/∂tω2 is the time derivative, and (28) gives(29)BTwest=−∫westωTAξωdl,BTeast=∫eastωTAξωdl,BTsouth=−∫southωTAηωdl,BTnorth=∫westωTAηωdl.A bounded energy for (28) required ∂/∂tω2≤0, as well as BTwest,BTeast,BTsouth,BTnorth≤0. In this instance, we show that the BC terms (29) were nonpositive definite. As BCs were the same except at the north boundary, we only included the treatment of the south and north boundaries:(30)BTsouth=∫southωTTAη+ΔAη+TAη+T+TAη−ΔAη−TAη−Tωdl=−∫southωTTAη+ΔAη+TAη+Tωdl≤0,(31)BTnorth=∫westωTRTAη−ΔAη+RTAη−T+TAη−ΔAη−TAη−Tωdl=ωTTAη−RΔAη+RT+ΔAη−TAη−Tωdl≤0.Equations (30) and (31) were nonpositive so that (28) satisfied the energy estimate. Using the SMF, the energy analysis could be simplified to a large extent as a generic system.


## 5. Semidiscrete Approximations

In this section, an energy-stable, semidiscrete approximation was derived. We used SBP finite-difference operators to approximate the spatial domain and imposed BCs using the SAT method, on an unstaggered grid.


### 5.1. Spatial Discrete Operators

To begin with, consider the uniform discretization of the domain,x∈xl,xr, with N + 1 grid points and step h > 0:(32)xi=ih,i=0,1,...,N+1,h=xr−xlN.The solution is stacked as a vectorU=uxL,uΔx,…,uxR, and its discrete form is V=v0,v1,…,vN+1. An SBP operator approximating ∂/∂x has the form(33)D=H−1Q,DV=Vx.The matrixD satisfies integration by parts if the following properties hold:(34)H=HT>0,Q+QT=EN−E0,where E0=e→0e→0T,EN=e→Ne→NT,e⟶0=10⋯0T,ande→N=00…1T. In this paper, 6th order SBP operators are used, as mentioned in [15].


### 5.2. SBP-SAT Methodology of SMF of Elastic Wave Equations

In this section, an SBP-SAT methodology of SMF of elastic wave equations is proposed and its numerical stability is proved. A computation domain was discretized usingNξ+1×Nη+1 grids. The numerical solution of SH wave equations on the grids is stacked as vectors:(35)ω=ω001,ω002,ω003,…,ωNξNη1,ωNξNη2,ωNξNη3T.For the P-SV wave equations, (35) then take the form(36)ω=ω001,ω002,ω003,ω004,ω005,…,ωNξNη1,ωNξNη2,ωNξNη3,ωNξNη4,ωNξNη5T.The SBP method in higher space dimensions is derived using the Kronecker product:(37)B⊗C=b0,0C⋯b0,nC⋮⋱⋮bm,0C⋯bm,nC,where B is an m+1×n+1 matrix and C is a p+1×q+1 matrix so that (37) is a m+1p+1×n+1q+1 matrix. (37) fulfills the rules of A⊗BC⊗D=AC⊗BD and A⊗BT=A⊗BT. Define I be the 3×3 identity matrix for the SH wave equation and 5×5 identity matrix for the P-SV wave equation, Iξ be the Nξ+1×Nξ+1, and Iη be the Nη+1×Nη+1 matrix; SBP operators in 2D are(38)Dξ=Dξ1D⊗Hη⊗I,Dη=Hξ⊗Dη1D⊗I.The coefficient matrices are determined in 2D yields:(39)Aξ=Iξ⊗Iη⊗Aξ1D,Aη=Iξ⊗Iη⊗Aη1D.Moreover, the inner product is defined as follows:(40)Hξη=Hξ⊗Hη⊗I.To combine with boundary conditions, we set(41)Ewest=elξ⊗Iη⊗I,Eeast=erξ⊗Iη⊗I,Esouth=Iξ⊗elη⊗I,Enorth=Iξ⊗erη⊗I.Here, the SBP-SAT discretization of (21) with the BCs (24) and (25) can be written as(42)∂ω∂t=AξDξω+AηDηω+SAT,where(43)SATwest=+Hξ−1EwestTAξwest−ΔAξwest−TAξwest−Tωwest,SATeast=−Hξ−1EeastTAξeast+ΔAξeast+TAξeast+Tωeast,SATsouth=+Hη−1EsouthTAηsouth−ΔAηsouth−TAηsouth−Tωsouth,SATnorth=−Hη−1EnorthTAηnorth+ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.For stability analysis, multiplying (42) by ωTHξη and adding the transpose lead to(44)∂∂tωξη2=BTwest+BTeast+BTsouth+BTnorth,where the boundary conditions are(45)BTwest=−ωwestTHηwestAξwestωwest+2ωwestTHηwestTAξwest−ΔAξ−TAξwest−Tωwest,BTeast=ωeastTHηeastAξeastωeast−2ωeastTHηeastTAξeast+ΔAξ+TAξeast+Tωeast,BTsouth=−ωsouthTHξsouthAηsouthωsouth+2ωsouthTHξsouthTAηsouth−ΔAη−TAηsouth−Tωsouth,BTnorth=ωnorthTHξnorthAηnorthωnorth−2ωnorthTHξnorthTAηnorth+ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.Similar to (30) and (31), we have(46)BTeast=−ωeastTHηeastTAξeast+ΔAξ+TAξeast+Tωeast+ωeastTHηeastTAξeast−ΔAξ−TAξeast−Tωeast,BTwest=−ωwestTHηwestTAξwest+ΔAξ+TAξwest+Tωwest+ωwestTHηwestTAξwest−ΔAξ−TAξwest−Tωwest,BTsouth=−ωsouthTHηsouthTAξsouth+ΔAξsouth+TAξsouth+Tωsouth+ωsouthTHηsouthTAξsouth−ΔAξsouth−TAξsouth−Tωsouth,BTnorth=ωnorthTHηnorthTAξnorth−RΔAηnorth+RT+ΔAηnorth−TAξnorth−Tωnorth−ωnorthTHξnorthTAηnorth+−RTAηnorth−ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.Since (46) are nonpositive definite matrices, ∂/∂tωξη2≤0 verifies the robustness.For large scale, as well as complex structural spatial domains, a computation block is not a preferable choice for numerical simulation. We divided the computation space into several blocks, with grids conforming at the block interface when interfaces were continuous, which had a form matching that in (47). Meanwhile, the multiblock SBP-SAT methodology of elastic wave equations was able to manage discontinuous structures such as cracks by using the interface definition. The form of two blocks is shown in Figure 2.(47)
SAT
east
1
=
−
H
ξ
−
1
1
E
east
1
T
A
ξ
east
+
1
Δ
A
ξ
east
+
1
T
A
ξ
east
+
1
T
ω
east
1
−
ω
west
2
,
SAT
west
2
=
+
H
ξ
−
1
2
E
west
2
T
A
ξ
west
−
2
Δ
A
ξ
west
−
2
T
A
ξ
west
−
2
T
ω
west
2
−
ω
east
1
,when the interface conditions were applied for simulating cracks that hadFigure 2
Multiblock model.(48)
SAT
east
1
=
−
H
ξ
−
1
E
east
T
A
ξ
east
+
Δ
A
ξ
east
+
T
A
ξ
east
+
T
ω
east
,
SAT
west
2
=
+
H
ξ
−
1
E
west
T
A
ξ
west
−
Δ
A
ξ
west
−
T
A
ξ
west
−
T
ω
west
.


## 5.1. Spatial Discrete Operators

To begin with, consider the uniform discretization of the domain,x∈xl,xr, with N + 1 grid points and step h > 0:(32)xi=ih,i=0,1,...,N+1,h=xr−xlN.The solution is stacked as a vectorU=uxL,uΔx,…,uxR, and its discrete form is V=v0,v1,…,vN+1. An SBP operator approximating ∂/∂x has the form(33)D=H−1Q,DV=Vx.The matrixD satisfies integration by parts if the following properties hold:(34)H=HT>0,Q+QT=EN−E0,where E0=e→0e→0T,EN=e→Ne→NT,e⟶0=10⋯0T,ande→N=00…1T. In this paper, 6th order SBP operators are used, as mentioned in [15].


## 5.2. SBP-SAT Methodology of SMF of Elastic Wave Equations

In this section, an SBP-SAT methodology of SMF of elastic wave equations is proposed and its numerical stability is proved. A computation domain was discretized usingNξ+1×Nη+1 grids. The numerical solution of SH wave equations on the grids is stacked as vectors:(35)ω=ω001,ω002,ω003,…,ωNξNη1,ωNξNη2,ωNξNη3T.For the P-SV wave equations, (35) then take the form(36)ω=ω001,ω002,ω003,ω004,ω005,…,ωNξNη1,ωNξNη2,ωNξNη3,ωNξNη4,ωNξNη5T.The SBP method in higher space dimensions is derived using the Kronecker product:(37)B⊗C=b0,0C⋯b0,nC⋮⋱⋮bm,0C⋯bm,nC,where B is an m+1×n+1 matrix and C is a p+1×q+1 matrix so that (37) is a m+1p+1×n+1q+1 matrix. (37) fulfills the rules of A⊗BC⊗D=AC⊗BD and A⊗BT=A⊗BT. Define I be the 3×3 identity matrix for the SH wave equation and 5×5 identity matrix for the P-SV wave equation, Iξ be the Nξ+1×Nξ+1, and Iη be the Nη+1×Nη+1 matrix; SBP operators in 2D are(38)Dξ=Dξ1D⊗Hη⊗I,Dη=Hξ⊗Dη1D⊗I.The coefficient matrices are determined in 2D yields:(39)Aξ=Iξ⊗Iη⊗Aξ1D,Aη=Iξ⊗Iη⊗Aη1D.Moreover, the inner product is defined as follows:(40)Hξη=Hξ⊗Hη⊗I.To combine with boundary conditions, we set(41)Ewest=elξ⊗Iη⊗I,Eeast=erξ⊗Iη⊗I,Esouth=Iξ⊗elη⊗I,Enorth=Iξ⊗erη⊗I.Here, the SBP-SAT discretization of (21) with the BCs (24) and (25) can be written as(42)∂ω∂t=AξDξω+AηDηω+SAT,where(43)SATwest=+Hξ−1EwestTAξwest−ΔAξwest−TAξwest−Tωwest,SATeast=−Hξ−1EeastTAξeast+ΔAξeast+TAξeast+Tωeast,SATsouth=+Hη−1EsouthTAηsouth−ΔAηsouth−TAηsouth−Tωsouth,SATnorth=−Hη−1EnorthTAηnorth+ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.For stability analysis, multiplying (42) by ωTHξη and adding the transpose lead to(44)∂∂tωξη2=BTwest+BTeast+BTsouth+BTnorth,where the boundary conditions are(45)BTwest=−ωwestTHηwestAξwestωwest+2ωwestTHηwestTAξwest−ΔAξ−TAξwest−Tωwest,BTeast=ωeastTHηeastAξeastωeast−2ωeastTHηeastTAξeast+ΔAξ+TAξeast+Tωeast,BTsouth=−ωsouthTHξsouthAηsouthωsouth+2ωsouthTHξsouthTAηsouth−ΔAη−TAηsouth−Tωsouth,BTnorth=ωnorthTHξnorthAηnorthωnorth−2ωnorthTHξnorthTAηnorth+ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.Similar to (30) and (31), we have(46)BTeast=−ωeastTHηeastTAξeast+ΔAξ+TAξeast+Tωeast+ωeastTHηeastTAξeast−ΔAξ−TAξeast−Tωeast,BTwest=−ωwestTHηwestTAξwest+ΔAξ+TAξwest+Tωwest+ωwestTHηwestTAξwest−ΔAξ−TAξwest−Tωwest,BTsouth=−ωsouthTHηsouthTAξsouth+ΔAξsouth+TAξsouth+Tωsouth+ωsouthTHηsouthTAξsouth−ΔAξsouth−TAξsouth−Tωsouth,BTnorth=ωnorthTHηnorthTAξnorth−RΔAηnorth+RT+ΔAηnorth−TAξnorth−Tωnorth−ωnorthTHξnorthTAηnorth+−RTAηnorth−ΔAηnorth+TAηnorth+−RTAηnorth−Tωnorth.Since (46) are nonpositive definite matrices, ∂/∂tωξη2≤0 verifies the robustness.For large scale, as well as complex structural spatial domains, a computation block is not a preferable choice for numerical simulation. We divided the computation space into several blocks, with grids conforming at the block interface when interfaces were continuous, which had a form matching that in (47). Meanwhile, the multiblock SBP-SAT methodology of elastic wave equations was able to manage discontinuous structures such as cracks by using the interface definition. The form of two blocks is shown in Figure 2.(47)
SAT
east
1
=
−
H
ξ
−
1
1
E
east
1
T
A
ξ
east
+
1
Δ
A
ξ
east
+
1
T
A
ξ
east
+
1
T
ω
east
1
−
ω
west
2
,
SAT
west
2
=
+
H
ξ
−
1
2
E
west
2
T
A
ξ
west
−
2
Δ
A
ξ
west
−
2
T
A
ξ
west
−
2
T
ω
west
2
−
ω
east
1
,when the interface conditions were applied for simulating cracks that hadFigure 2
Multiblock model.(48)
SAT
east
1
=
−
H
ξ
−
1
E
east
T
A
ξ
east
+
Δ
A
ξ
east
+
T
A
ξ
east
+
T
ω
east
,
SAT
west
2
=
+
H
ξ
−
1
E
west
T
A
ξ
west
−
Δ
A
ξ
west
−
T
A
ξ
west
−
T
ω
west
.


## 6. Numerical Experiments

### 6.1. Multiblock Condition of P-SV Wave

In this section, we consider a multiblock condition that can divide the computation domain into appropriate scale blocks. Each block was able to transform data through an adjacent boundary, which had the ability to effect computations in parallel. To demonstrate the potential of this multiblock treatment, we imposed nonreflecting BCs at boundaries, and continuous interface conditions for blocks. The properties in blocks werevp = 1800 m/s, vs = 1000 m/s, and ρ=1000kg/m3. The computation domain was bounded by 0,3000m×0,3000m, which involved nine blocks with Nξ=100andNη=100, as shown in Figure 3. The time step is set to dt=10−3s, and the initial data are(49)gt=exp−x−15002+y−150028000.Figure 3
Illustrative model containing multiblock.Using the SBP-SAT methodology for the SMF of elastic wave equations, we obtained amplitude diagrams for different times. Thex- and y-direction wave fields showed different waveforms in the P-SV mode since the initial data had been applied in the x direction. P-wave propagated from block five to other blocks at time = 0.3 s (Figure 4). This indicated that the P-SV wave field motivated in block five showed little change when propagating to other blocks.Figure 4
Velocity distribution inx and z directions at time = 0.3 s. (a) x direction (b) y direction.
(a)
(b)P-wave propagated to the nonreflection boundary at time = 0.9 s (Figure5), and SV-wave propagated to the nonreflection boundary at time = 1.5 s (Figure 6). This indicated the feasibility of unbounded domain simulation.Figure 5
Velocity distribution inx and y directions at time = 0.9 s. (a) x direction (b) y direction.
(a)
(b)Figure 6
Velocity distribution inx and y directions at time = 1.5 s. (a) x direction (b) y direction.
(a)
(b)The treatment of BCs is always more complex when applying the higher-order difference method. Generally, the accuracy of BCs is lower than that for the inner medium. Therefore, the stability problem often occurs at the boundary position. This problem becomes more prominent when multiple computational domains are combined. When doing so, we divided the domain into nine subblocks and 36 boundaries, which posed a significant challenge to the stability and impact of the numerical simulation. Stability is demonstrated in Section5, as part of this numerical experiment. Meanwhile, when the elastic wave propagated to the inner boundary of the medium, there was no disturbance field; that is, the block form had little influence on the wave field, which indicated that the multiblock simulation of the wave equation was feasible. The SBP-SAT methodology applied in this work clearly had a positive promoting effect on numerical simulations, which can significantly expand the simulation scale.To investigate the connection between the multiblock and propagation of the elastic wave, we next present the amplitude of velocity in thex direction at time = 0.6 s in block two independently (Figure 7). We verified that continuous interface conditions had little impact on the overall wave field when waves propagated to the interface boundaries. Meanwhile, the values of the boundary nodes of the two adjacent computational domains coincided with one another (Figure 8). It is shown that a multiblock domain can enlarge the simulation scale and significantly improve computational efficiency under reasonable settings.Figure 7
Velocity distribution in thex direction at time = 1.5 s.Figure 8
Velocity amplitude in thex direction.


### 6.2. Crack Structure for SH Wave

When the inner computational domain was considered as a whole, the boundary types in the subcomputational region could also be diversified, which provided convenience for simulating complex structures. For example, the traditional finite-difference method is required to perform complex operations or approximations in order to realize a discontinuous medium; for example, filling intermittent parts with low-speed media. In order to visually demonstrate the propagation process of elastic waves in cracked media, SH waves were utilized since waveform transformation was not involved. The properties of the medium werevs = 100 m/s, ρ=2kg/m3, and G=104Pa. The computational domain was bounded by 0,300m×0,300m, which involved nine blocks with Nξ=100andNη=100. The time step was set to dt=10−3s, and the initial data were similar to that in (49).When an SH wave propagated to the crack structure, it generated reflected and diffracted waves, which caused continuous disturbance of the wave field in the medium. The superposition and mutual influence of the wave field had a significant impact on wave propagation characteristics. Meanwhile, the amplitude of diffracted waves decreased in an obvious manner (Figure9). This helped us to explain that in the case of real earthquakes, seismic wave attenuation factors must be considered based on the impact of cracks. Moreover, the method proposed in this paper is more convenient for realizing a crack structure that is able to avoid the amplification of errors by low-speed media filling.Figure 9
Velocity distribution in thez direction (a) at time = 0.7 s, (b) at time = 0.9 s, (c) at time = 1.2 s, (d) at time = 1.4 s, (e) at time = 1.5 s, and (f) at time = 1.8 s.
(a)
(b)
(c)
(d)
(e)
(f)


### 6.3. Curve Crack Structure in Curvilinear Domain

In this numerical experiment, we considered curve crack structure in a curvilinear domain. The interface BC was set between block five and block eight, and a crack with no contact between block two and block five was employed to show the advantages of this method in modeling and numerical simulations (Figure10). The remaining settings were the same as in the previous example.Figure 10
Illustrative model containing curve cracks.In this test, the wave field was more complex, which caused persistent disturbances by reflected and diffracted waves, but without spurious waves being generated by curve mesh conditions (Figure11). We observed long-time disturbances and a reduced maximum of amplitude, which provided references for the seismic and postseismic analysis of the structure. This test also characterized the applicability and feasibility of SBP-SAT of the SMF of elastic wave equations to simulate structures with complex shapes.Figure 11
Velocity distribution in thez direction (a) at time = 0.6 s, (b) at time = 1.0 s, (c) at time = 1.4 s, (d) at time = 1.8 s, (e) at time = 2.4 s, and (f) at time = 3.0 s.
(a)
(b)
(c)
(d)
(e)
(f)


### 6.4. Large-Scale Simulation for Complex Terrain

To demonstrate the full potential of the SBP-SAT methodology of the SMF in seismology, we next considered specific terrain in Lushan, China. Elevation data were derived from a profile of Baoxing seismic station, as shown in Figure12, for an area 798 km in length. Since multiblock treatment involving complex terrain was the main focus of the numerical simulation, we employed homogeneous and isotropic parameters where vp = 6400 m/s and vs = 3695 m/s and 19×6 blocks with 201×101 grids discretized in each.Figure 12
Velocity distribution in thez direction.A method was developed combining PML [48] with the SBP-SAT methodology of the SMF of elastic wave equations, and a general curvilinear geometry was presented by this test. We focused on elastic wave propagation in the medium and the stability of the simulation. The initial data were set similar to that in (49).The wave field at 5 s (Figure13) showed that hypocentral P-wave propagated to the free surface, while Figure 14 showed the wave field at 25 s, once elastic waves had been reflected and scattered by the irregular surface, causing the wave field to be extremely complex. The wave field at 50 s and 70 s is shown in Figures 15 and 16, respectively. As can be seen, the numerical simulation led to a converging result. Since the PML methodology was used for the curve mesh, the absorption effect at ABCs was additionally strengthened and did not undermine the robust nature of the system.Figure 13
Velocity distribution in thex direction at time = 5 s.Figure 14
Velocity distribution in thex direction at time = 25 s.Figure 15
Velocity distribution in thex direction at time = 50 s.Figure 16
Velocity distribution in thex direction at time = 70 s.Figures15 and 16 show the wave field at 50 s and 70 s. The terrain was generally mountainous on the left and relatively flat on the right so that the wave field was more complex on the left due to the scattered wave. Figures 17 and 18 show the wave field at 100 s and 175 s. A long term and large-scale simulation with a stable wave field highlights the superior nature of using a multiblock approach in order to divide the large model into appropriate subblocks, which can improve simulation efficiency. Meanwhile, the continuous interface BCs had little impact on the wave field.Figure 17
Velocity distribution in thex direction at time = 100 s.Figure 18
Velocity distribution in thex direction at time = 175 s.


## 6.1. Multiblock Condition of P-SV Wave

In this section, we consider a multiblock condition that can divide the computation domain into appropriate scale blocks. Each block was able to transform data through an adjacent boundary, which had the ability to effect computations in parallel. To demonstrate the potential of this multiblock treatment, we imposed nonreflecting BCs at boundaries, and continuous interface conditions for blocks. The properties in blocks werevp = 1800 m/s, vs = 1000 m/s, and ρ=1000kg/m3. The computation domain was bounded by 0,3000m×0,3000m, which involved nine blocks with Nξ=100andNη=100, as shown in Figure 3. The time step is set to dt=10−3s, and the initial data are(49)gt=exp−x−15002+y−150028000.Figure 3
Illustrative model containing multiblock.Using the SBP-SAT methodology for the SMF of elastic wave equations, we obtained amplitude diagrams for different times. Thex- and y-direction wave fields showed different waveforms in the P-SV mode since the initial data had been applied in the x direction. P-wave propagated from block five to other blocks at time = 0.3 s (Figure 4). This indicated that the P-SV wave field motivated in block five showed little change when propagating to other blocks.Figure 4
Velocity distribution inx and z directions at time = 0.3 s. (a) x direction (b) y direction.
(a)
(b)P-wave propagated to the nonreflection boundary at time = 0.9 s (Figure5), and SV-wave propagated to the nonreflection boundary at time = 1.5 s (Figure 6). This indicated the feasibility of unbounded domain simulation.Figure 5
Velocity distribution inx and y directions at time = 0.9 s. (a) x direction (b) y direction.
(a)
(b)Figure 6
Velocity distribution inx and y directions at time = 1.5 s. (a) x direction (b) y direction.
(a)
(b)The treatment of BCs is always more complex when applying the higher-order difference method. Generally, the accuracy of BCs is lower than that for the inner medium. Therefore, the stability problem often occurs at the boundary position. This problem becomes more prominent when multiple computational domains are combined. When doing so, we divided the domain into nine subblocks and 36 boundaries, which posed a significant challenge to the stability and impact of the numerical simulation. Stability is demonstrated in Section5, as part of this numerical experiment. Meanwhile, when the elastic wave propagated to the inner boundary of the medium, there was no disturbance field; that is, the block form had little influence on the wave field, which indicated that the multiblock simulation of the wave equation was feasible. The SBP-SAT methodology applied in this work clearly had a positive promoting effect on numerical simulations, which can significantly expand the simulation scale.To investigate the connection between the multiblock and propagation of the elastic wave, we next present the amplitude of velocity in thex direction at time = 0.6 s in block two independently (Figure 7). We verified that continuous interface conditions had little impact on the overall wave field when waves propagated to the interface boundaries. Meanwhile, the values of the boundary nodes of the two adjacent computational domains coincided with one another (Figure 8). It is shown that a multiblock domain can enlarge the simulation scale and significantly improve computational efficiency under reasonable settings.Figure 7
Velocity distribution in thex direction at time = 1.5 s.Figure 8
Velocity amplitude in thex direction.


## 6.2. Crack Structure for SH Wave

When the inner computational domain was considered as a whole, the boundary types in the subcomputational region could also be diversified, which provided convenience for simulating complex structures. For example, the traditional finite-difference method is required to perform complex operations or approximations in order to realize a discontinuous medium; for example, filling intermittent parts with low-speed media. In order to visually demonstrate the propagation process of elastic waves in cracked media, SH waves were utilized since waveform transformation was not involved. The properties of the medium werevs = 100 m/s, ρ=2kg/m3, and G=104Pa. The computational domain was bounded by 0,300m×0,300m, which involved nine blocks with Nξ=100andNη=100. The time step was set to dt=10−3s, and the initial data were similar to that in (49).When an SH wave propagated to the crack structure, it generated reflected and diffracted waves, which caused continuous disturbance of the wave field in the medium. The superposition and mutual influence of the wave field had a significant impact on wave propagation characteristics. Meanwhile, the amplitude of diffracted waves decreased in an obvious manner (Figure9). This helped us to explain that in the case of real earthquakes, seismic wave attenuation factors must be considered based on the impact of cracks. Moreover, the method proposed in this paper is more convenient for realizing a crack structure that is able to avoid the amplification of errors by low-speed media filling.Figure 9
Velocity distribution in thez direction (a) at time = 0.7 s, (b) at time = 0.9 s, (c) at time = 1.2 s, (d) at time = 1.4 s, (e) at time = 1.5 s, and (f) at time = 1.8 s.
(a)
(b)
(c)
(d)
(e)
(f)


## 6.3. Curve Crack Structure in Curvilinear Domain

In this numerical experiment, we considered curve crack structure in a curvilinear domain. The interface BC was set between block five and block eight, and a crack with no contact between block two and block five was employed to show the advantages of this method in modeling and numerical simulations (Figure10). The remaining settings were the same as in the previous example.Figure 10
Illustrative model containing curve cracks.In this test, the wave field was more complex, which caused persistent disturbances by reflected and diffracted waves, but without spurious waves being generated by curve mesh conditions (Figure11). We observed long-time disturbances and a reduced maximum of amplitude, which provided references for the seismic and postseismic analysis of the structure. This test also characterized the applicability and feasibility of SBP-SAT of the SMF of elastic wave equations to simulate structures with complex shapes.Figure 11
Velocity distribution in thez direction (a) at time = 0.6 s, (b) at time = 1.0 s, (c) at time = 1.4 s, (d) at time = 1.8 s, (e) at time = 2.4 s, and (f) at time = 3.0 s.
(a)
(b)
(c)
(d)
(e)
(f)


## 6.4. Large-Scale Simulation for Complex Terrain

To demonstrate the full potential of the SBP-SAT methodology of the SMF in seismology, we next considered specific terrain in Lushan, China. Elevation data were derived from a profile of Baoxing seismic station, as shown in Figure12, for an area 798 km in length. Since multiblock treatment involving complex terrain was the main focus of the numerical simulation, we employed homogeneous and isotropic parameters where vp = 6400 m/s and vs = 3695 m/s and 19×6 blocks with 201×101 grids discretized in each.Figure 12
Velocity distribution in thez direction.A method was developed combining PML [48] with the SBP-SAT methodology of the SMF of elastic wave equations, and a general curvilinear geometry was presented by this test. We focused on elastic wave propagation in the medium and the stability of the simulation. The initial data were set similar to that in (49).The wave field at 5 s (Figure13) showed that hypocentral P-wave propagated to the free surface, while Figure 14 showed the wave field at 25 s, once elastic waves had been reflected and scattered by the irregular surface, causing the wave field to be extremely complex. The wave field at 50 s and 70 s is shown in Figures 15 and 16, respectively. As can be seen, the numerical simulation led to a converging result. Since the PML methodology was used for the curve mesh, the absorption effect at ABCs was additionally strengthened and did not undermine the robust nature of the system.Figure 13
Velocity distribution in thex direction at time = 5 s.Figure 14
Velocity distribution in thex direction at time = 25 s.Figure 15
Velocity distribution in thex direction at time = 50 s.Figure 16
Velocity distribution in thex direction at time = 70 s.Figures15 and 16 show the wave field at 50 s and 70 s. The terrain was generally mountainous on the left and relatively flat on the right so that the wave field was more complex on the left due to the scattered wave. Figures 17 and 18 show the wave field at 100 s and 175 s. A long term and large-scale simulation with a stable wave field highlights the superior nature of using a multiblock approach in order to divide the large model into appropriate subblocks, which can improve simulation efficiency. Meanwhile, the continuous interface BCs had little impact on the wave field.Figure 17
Velocity distribution in thex direction at time = 100 s.Figure 18
Velocity distribution in thex direction at time = 175 s.


## 7. Conclusions and Future Work

The primary motivation of conducting this work was to derive a high-order accurate SBP-SAT approximation of the SMF of elastic wave equations, which can realize elastic wave propagation in a medium with complex cracks and boundaries. Based on the SMF, the stability of the algorithm can be proven more concisely through the energy method, which infers that the formula calculation and implementation of the program compilation for discretizing of the elastic wave equation will be more convenient, in addition to being able to expand the dimension of elastic waves. The SBP-SAT discretization of the SMF of elastic wave equations appears to be robust, despite the fact that an area is divided into multiple blocks, and the number of BCs is increased. The method developed in the present study was successfully applied to several simulations involving cracks and multiple blocks. It was also successfully applied to actual elevation data in Lushan, China, indicating that our method has the potential for simulating actual earthquakes.In follow-up work, we aim to build on the current model with more accurate nonreflecting BCs, without undermining the model’s stability. Meanwhile, we also hope to include other SBP operators, such as upwind or staggered and upwind operators. Numerical simulation indicates that the method has wide application prospects. It can be challenging to model, discretize, and simulate actual elevation and geological structure, and strict requirements are necessary to ensure the stability of such a method.



---

*Source: 8401537-2020-01-09.xml*