A GFDM with PML for seismic wave equations in heterogeneous media

Journal of Computational and Applied Mathematics 252 (2013) 40–51

Contents lists available at SciVerse ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

A GFDM with PML for seismic wave equations in heterogeneous media J.J. Benito a , F. Ureña b,∗ , L. Gavete c , E. Salete a , A. Muelas a a

Departamento de Construcción y Fabricación, Universidad Nacional de Educación a Distancia, Spain

b

Departamento de Matemática Aplicada, Universidad de Castilla-La Mancha, Spain

c

Departamento de Matemática Aplicada a los Recursos Naturales, Universidad Politécnica de Madrid, Spain

article

info

Article history: Received 15 September 2011 Received in revised form 29 July 2012 MSC: 65M06 65M12 74S20 80M20 Keywords: Meshless methods Generalized finite difference method Moving least squares Seismic waves Perfectly matched layer

abstract The interior of the Earth is heterogeneous with different materials and may have complex geometry. The free surface can also be uneven. Therefore, the use of a meshless method ‘‘with the possibility of using an irregular grid-point distribution’’ can be of interest for modelling this kind of problem. This paper shows the application of GFDM to the problem of seismic wave propagation in 2-D. To use this method in unbounded domains one must truncate the computational grid-point avoiding reflection from the edges. PML absorbing boundary condition has then been included in the numerical model proposed in this work. © 2012 Elsevier B.V. All rights reserved.

1. Introduction During recent years, meshless methods have emerged as a class of effective numerical methods which are capable of avoiding the difficulties encountered in conventional computational mesh based methods. An important path in the evolution of meshless methods has been the development of the Generalized Finite Difference Method (GFDM), also called the meshless finite difference method. The bases of the GFD were published in the early seventies. Jensen [1] was the first to introduce a fully arbitrary mesh. He considered Taylor’s series expansions interpolated on six-node stars in order to derive the finite difference (FD) formulae approximating derivatives up to the second order. Perrone and Kao [2] suggested that additional nodes in the six-point scheme should be considered and an averaging out process for the generalization of finite difference coefficients applied. The idea of using an eight node star and weighting functions, to obtain finite difference formulae for irregular meshes, was first put forward in [3] using moving least squares (MLS) interpolation, and an advanced version of the GFDM was given in [4]. Benito et al. [5] reported that the solution of the generalized finite difference method depends on the number of nodes in the cloud, the relative coordinates of the nodes with respect to the star node and on the weight function employed.

∗

Corresponding author. E-mail addresses: [email protected] (J.J. Benito), [email protected] (F. Ureña), [email protected] (L. Gavete), [email protected] (E. Salete), [email protected] (A. Muelas). 0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.08.007

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

41

An h-adaptive method in GFDM is described in [6–8]. In this paper, GFDM is applied to the solution of the problem of seismic wave propagation using regular or irregular meshes. The Finite Difference Method (FDM) is a robust numerical method applicable to structurally complex media. Due to its relative accuracy and computational efficiency it is the dominant method in modelling earthquake motion (see [9,10]). The perfectly matched layer (PML) absorbing boundary performs more efficiently and more accurately than most traditional or differential equation-based absorbing boundaries [11–14]. The paper is organized as follows. Section 1 is an introduction. Section 2 describes the GFDM obtaining the explicit generalized difference schemes for seismic wave propagation. Section 3 describes the recursive equations with PML in x-direction and y-direction. Some numerical results are included in Section 4 for homogeneous and heterogeneous domains. Finally, in Section 5 some conclusions are given. 2. Explicit generalized difference schemes for the seismic wave propagation problem 2.1. Equation of motion for a perfectly elastic, homogeneous and isotropic medium The equations of motion for a perfectly elastic, homogeneous, isotropic medium in 2-D are

 2 2 2 2 ∂ Ux (x, y, t )  2 ∂ Ux (x, y, t ) 2 ∂ Ux (x, y, t ) 2 2 ∂ Uy (x, y, t )   = α + β + (α − β )  ∂t2 ∂ x2 ∂ y2 ∂ x∂ y 2 2 2 2  ∂ Uy (x, y, t ) ∂ U y ( x, y , t ) ∂ Uy (x, y, t ) ∂ Ux (x, y, t )   = β2 + α2 + (α 2 − β 2 )  ∂t2 ∂ x2 ∂ y2 ∂ x∂ y

(1)

with the initial conditions

 Ux (x, y, 0) = f1 (x, y);  ∂ Ux (x, y, 0) = f3 (x, y); ∂t

U y ( x, y , 0 ) = f 2 ( x, y )

∂ Uy (x, y, 0) = f 4 ( x, y ) ∂t

(2)

and the boundary condition

 ∂ Ux (x0 , y0 , t )   = g1 (t ) a1 Ux (x0 , y0 , t ) + b1 ∂n  ∂ U (x , y , t )  a2 Uy (x0 , y0 , t ) + b2 y 0 0 = g2 (t ) ∂n

on Γ

(3)

where a1 , a2 , b1 , b2 are constants, and f1 (x, y), f2 (x, y), f3 (x, y), f4 (x, y), g1 (t ), g2 (t ) are known functions,

   λ + 2µ   α =  ρ    µ   β =  ρ ρ is the density, λ and µ are Lamé elastic coefficients and Γ is the boundary of Ω . 2.2. Heterogeneous formulation If density and Lame’s coefficients λ and µ are functions of spatial coordinates, in agreement with [9], in seismological problems the homogeneous approach can be complicated and, then, the so-called heterogeneous formulation is preferred. Therefore, the same formulae are used for all points in the domain and material discontinuities are accounted for spatial variation of the material parameters. 2.3. A GFDM explicit scheme The aim is to obtain explicit linear expressions for the approximation of partial derivatives in the points of the domain. First of all, an irregular grid or cloud of points is generated in the domain Ω ∪ Γ . On defining the central node with a set of nodes surrounding that node, the star then refers to a group of established nodes in relation to a central node. Every node in the domain has an associated star assigned to it.

42

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

This scheme uses the central-difference form for the time derivative

 2 unx,+01 − 2unx,0 + unx,−01 ∂ Ux (x0 , y0 , n1t )    =  ∂t2 (1t )2  uny,+01 − 2uny,0 + uny,−01  ∂ 2 Uy (x0 , y0 , n1t )   = . ∂t2 (1t )2

(4)

Following [5,6,8], the explicit difference formulae for the spatial derivatives are obtained

 2 N  ∂ Ux (x0 , y0 , n1t )  n  = − m u + mj unx,j  0 x , 0   ∂ x2  j = 1     N    ∂ 2 Uy (x0 , y0 , n1t )  n  = − m u + mj uny,j  0 y,0  2  ∂ x  j =1     N  2   ∂ Ux (x0 , y0 , n1t )  n  = −η u ηj unx,j +  0 x , 0   ∂ y2 j=1 N    ∂ 2 Uy (x0 , y0 , n1t )   = −η0 uny,0 + ηj uny,j   2  ∂ y  j=1     N  2   ∂ Ux (x0 , y0 , n1t )  n  = −ζ u + ζj unx,j  0 x ,0  ∂ x∂ y   j =1     N  2   ∂ Uy (x0 , y0 , n1t )   = −ζ0 uny,0 + ζj uny,j  ∂ x∂ y j =1

(5)

where N is the number of nodes in the star whose central node has the coordinates (x0 , y0 ) (in this work N = 8 and the star nodes are selected by using the four quadrant criteria [5]). m0 , η0 , ζ0 are the coefficients that multiply the approximate values of the functions Ux and Uy at the central node for the time n1t (unx,0 and uny,0 respectively) in the generalized finite difference explicit expressions for the space derivatives. mj , ηj , ζj are the coefficients that multiply the approximate values of the functions Ux and Uy at the rest of the star nodes for the time n1t (unx,j and uny,j respectively) in the generalized finite difference explicit expressions for the space derivatives. The replacement in Eq. (1) of the explicit expressions obtained for the partial derivatives leads to the explicit difference scheme

      N N    n+1 n −1 n 2 2 n n 2 n n   ux,0 = 2ux,0 − ux,0 + (1t ) α −m0 ux,0 + m j ux , j + β −η0 ux,0 + η j ux , j    1 1        N    2 2 n n  + (α − β ) −ζ0 uy,0 + ζ j uy , j    1

      N N     n+1 n−1 n 2 2 n n 2 n n  uy,0 = 2uy,0 − uy,0 + (1t ) β −m0 uy,0 + mj uy,j + α −η0 uy,0 + ηj uy,j     1 1       N     2 2 n n  + (α − β ) −ζ0 ux,0 + ζ j ux , j . 

(6)

1

As we are using an explicit method [15], we have studied the stability and obtained the star stability condition in [16]. In this paper, the condition for the stability of the star has been established as:

   1t < 

4

(α 2 + β 2 )[(|m0 | + |η0 |) +

 . (m0 + η0 )2 + ζ02 ]

We have also studied the dispersion for the phase and group velocities in [16] in depth.

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

43

3. Recursive equations

3.1. Recursive equations with PML in x-direction For computational convenience, we split the second order equations of motion Eq. (1) into five coupled first order equations by introducing the new field variables γxx , γxy , γyy

 ∂ Ux (x, y, t ) ∂γxx (x, y, t ) ∂γxy (x, y, t )  = + ρ   ∂ t ∂x ∂y      ∂γ ( x , y , t ) ∂γ ( ∂ U ( x , y , t )  xy yy x, y, t ) y   = + ρ   ∂t ∂x ∂y     ∂γ (x, y, t ) ∂ U ( x , y , t ) ∂ Uy (x, y, t ) xx x = (λ + 2µ) +λ  ∂t ∂x ∂y      ∂ U ( x , y , t ) ∂ U ( x , y , t ) ∂γ ( x , y , t ) x y   xy =µ +µ   ∂ t ∂ y ∂x      ∂γ ( x , y , t ) ∂ U ( x , y , t ) ∂ Uy (x, y, t )  x  yy =λ + (λ + 2µ) . ∂t ∂x ∂y

(7)

We shall make two simplifications and assume that the space far from the region of interest is homogeneous, linear and time invariant. Then, under these assumptions, the radiating solution in infinite space must be (superposition of plane waves):

ω(x, t ) = W (x, t )ei(κ·x−wt ) .

(8)

As w is an analytic function of x, we can analytically continue it, evaluating the solution at complex values of x. Then, the solution is not changed in the region of interest and the reflections are avoided:

 Ux (x, y, t ) = ux (x, y)e−iwt ⇒ U˙ x (x, y, t ) = −iw ux (x, y)e−iwt = −iw Ux (x, y, t )     −i w t  ⇒ U˙ y (x, y, t ) = −iw uy (x, y)e−iwt = −iw Uy (x, y, t )  Uy (x, y, t ) = uy (x, y)e γxx (x, y, t ) = Γxx (x, y)e−iwt ⇒ γ˙xx (x, y, t ) = −iwΓxx (x, y)e−iwt = −iwγxx (x, y, t )    γxy (x, y, t ) = Γxy (x, y)e−iwt ⇒ γ˙xy (x, y, t ) = −iwΓxy (x, y)e−iwt = −iwγxy (x, y, t )    γyy (x, y, t ) = Γyy (x, y)e−iwt ⇒ γ˙yy (x, y, t ) = −iw Γyy (x, y)e−iwt = −iwγyy (x, y, t ).

(9)

Thus, we have a complex coordinate x˜ = x + if .

(10)

As this complex coordinate is inconvenient, we have a change of variables in this region (PML)

  df ∂ x˜ = 1 + i ∂ x. dx

(11)

In order to have an attenuation rate in the PML independent of frequency (ω), we have chosen df dx

=

δx (x) ω

(12)

where ω is the angular frequency and δx is a function of x. PML x-dir can be conceptually assumed up by a single transformation of the original equation. Then wherever an x derivative appears in the wave equations, it is replaced in the form

∂ 1 ∂ → . δ (x) ∂x 1 + i xw ∂ x

(13)

44

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

The equations are frequency-dependent, and to avoid it a solution is to use an auxiliary differential equation (ADE) approach in the implementation of PML. The following equations are obtained

   1 ∂γxx (x, y, t ) ∂γxy (x, y, t ) ∂ Ux (x, y, t )   = + + ψ1 (x, y, t ) − δx Ux (x, y, t )   ∂t ρ ∂x ∂y        1 ∂γxy (x, y, t ) ∂γyy (x, y, t ) ∂ Uy (x, y, t )    = + + ψ2 (x, y, t ) − δx Uy (x, y, t )   ∂t ρ ∂x ∂y      ∂ Ux (x, y, t ) ∂ Uy (x, y, t )  ∂γxx (x, y, t )  = (λ + 2µ) +λ + ψ3 (x, y, t ) − δx γxx (x, y, t )    ∂t ∂x ∂y      ∂γxy (x, y, t ) ∂ Ux (x, y, t ) ∂ Uy (x, y, t )   =µ +µ + ψ4 (x, y, t ) − δx γxy (x, y, t )    ∂t ∂y ∂x      ∂ Ux (x, y, t ) ∂ U y ( x, y , t ) ∂γ (x, y, t )   yy =λ + (λ + 2µ) + ψ5 (x, y, t ) − δx γyy (x, y, t )  ∂t ∂x ∂y  ∂ψ1 (x, y, t ) δx ∂γxy (x, y, t )    =   ∂ t ρ ∂y      ∂ψ2 (x, y, t ) δx ∂γyy (x, y, t )   =    ∂ t ρ ∂y      ∂ Uy (x, y, t ) ∂ψ3 (x, y, t )   = λδx   ∂ t ∂y      ∂ψ4 (x, y, t ) ∂ Ux (x, y, t )    = µδx   ∂ t ∂y      ∂ψ5 (x, y, t ) ∂ Uy (x, y, t )   = (λ + 2µ)δx ∂t ∂y

(14)

where the five last equations in Eq. (14) are ADE approach and the new field variables

 1 δx ∂γxy (x, y, t )  ψ1 (x, y, t ) = i   ρ ω ∂y      1 δx ∂γyy (x, y, t )    ψ2 (x, y, t ) = i   ρ ω ∂y     δx ∂ Uy (x, y, t ) ψ3 (x, y, t ) = iλ  ω ∂y      δx ∂ Ux (x, y, t )   ψ4 (x, y, t ) = iµ   ω ∂y       ψ (x, y, t ) = i(λ + 2µ) δx ∂ Uy (x, y, t ) . 5 ω ∂y

(15)

3.1.1. A scheme in GDFM for domain interest Following [5,6,16], the explicit difference formulae for the first order spatial derivatives of a function are

 N  ∂ Ux (x0 , y0 , n1t )   = −m1,0 unx,0 + m1,j unx,j    ∂x j=1 N    ∂ Ux (x0 , y0 , n1t )   = −m2,0 unx,0 + m2,j unx,j  ∂y j=1

and similarly for first spatial derivatives of the functions: Uy , γxx , γxy , γyy , ψ1 , ψ2 , ψ3 , ψ4 , ψ5 .

(16)

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

45

Substituting Eq. (16) into Eq. (7) the scheme in GFDM for elastic part is obtained

   N N   1t  n +1 n n n n n  −m1,0 γxx,0 + m1,j γxx,j − m2,0 γxy,0 + m2,j γxy,j ux , 0 = ux , 0 +    ρ  j =1 j =1        N N     n n n un+1 = un + 1t −m γ n +  m1,j γxy,j − m2,0 γyy,0 + m2,j γyy,j 1,0 xy,0  y,0 y,0  ρ  j =1 j=1           N N    γxxn+,01 = γxxn ,0 + 1t (λ + 2µ) −m1,0 unx,0 + m1,j unx,j + λ −m2,0 uny,0 + my,j un2,j   j = 1 j = 1           N N     n + 1 n n n n n   γxy,0 = γxy,0 + 1t µ −m2,0 ux,0 + m2,j ux,j + µ −m1,0 uy,0 + m1,j uy,j    j =1 j =1           N N     n + 1  n n n n n  m1,j ux,j + (λ + 2µ) −m2,0 uy,0 + m2,j uy,j . γyy,0 = γyy,0 + 1t λ −m1,0 ux,0 + j =1

(17)

j =1

3.1.2. A scheme in GDFM for PML part Substituting Eqs. (16) into Eqs. (14) the scheme in GFDM for PML part is obtained

   N N   1t n+1 n n n n n   −m1,0 γxx,0 + m1,j γxx,j − m2,0 γxy,0 + m2,j γxy,j + 1t [ψ1n,0 − δx unx,0 ] ux , 0 = ux , 0 +   ρ  j =1 j =1        N N     1 t  n +1 n n n n n n n  m1,j γxy m2,j γyy  ,j − m2,0 γyy,0 + ,j + 1t [ψ2,0 − δx uy,0 ] uy,0 = uy,0 + ρ −m1,0 γxy,0 +   j = 1 j = 1           N N     n + 1 n n n n n  γxx,0 = γxx,0 + 1t (λ + 2µ) −m1,0 ux,0 + m1,j ux,j + λ −m2,0 uy,0 + m2,j uy,j     j=1 j=1     + 1t [ψ3n,0 − δx γxxn ,0 ]            N N     n +1 n n n n n   γxy,0 = γxy,0 + 1t µ −m2,0 ux,0 + m2,j ux,j + µ −m1,0 uy,0 + m1,j uy,j + 1t [ψ4n,0 − δx γxyn ,0 ]    j = 1 j = 1           N N     γ n+1 = γ n + 1t λ −m un + n n n  m1,j ux,j + (λ + 2µ) −m2,0 uy,0 + m2,j uy,j 1 ,0 x ,0  yy,0   yy,0 j =1 j =1 + 1t [ψ5n,0 − δx γyyn ,0 ]        N    1t  n+1 n n n  δx −m2,0 γxy,0 + m2,j γxy,j ψ1,0 = ψ1,0 +    ρ  j =1        N   1t  n +1 n n n   ψ2,0 = ψ2,0 + δx −m2,0 γyy,0 + m2,j γyy,j   ρ  j =1        N    ψ n+1 = ψ n + λ1t δ −m un + n  m2,j uy,j x 2 ,0 y ,0  3,0   3 ,0  j =1       N      ψ4n,+0 1 = ψ4n,0 + µ1t δx −m2,0 unx,0 + m2,j unx,j     j = 1        N    n + 1 n n n  ψ5,0 = ψ5,0 + (λ + 2µ)1t δx −m2,0 uy,0 + m2,j uy,j . j =1

(18)

46

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

3.2. Recursive equations with PML in x-direction and y-direction In this case

 ∂     → ∂x  ∂    → ∂y

  ∂ δ −1 1+i ∂x ω   ∂ δ −1 1+i . ∂y ω

(19)

We obtain

   ∂ Ux (x, y, t ) 1 ∂γxx (x, y, t ) ∂γxy (x, y, t )  = + − δ Ux (x, y, t )   ∂t ρ ∂x ∂y        ∂ Uy (x, y, t ) 1 ∂γxy (x, y, t ) ∂γyy (x, y, t )   = + − δ U y ( x, y , t )   ∂t ρ ∂x ∂y    ∂γxx (x, y, t ) ∂ Ux (x, y, t ) ∂ Uy (x, y, t ) = (λ + 2µ) +λ − δγxx (x, y, t )  ∂ t ∂ x ∂y     ∂τxy (x, y, t ) ∂ Ux (x, y, t ) ∂ Uy (x, y, t )   =µ +µ − δγxy (x, y, t )    ∂ t ∂ y ∂x      ∂γyy (x, y, t ) = λ ∂ Ux (x, y, t ) + (λ + 2µ) ∂ Uy (x, y, t ) − δγ (x, y, t ). yy ∂t ∂x ∂y

(20)

3.2.1. A scheme in GDFM for domain of interest The equations of the elastic part are given by Eq. (17). 3.2.2. A scheme in GDFM for PML part Substituting Eqs. (16) into Eqs. (20) the scheme in GFDM for PML part is obtained

   N N   1t  n+1 n n n n n   m1,j γxx,j − m2,0 γxy,0 + m2,j γxy,j − 1t δ unx,0 ux,0 = ux,0 + ρ −m1,0 γxx,0 +   j = 1 j = 1       N N     1t n+1  n n n n n uy,0 = uy,0 + m1,j γxy,j − m2,0 γyy,0 + m2,j γyy,j − 1t δ uny,0 −m1,0 γxy,0 +   ρ   j=1 j =1         N N    n +1 n n n n n γxx,0 = γxx,0 + 1t (λ + 2µ) −m1,0 ux,0 + m1,j ux,j + λ −m2,0 uy,0 + my,j uy,j − 1t δγxxn ,0   j =1 j =1          N N     n + 1 n n n n n   m2,j ux,j + µ −m1,0 uy,0 + m1,j uy,j − 1t δγxyn ,0 γxy,0 = γxy,0 + 1t µ −m2,0 ux,0 +    j=1 j =1          N N     n +1 n n n n n   m1,j ux,j + (λ + 2µ) −m2,0 uy,0 + m2,j uy,j − 1t δγyyn ,0 . γyy,0 = γyy,0 + 1t λ −m1,0 ux,0 + j =1

(21)

j =1

4. Numerical results 4.1. Results for several wavelengths Table 1 shows the values of global error for several analytical solutions of the problem defined by Eq. (1), when we assign different value parameters A and B in Eq. (25). The domain is Ω = [0, 1] × [0, 1] ⊂ R2 with Dirichlet boundary conditions and initial conditions given by:

 Uy (x, y, 0) = cos x cos y Ux (x, y, 0) = sin x sin y; ∂ Uy (x, y, 0) ∂ Ux (x, y, 0)  = 0; = 0. ∂t ∂t

(22)

We use an irregular mesh with 121 nodes (see Fig. 1) with an Index of Irregularity of IIC = 0.8944 and 1t = 0.01. The solution included in Table 1 is for time step 500.

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

47

Fig. 1. Dispersion of the waves P in the irregular mesh with IIC = 0.8944. Table 1 Global errors for an irregular mesh with 121 nodes, IIC = 0.8944, n = 500 and 1t = 0.01. Analytical sol. √

A A A A A A

= √ ,B = = √ 2, B = 1 = 2 2, B = √ 2 = 4π, B = 2√2π = 8π, B = 4 √ 2π = 16π, B = 8 2π 2 2

1 2

Error global Ux

Error global Uy

1.232 × 10 1.646 × 10−6 4.081 × 10−4 9.188 × 10−2 3.035 × 10−1 4.942 × 10−1

2.443 × 10−5 1.778 × 10−6 2.001 × 10−4 8.647 × 10−2 3.275 × 10−1 4.999 × 10−1

−5

The weighting function is

w(hjx , hjy ) = 

1 h2jx + h2jy

(23)

3

and the criterion for the selection of star nodes is the quadrant criterion [5–7]. The global error is evaluated for the last time step considered, using the following formula



NT 

(sol(j)−exac(j))2

j=1

Global error =

NT

(24)

|exacmax |

where sol(j) is the GFDM solution at the node j, exac(j) is the exact value of the solution at the node j, exacmax is the maximum value of the exact solution in the cloud of nodes considered and NT is the total number of nodes of the domain. The analytical solutions are: Ux (x, y, t ) = cos(Aβ t ) sin(Bx) sin(By);

Uy (x, y, t ) = cos(Aβ t ) cos(Bx) cos(By).

(25)

4.2. Dispersion and irregularity Fig. 1 shows the dispersion of the waves P in each node of the irregular mesh with IIC = 0.8944. 4.3. GFDM with PML Let us solve Eq. (1), in Ω = [0, 2] × [0, 1] ⊂ R2 , with homogeneous Dirichlet boundary conditions and the initial conditions are given √ by Eq. (22), using the regular mesh with 861 nodes (see Fig. 2), the analytical solutions are given by Eq. (25) with A = 2 and B = 1 (see Fig. 3). The weighting function is given by Eq. (23) and the criterion for the selection of star nodes is the quadrant criterion. Fig. 5 depicts the approximated solution of ux , after 100 time steps, with PML in x-direction and y-direction for 1.4 ≤ x ≤ 2 and 0.6 ≤ y ≤ 1 (see Fig. 4) and homogeneous boundary conditions. Figs. 7–9 show the approximated solution of ux , after 5, 20 and 100 time steps respectively, with PML in x-direction and y-direction for ≤ x ≤ 0.6 and 0 ≤ y ≤ 0.2, 0.8 ≤ y ≤ 1 (see Fig. 6). Fig. 11 shows the approximated solution of ux , after 100 time steps, with PML in x-direction and y-direction for ≤ x ≤ 0.6 and 0 ≤ y ≤ 0.2, 0.8 ≤ y ≤ 1 and for a heterogeneous region of interest defined, by zi (i = 1, 2, 3) in Fig. 10, and with the properties defined in Table 2.

48

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

Fig. 2. Regular mesh (861 nodes).

Fig. 3. Exact solution Ux without PML. 1

0.8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 4. Regular mesh with PML region.

Fig. 5. Approximated solution Ux with PML.

5. Conclusions This paper shows a scheme in generalized finite differences to the solution of the problem of seismic wave propagation in 2-D in homogeneous and heterogeneous media, including PML absorbing boundaries in the model. In the paper, we analyse error variation with the number of nodes on the wave-length. From these first analyses we obtain errors of the order of one thousand with variations that comprise around twelve nodes per wave-length and meshes which are not too irregular (with an IIC value of nearly one).

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

49

1

0.8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 6. Regular mesh with PML region.

Fig. 7. Approximated solution Ux with PML, after 5 time steps.

Fig. 8. Approximated solution Ux with PML, after 20 time steps.

Fig. 9. Approximated solution Ux with PML, after 100 time steps.

In general, a regular mesh produces the least error. This error is similar to the error obtained using the standard Finite Difference Method. For this reason it is advisable to use a regular mesh as a base and change the position of some nodes or include more nodes if the shape of the domain or the solution requires it. The error increases when the irregularity of the

50

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51 1

0.8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 10. Regular mesh with PML region and heterogeneous dates.

Fig. 11. Approximated solution Ux with PML. Table 2 Properties of the region.

z1 z2 z3

λ

µ

ρ

α

β

0.5 0.7 0.8

0.25 0.4 0.5

1.0 0.2 1.1

1.0 2.738 1.279

0.5 1.414 0.674

mesh is greater. In the event of this happening we have to use smaller time-steps, therefore the introduction of irregularities in the mesh should be carried out in a suitable manner and it must be totally justified. The formulation of the PML is compatible with GFDM. In this paper, we have obtained the recursive equations for one and two directions. Numerical results confirm that PML in absorbing outgoing waves performs extraordinarily as can be appreciated in the academic test included in this paper, where it is possible to see how displacements rapidly decrease in these zones. Finally, we have used heterogeneous formulation for heterogeneous problem because homogeneous ones, following [9], can be complicated in seismological practice and we have obtained good results dealing with it. Acknowledgement The authors acknowledge the support from Ministerio de Ciencia e Innovación of Spain, project CGL2008-01757/CLI. References [1] P.S. Jensen, Finite difference technique for variable grids, Computers & Structures 2 (1972) 17–29. [2] N. Perrone, R. Kao, A general finite difference method for arbitrary meshes, Computers & Structures 5 (1975) 45–58. [3] T. Liszka, J. Orkisz, The finite difference method at arbitrary irregular grids and its application in applied mechanics, Computers & Structures 11 (1980) 83–95. [4] J. Orkisz, Finite difference method (Part, III), in: M. Kleiber (Ed.), Handbook of Computational Solid Mechanics, Spriger-Verlag, Berlin, 1998. [5] J.J. Benito, F. Ureña, L. Gavete, Influence several factors in the generalized finite difference method, Applied Mathematical Modelling 25 (2001) 1039–1053. [6] J.J. Benito, F. Ureña, L. Gavete, R. Alvarez, An h-adaptive method in the generalized finite difference, Computer Methods in Applied Mechanics and Engineering 192 (2003) 735–759. [7] J.J. Benito, F. Ureña, L. Gavete, B. Alonso, Application of the generalized finite difference method to improve the approximated solution of PDEs, Computer Modelling in Engineering & Sciences 38 (2009) 39–58.

J.J. Benito et al. / Journal of Computational and Applied Mathematics 252 (2013) 40–51

51

[8] J.J. Benito, F. Ureña, L. Gavete, Leading-Edge Applied Mathematical Modelling Research, Nova Science Publishers, New York, 2008 (Chapter 7). [9] P. Moczo, Introduction to Modelling Seismic Wave Propagation by Finite-Difference Method, in: Lectures Notes, Kyoto, 1998. [10] P. Moczo, J. Kristek, M. Galis, P. Pazak, M. Balazovjech, The finite-difference and finite-element modeling of seismic wave propagation and earthquake motion, Acta Physica Slovaca 57 (2) (2007) 177–406. [11] J.P. Berenger, A perfectly matched layer for the absortion of electromagnetic, Journal of Computational Physics 114 (1994) 185–200. [12] W.C. Chew, Q.H. Liu, Perfectly matched layer for elastodynamics: a new absorbing boundary condition, Journal of Computational Acoustics 4 (1996) 341–359. [13] S.G. Jonshon, Notes on perfectly matched layers (PMLs). Courses 18.369 and 18.336 at MIT, 2008. [14] E.A. Skelton, S.D.M. Adams, R.V. Craster, Guided elastic waves and perfectly matched layers, Wave Motion 44 (2007) 573–592. [15] J.J. Benito, F. Ureña, L. Gavete, B. Alonso, Solving parabolic and hyperbolic equations by generalized finite difference method, Journal of Computational and Applied Mathematics 209 (2) (2007) 208–233. [16] F. Ureña, J.J. Benito, E. Salete, L. Gavete, A note on the application of the generalized finite difference method to seismic wave propagation in 2-D, Journal of Computational and Applied Mathematics 236 (12) (2011) 3016–3025.