Open access

Perfectly Matched Layer for Finite Element Analysis of Elastic Waves in Solids

Written By

Koji Hasegawa and Takao Shimada

Published: 10 October 2012

DOI: 10.5772/46162

From the Edited Volume

Finite Element Analysis - Applications in Mechanical Engineering

Edited by Farzad Ebrahimi

Chapter metrics overview

6,022 Chapter Downloads

View Full Metrics

1. Introduction

Numerical analysis of scattering and propagation of elastic waves in solids gives insight into physical phenomena under operation of ultrasonic devices such as electromechanical filters and resonators, nondestructive testing with ultrasonic waves and seismic prospecting. To take anisotropy of solids and complex structures of composite solids into account, commercial simulator based on the finite element method are available. Numerical models for finite element analysis (FEA) must be bounded and infinite half spaces of models should be replaced with finite domains and absorbing boundary conditions.

Perfectly matched layer (PMLs) is one of popular absorbing boundary conditions for truncating the computational domain of open regions without reflection of oblique incident waves. In 1994, Berenger invented a PML for electromagnetic waves in the finite difference time domain (FD-TD) method by a splitting field method.[1] Because fields in Berenger’s PML do not satisfy the Maxwell’s equations, two concepts have been introduced for implementation in the finite element method (FEM) of electromagnetic wave problems: the analytic continuation or the complex coordinate stretching[2, 3](CCS) and anisotropic PMLs.[4] Nowadays PMLs for electromagnetic waves are widely used in the FD-TD method and the FEM.

Extension of PMLs to elastic waves in isotropic solids in the Cartesian coordinate first appeared in 1996.[5, 6] In the cylindrical and spherical coordinates, PMLs were presented by using splitting field method in isotropic solids in 1999[7] and by using analytic continuation in anisotropic solids in 2002.[8] Recently validity and usefulness of PMLs derived from the analytic continuation in piezoelectric solids was demonstrated. [9, 10, 11] Hastings et. al.[5] reported better performance of PMLs by the FD-TD method than the second-order absorbing boundary condition (ABC) of Peng and Toksöz: in the range of the incident angle from 0° to 80°, reflection powers of S- or P-wave, which is excited by a pure S- or P-wave line source and propagating in a two-dimensional infinite isotropic solid modeled by a rectangular solid with its opposite sides attached sponge mediums and other sides imposed ABC or loaded PMLs, from the PML side are suppressed below -45 and -80 dB with 8 and 16 grid spaces of the PML region, respectively. On the other hand, reflection power level at the computational domain edge imposed ABC is in the range of -90 dB to -10 dB. This implies that PMLs yield more superior approximation of perfect matching than the ABC with thickening PML and increasing the number of grid spaces.

We recommend that readers who are unfamiliar with PMLs consult Basu and Chopra[12] about explanations and finite element (FE) implementation of PMLs for time-harmonic elastodynamics, Michler et al.[13] about derivation of material constants of PML in FE method by analytic continuation, and Taflove and Hagness[14] about PMLs for electromagnetic waves in FD-TD method.

Although PML is one of attractive artificial materials, two questions of PMLs derived from the analytic continuation are left: why are the particle displacements in the complex coordinate identical to those in the real coordinate and why must we multiply stress tensors by the Jacobian of the coordinate transformation?

For replying to the questions, we will examine a derivation of PMLs for elastic waves in the Cartesian, the cylindrical and the spherical coordinates from the differential form on manifolds. Our results reveal that the components of stress tensors and the particle displacement vectors in the analytic continuation are not transformed to the real space.[15] In addition, the rule for determining PML parameters in the Cartesian coordinate holds in the cylindrical and spherical coordinates.[16]

Mathematical models of PMLs, which are given by differential equations and boundary conditions, are exactly perfect matching medium. In numerical models, however, discretizing PMLs changes phase velocities of propagating waves and generates reflection waves from the PML region.[17] Furthermore, approximation of infinite regions with finite thick layers also generates reflection waves from the PML terminal.[1, 17, 18]

Estimating matching performance and optimizing parameters of PMLs in a numerical domain are required before solving problems. Chew and Jin investigated dependence of PML’s performance on attenuation parameters of FE analysis of electromagnetic wave problems.[18] For FD-TD method Collino and Monk also carried out such an investigation.[19] Recently, Bermúdez et al. investigated absorbing functions for time harmonic Helmholtz equations in the Cartesian and cylindrical coordinates under the condition of ignoring reflection caused by FE-discretization and showed the advantages of non-integrable absorbing functions over conventional functions of power series.[20, 21]

Most of these investigations of optimizing attenuation parameters of PML employed numerical analysis of scattering problems in the two dimensions such as plane or cylindrical wave scattering problems. For tackling optimization problem of PML parameters, plane wave scattering problem is appropriate because required resource of computation is small. For FEA, Chew and Jin[18] modeled scattering of plane waves as electromagnetic field analysis in the thick layer in the one dimension. But this model has not been applied to elastic wave scattering.

In this chapter, we also examine PML performance of FE-models in the frequency range with scattering problems of elastic waves in an isotropic solid as field analysis in the thick layer in the one dimension. To the best of our knowledge, quantification of reflection power generated by FE-discretization has not attracted attention. Recently, for electromagnetic waves, we reported that the reflection power caused by discretization can be computed by the equivalent transmission line with its impedance and propagation constant determined by discretized wave numbers.[15] Because, for elastic waves, the transfer matrix is popular, we explain the reflection from PMLs by the transfer matrix of elastic waves, and confirm that numerical results of FE-models may be predicted by replacement of propagation constants of elastic waves in PML with discretized wave numbers.

Advertisement

2. Derivation of perfectly matched layers for elastic waves by using complex coordinate stretching and differential form

2.1. Differential form

A particle displacement vectoru, particle velocity vectorv, density of momentaP, stress tensor T¯¯ and displacement gradient tensor F¯¯ are given as follows:

u=uixi,E1
v=vixi,E2
P=13!Pαβγixidxαdxβdxγ,E3
T¯¯=12Tαβixidxαdxβ,E4
F¯¯=Fαixidxα=du,E5

where xi and dxi(i=0,1,2) are the contravariant and covariant basis vectors, and represent the tensor product and the cross product, respectively. dis the exterior differential operator. Newton’s equation of motion is

dT¯¯=Pt.E6

Changing the coordinate gives relations of tensor components: for a tensor with a tensor type of the contravariant of rank 1 and the covariant of rank q, V=VXα1αqiXidXα1dXαq=Vxβ1βqkxkdxβ1dxβq, the relation of tensor components is

VXα1αqi=Xixkxβ1Xα1xβqXαqVxβ1βqk.E7

Using CCS[2, 3, 8] given by Xi=s˜i(τ)dτ=s˜iR(τ)+js˜iI(τ)dτ with the two real functions s˜iR(τ) ands˜iI(τ), we have the relation

VXα1αqi=s˜i(xi)[s˜α1(xα1)s˜αq(xαq)]1Vxα1αqi.E8

Here j is the imaginary unit.

2.2. PMLs in the Cartesian, the cylindrical and the spherical coordinates

In the complex coordinate stretching (CCS), we consider that the real coordinate (x0,x1,x2) is(x,y,z), (r,θ,z)or (r,θ,ϕ) for the Cartesian, the cylindrical or the spherical coordinate, respectively. Assuming that the same constitutive equations in the real Cartesian, cylindrical and spherical coordinate exist in the complex coordinate, (X0,X1,X2), we have

Pc=ρvc,Tijc=CijklSklc=Cijkl(Fklc+Flkc)/2E9
=CijklFklc.E10

Here, the superscript c denotes the value in the complex coordinate and the mass density ρ and the stiffness Cijkl(i,j,k,l=X0,X1,X2) are the values corresponding to original material parameters, mass density and stiffness constants, of its PML in the Cartesian, the cylindrical and the spherical coordinates. Using eq. (8) to eqs.(2)-(5), and recalling that the base differentials dξ,dη and dζ of the general orthogonal coordinate system (ξ,η,ζ) are dual to the unit vectors ξ^hξ,η^hη andζ^hζ, we have

vic=sivi(nosummation),E11
Pic=sis0s1s2Pi(nosummation),E12
Tijc=sisjs0s1s2Tij(nosummation),E13
Fijc=sisjFij(nosummation).E14

Here si=hichirs˜i with hir and hic being scale factors of general orthogonal coordinate systems (x0,x1,x2) and(X0,X1,X2), respectively. Note that the scale factors hi are given by follows: in the cylindrical coordinate(r,θ,z)h0=1,h1=r,h2=1, and in the spherical coordinate(r,θ,ϕ)h0=1,h1=r ,h2=rsinθ. In addition, in the Cartesian coordinate,h0=h1=h2=1.

The quotient rule and eqs. (9)- (14) yield PML material constants: the mass density ρPML is

ρPML=s0s1s2ρE15

and the stiffness is

CijklPML=s0s1s2sksisjslCijkl(nosummation).E16

Here, s0=s˜0, s1=s˜1Rr, s2=s˜2in the cylindrical coordinate system (r,θ,z) with its complex coordinate(R,Θ,Z), ands0=s˜0, s1=s˜1Rr, s2=s˜2RsinΘrsinθin the spherical coordinate system (r,θ,ϕ) with its complex coordinate(R,Θ,Φ). In addition, si=s˜iin the Cartesian coordinate system.

Eqs. (15) and (16) show that PML parameters for elastic waves in solids in the cylindrical and spherical coordinates may be calculated by the same procedure in the Cartesian coordinate.

2.3. Derivation of PML constants by the analytic continuation

For simplicity, we present a procedure of deriving material constants in only cylindrical coordinates by the analytic continuation[8] below. Note that in spherical coordinates the same procedure may be applied. We recommend that the reader who is interesting in the procedure consult Zheng and Huang.[8]

First we consider Newton’s equation of motion. In a cylindrical coordinate(r,θ,z), the governing equations are

ρω2ur=1rr(rTrr)+1r(Trθθ+r^θ^θTθθ)+Trzz,E17
ρω2uθ=1rr(rTθr)+1r(θ^r^θTrθ+Tθθθ)+Tθzz,E18
ρω2uz=1rr(rTzr)+1rTzθθ+Tzzz.E19

Here, we use phasor notation. The time dependences of the fields are exp(jωt) where ω is angular frequency. Applying CCS with a complex coordinate(R,Θ,Z), multiplying the CCS equations by s0s1s2 and using the assumption ofuic=ui, R^=r^, Θ^=θ^andZ^=z^, we get the governing equations in the PML:

ρPMLAω2ur=1rr(rs1s2Trrc)+1r((s0s2Trθc)θ+r^θ^θ(s0s2Tθθc))+(s0s1Trzc)z,E20
ρPMLAω2uθ=1rr(rs1s2Tθrc)+1r(θ^r^θ(s0s2Trθc)+(s0s2Tθθc)θ)+(s0s1Tθzc)z,E21
ρPMLAω2uz=1rr(rs1s2Tzrc)+1r(s0s2Tzθc)θ+(s0s1Tzzc)z.E22

Here, the mass density of PML is defined asρPMLA=s0s1s2ρ. When we rewrite components of stress tensors in the PMLs asTijPMLA=s0s1s2sjTijc, we may identify eqs.(20) ~ (22) to eqs.(17) ~ (19).

Next we consider the displacement gradient uc=[Γklc] in the complex coordinate(R,Θ,Z). Using definition du=u(x^ihidxi) and the assumption ofuic=ui, R^=r^, Θ^=θ^, Z^=z^, and applying CCS with a complex coordinate(R,Θ,Z), we have the relation:

uc=[uRcR1R(uRcΘ+R^Θ^ΘuΘc)uRcZuΘcR1R(Θ^R^ΘuRc+uΘcΘ)uΘcZuZcR1RuZcΘuZcZ]=[1s0urr1s11r(urθ+r^θ^θuθ)1s2urz1s0uθr1s11r(θ^r^θur+uθθ)1s2uθz1s0uzr1s11ruzθ1s2uzz].E23

Hence we haveΓklc=1slΓkl.

Using the quotient rule and the constitutive equation, Tijc=CijklΓklc, we get the constitutive equation of the PML in the real coordinate(r,θ,z):TijPMLA=s0s1s2sjslCijklΓkl. Therefore, we may define the stiffness of the PML derived by the analytic continuation: CijklPMLA=s0s1s2sjslCijkl

2.4. Comparison with PML material constants derived from differential forms and the analytic continuation

By the analytic continuation, Zheng and Huang[8] derived the mass density and stiffness of PML in the cylindrical and spherical coordinates: ρPMLA=s0s1s2ρandCijklPMLA=s0s1s2sjslCijkl. The mass density agree with our result, eq. (15), because multiplying the stress tensors by the Jacobian of the coordinate transformation, s0s1s2, adjusts the mass density. We note that the form of eq. (15) is also derived from eq. (6) with the tensor type of mass density being covariant of rank 3, i.e. 3-form. The stiffness is different from eq. (16) because in the analytic continuation, the manipulation of the coordinate transformation corresponding to the part of stress tensor and the particle displacement vector, contravariant of rank 1, is excluded. This fact can be confirmed by the derivation procedure presented in the previous section for the cylindrical coordinate:we put TijPMLA=s0s1s2sjTijc and use the assumptionui=uic.

To show a difference between PML material constants, we consider an isotropic solid with following stiffness constants in the Cartesian coordinate(x0,x1,x2):Cijkl=λδijδkl+μ(δikδjl+δilδjk). Here, λand μ are the Lamé constants of an isotropic solid, the subscripts i,j,k and l denote the xi-, xj-, xk- and xl-axis, respectively, and δij is the Kronecker delta. Components of stiffness tensors derived from the differential form and analytic continuation, CijklPMLandCijklPMLA, respectively, are given by

CijklPML=[(λsi2δijδkl+μsj2δikδjl+μsi2δilδjk]s0s1s2(nosummation),E24
CijklPMLA=[λsiskδijδkl+μsj2δikδjl+μsisjδilδjk]s0s1s2(nosummation).E25

Table 1 shows all components of the stress tensor computed with CijklPML andCijklPMLA. CijklPMLAgives T˜ijT˜ji(ij) and we predict that rotational forces may be observed. WithCijklPML, however, we have a symmetric stress tensor,Tij=Tji(ij).

Table 1.

Components of a stress tensor in a PML material of an isotropic solid in the Cartesian coordinate.

Advertisement

3. Reflection from PMLs discretized for finite element models in the frequency domain

We consider a plane elastic wave propagating in a half infinite isotropic solid attached with its PML backed with a vacuum region as shown in Fig.1. Here θ is the incident angle, θpand θs are propagation angles of P-waves and SV- or SH-waves, Lis thickness of the PML, kiand kr,m(m=0,1,2) are wave vectors of the incident wave and reflected P-, SV- and SH-waves, respectively.

We use the phasor notation and assume that the time dependences of all fields areexp(jωt), where j is the imaginary unit and ω is the angular frequency.

Figure 1.

Reflection by the plane boundary between an isotropic solid and its PML.

When the stiffness component of the isotropic solid Cijkl(i,j,k,l=x,y,z) is given by Cijkl=λδijδkl+μ(δikδjl+δjlδik) where λ and μ are the Lamé constants and δij is the Kronecker delta, the stiffness component of its PML CijklPML is

CijklPML=(λsi2δijδkl+μsj2δikδjl+μsi2δilδjk)sxsysz.E26

Here si(i=x,y,z) is a coordinate stretching factor of i-direction.[2] The mass density of the PML ρPML is given by

ρPML=sxsyszρ.E27

Here ρ is the mass density of the isotropic solid. For examining absorbing performance of PMLs in the x direction, taking an assumption of considering fields being consisted by plane waves propagating on the xy plane, we have a differential equation in one variablex: from Newton’s equation of motion and constitutive equation we get the differential equation in the PML

CijklPMLxj(ukxl)=ω2ρPMLuiE28

where ui is the component of the particle displacement in the i-direction (i=x,y,z).

In this case, we may choose the coordinate stretching factor as follows:

sx=1jsxI(x),sy=sz=1.E29

Here sxI(x) is the imaginary part of sx and therefore a real function, which controls absorbing performance of propagating waves in PMLs.

Boundary conditions at the interface of isotropic solid and PML, x=0, are the nonslip condition and the continuous condition of the normal component of the stress:[15]

ui(0)=siui(+0),E30
Tix(0)=sisxsxsyszTix(+0),i=x,y,z.E31

At the terminal of PML, x=L, the boundary condition is

sisxsxsyszTix(L)=0,i=x,y,z.E32

3.1. Numerical procedure

3.1.1. Finite element analysis

Because finite element formulation of a thick plate with line elements as shown in Fig. 2 is well known and we use COMSOL MultiPhysics for FEA, we explain the Robin condition at x=0 and a formula for reflection coefficients.

Figure 2.

Line element with (n+1)-nodes.

In the half isotropic solid, the field distribution, components of the particle displacement and stress, can be expressed by superposition of incident and reflected plane waves:

u(r)=lRlul(r)ejkr,lr+uiejkir.E33

WhereRl, kr,l, ul(l=0,1,2)and ui are reflection constants, the wave vector, the particle displacement vectors of reflection P-, SV- and SH-waves and the incident wave respectively, which are given by the solutions of the Christoffel equation for the isotropic solid. Whenr=0, we have

[u(0)]=[L][R0R1R2]+[ui(0)],E34
L]=[x^u0x^u1x^u2y^u0y^u1y^u2z^u0z^u1z^u2],E35
[u(x)]=[x^u(x)y^u(x)z^u(x)]TE36

where the superscript T denotes transpose and i^(i=x,y,z) is the unit vector of the i-direction. Derivative of (33) with respect to x is

xu=l(jkr,lx^)Rlulejkr,lr+(jkix^)uiejkir,E37

and for r=0 we have

[R0R1R2]=[K]1(jkix^[x^ui(0)y^ui(0)z^ui(0)]+[x^ux|x=0y^ux|x=0z^ux|x=0]).E38

Here

[K]=j[kr,0x^x^u0kr,1x^x^u1kr,2x^x^u2kr,0x^y^u0kr,1x^y^u1kr,2x^y^u2kr,0x^z^u0kr,1x^z^u1kr,2x^z^u2].E39

Using eqs.(34) and (38), we have

[u(0)]=[L][K]1[x^ux|x=0y^ux|x=0z^ux|x=0]+([I]+jkix^[L][K]1)[x^ui(0)y^ui(0)z^ui(0)].E40

Equation (31) yields

[u]x|x=0=[Ci1]1[Ct1][u]x|x=+0,E41
Ci1]=[λ+2μ000μ000μ],E42
[Ct1]=[λ+2μ000μ000μ],E43
Ct1]=[λ+2μ000sysxμ000szsxμ].E44

Substituting eqs. (30) and (41) into eq. (40), we get the Robin condition:

[L][K]1[Ci1]1[Ct1][u]x|x=+0[s][u(+0)]=([I]+jkix[L][K]1)[ui(0)].E45

After we solve the distributions of the particle displacements in the PML by COMSOL MultiPhysics, the reflection coefficients Rl(l=0,1,2) are computed with the following relation derived from eq. (34):

[R0R1R2]=[L]1([sx000sy000sz][u(+0)][ui(0)]).E46

Here u(+0) and ui(0) are particle displacements at PML’s incident side given by FEA solution and known incident field vector of displacements.

3.1.2. Discretized wave number

The finite element approximation of the propagating elastic fields changes the propagation constant given by the Christoffel equation, which is called the intrinsic wave numberβPML, to discretized wave numberβ˜PML. Table 2 shows discretized wave numbers for nodal finite elements with the polynomial interpolate function as shown in Fig.2 after Scott[22]. Here β and h are the x-component of the intrinsic wave number of P-, SV- or SH-wave propagating in the PML and equal interval between nodes, respectively. Figure 3 shows the difference of the discretized wave number and the intrinsic wave number as the function of the x-axis propagation constant.

Table 2.

Discretized wave number in PML.

Figure 3.

Phase error and attenuation error as a function of /(βh) for 1st-, 2nd-, 3rd-, and 4th order elements.

3.1.3. Transfer matrix analysis

Because the structure shown in Fig.1 is a layered structure where the propagation constants in the isotropic solid and its PML are given as the intrinsic wave numbers and discretized wave numbers respectively, we can compute the reflection coefficient by the transfer matrix.

In this section, we consider the fields composed of P- and SV-waves propagating on the x-y plane with the same y-component ky of the P- and SV-wave numbers only since SH-waves are not coupled with P- or SV-waves and SH-wave scattering problem is straightforward. Assuming that field distributions do not vary in the z-direction, we have the particle displacements in the solid:[15]

ux,m=ejkyy(i=14Ai,mfx,isxejkx,isxx),E47
uy,m=ejkyy(i=14Ai,mfy,iejkx,isxx).E48

Here, kx,iis the x-component of the intrinsic wave number for the isotropic solid and the discretized wave number for the PML, the subscripts i = 1 and i = 3 denote P-waves propagating to +x- and -x-direction respectively, i = 2 and i = 4 denote SV-waves propagating to +x- and -x-direction respectively, and Ai,m(i=1,2,3,4,m=0,1) is the amplitude at x=0 in the isotropic solid (m=0) or PML (m=1). fx,iand fy,i are shown in Table 3. Here θp and θs are angles between the x-direction and the wave vectors of P-waves or SV-waves as shown in Fig. 1. In the isotropic region, we setsx=1.

Table 3.

Displacement directions of P- and SV-waves, fx,i and fy,i.

Using the boundary conditions at x=0 and x=L, and eliminatingAi,1, we get the relation

[A3,0A4,0]=[X31X32X41X42][X11X12X21X22]1[A1,0A2,0]E49

where [X] is the square matrix with four columns and rows given by

[X]=[Y0]1[s][Y1][T(L)][Y1]1[s]1.E50

Here,

[s]=[sx0000sy0000sxsysz00001sz],E51
[Y0]=[Y11,0Y12,0Y13,0Y14,0Y21,0Y22,0Y23,0Y24,0Y31,0Y32,0Y33,0Y34,0Y41,0Y42,0Y43,0Y44,0],E52
Y1i,0=fx,i,E53
Y2i,0=fy,i,E54
Y3i,0=j(kx,i(λ+2μ)fx,i+kyλfy,i),E55
Y4i,0=j(kyμfx,i+kx,iμfy,i),E56
[Y1]=[Y11,1Y12,1Y13,1Y14,1Y21,1Y22,1Y23,1Y24,1Y31,1Y32,1Y33,1Y34,1Y41,1Y42,1Y43,1Y44,1],E57
Y1i,1=fx,i/sx,E58
Y2i,1=fy,i,E59
Y3i,1=j(kx,isx1sx(λ+2μ)fx,isx+ky1sxλfy,i),E60
Y4i,1=j(kysxμfx,isx+kx,isx1sxμfy,i),E61
T(x)]=[ejkx,1sxx0000ejkx,2sxx0000ejkx,3sxx0000ejkx,4sxx].E62

The reflection coefficients at the boundaryx=0, we obtain (A3,0/A1,0 and A4,0/A1,0) with A2,0=0 when the incident wave is the P-wave and in the case of SV-wave incidence we obtain (A3,0/A2,0 and A4,0/A2,0) withA1,0=0.

3.2. Computed results

Figure 4 shows the computed results of the reflection coefficient dependence on 2π/(βh) in the case of the SH-wave incidence with incident angleθ=0, the attenuation coefficient sxI(x)=0.1 and normalized thicknessksL=24π. Here ks is the intrinsic wave number of the SH-wave in the isotropic solid. Decreasing the interval between adjacent nodal pointsh, the reflection coefficient approaches the value estimated by the truncation effect which caused by the reflection at the PML end terminal and can be estimated by attenuated waves in the PML, 20log10(exp(2ksLs2I))=20×4.8πlog10e=131dB. A higher order element causes lower reflection because of a better approximation of the intrinsic wave number. Figure 5 shows dependence on the incident angle. Smaller intervals of finite element nodes, ksh=0.1π, gives a better approximation thanksh=0.2π. Increasing the incident angle, β decreases and the approximation of the intrinsic wave number with discretized wave number becomes better. Hence, the reflection by FE-discretization decreases. However, incident angle becomes larger than the angle such as about 63.5 degrees for 1st order element, reflection increases because decreasing β yields decreasing of wave attenuation in PMLs and the reflection by the truncation effect increases. Figures 4 and 5 show that the results of the transfer matrix agree well those of FEA and we confirm that the reflection of the FE-model of the PML may be explained with the discretized wave number and the truncation effect.

Figure 4.

Dependence of SH-wave perpendicular incidence on 2π/(βh) for ksL=24π and sxI(x)=0.1.

Figure 5.

Dependence of reflection coefficients of the SH-wave perpendicular incidence with the attenuation coefficient sxI(x)=0.1 and normalized thickness ksL=24π on 2π/(βh). Here N is the number of FEs.

Figure 6.

Dependence of reflection coefficients on 2π/(βh) in the case of perpendicular P-wave incident, ksL=24π and sxI(x)=0.1.

Figure 7.

Dependence on P-wave incident angle θi for ksh=0.2π and ksL=24π.

We consider an isotropic solid and its PML with the Poisson ratio σ = 0.3 in this section.

Next, we consider P-wave or SV-wave scattering problems. Figure 6 shows the computed result of the reflection coefficient dependence on 2π/(βh) in case of the P-wave perpendicular incidence with the attenuation coefficient sxI(x)=0.1 and normalized thicknessksL=24π. Because the result of the SV-wave perpendicular incidence is the same result of the SH-vave perpendicular incidence owing to a symmetry of the problem, Fig. 4 also shows the result of SV-wave incidence with the attenuation coefficient sxI(x)=0.1 and normalized thicknessksL=24π. Both cases also approaches the value estimated by the truncation effect, -70.0dB and -131dB. Note that the wave number of the P-wave is μ/(λ+2μ)=2/7 times of the SV-wave wave number. Dependencies of P- and SV-wave reflections on P- and SV-wave incident angle are shown in Figs. 7 and 8, respectively. Reflection coefficients of incident waves computed by the transfer matrix except the range that is larger than the critical angle, about 32.3 degrees, of SV-wave incidence are good agreement with the results of FEA. However, reflection coefficients of converted waves from incident waves are smaller for the SV-wave excited by the incident P-wave and larger for the P-wave than the results of FEA. We still can not explain this discrepancy.

Figure 8.

Dependence on SV-wave incident angle θi for ksh=0.2π and ksL=24π.

Increasing sxI, the reflection coefficient of P-wave in case of the P-wave incidence decreases in the incident angle range that is larger than 59 degrees. In the lower range, the reflection does not decrease because FE-discretization effect dominates the reflection. In the case of SV-wave incidence, we confirm that the P-wave converted from the SV-wave is amplified in the incident angle range that is lager than the critical angle when P-wave’s wave number is zero in the isotropic region because of PML’s intrinsic characteristics for non-propagating waves.

Advertisement

4. Conclusions

In this chapter, first, PMLs in the Cartesian, the cylindrical and the spherical coordinates for elastic waves in solids were derived from differential forms on manifolds. Our results show that PML parameters in any orthogonal coordinate system for elastic waves in solids may be determined by the same procedure in the Cartesian coordinates. Next, scattering of elastic waves in an isotropic solid was analyzed by field analysis in the thick layer in the one dimension. Numerical results show that the reflection from PMLs by the transfer matrix of elastic waves approximates the numerical results of FE-models successfully. We concluded that the reflection by FE discritization may be explained by FE-approximation of the intrinsic wave number.

References

  1. 1. J.P. Berenger1994A perfectly matched layer for the absorption of electromagnetic wavesJ. Comput. Phys. : 185200
  2. 2. W.C.Chew and W.H. Weedon1994A 3D perfectly matched medium for modified Maxwell’s equations with stretched coordinates.Microwave Opt. Tech. Lett. : 599604
  3. 3. F.L.Teixeira and W.C.Chew2000Complex space approach to perfectly matched layers: a review and some new developmentsInt. J. Numer. Model. : 441455
  4. 4. SacksZ. S.KingslandD. M.R.Lee and J.F.Lee (1995A perfectly matched anisotropic absorber for use as an absorbing boundary conditionIEEE Trans. Antennas Propag. : 14601463
  5. 5. F.D.Hastings, J.B.Schneider and S.L.Broschat1996Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagationJ Acoust. Soc. Am. : 30613069
  6. 6. W.C.Chew and Q.H.Liu1996Perfectly matched layers for elastodynamics: a new absorbing bounary condition. J. Comput. Acoust. : 341359
  7. 7. Q.H.Liu1999Perfectly matched layers for elastic waves in cylindrical and spherical coordinatesJ. Acoust. Soc. Am. : 20752084
  8. 8. Y.Zheng and X.Huang2002Anisotropic perfectly matched layers for elastic waves in cartesian and curvilinear coordinates. MIT Earth Resources Laboratory Industry Consortium Meeting.
  9. 9. S.Ballandras,D.Gachon, J.Masson and W.Daniau2007Development of absorbing conditions for the analysis of finite dimension elastic wave-guides. Proc. Int. Freq. Contr. Symp.: 729732
  10. 10. MayerM.ZaglmayrS.WagnerK.SchöberlJ.2007Perfectly matched layer finite element simulation of parasitic acoustic wave radiation in microacoustic devices. Proc. IEEE Ultrason. Symp.: 702706
  11. 11. LiY.MatarO. B.PreobrazhenskyV.PernodP.2008Convolution-Perfectly Matched Layer(C-PML) absorbing boundary condition for wave propagation in piezoelectric solid.Proc. IEEE Ultrason. Symp.: 15681571
  12. 12. U.Basu and A.K.Chopra2003Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementationComput. Methods Appl. Mech. Eng. : 13371375
  13. 13. C.Michler, L.Demkowicz, J.Kurtz and D.Pardo2007Improving the performance of perfectly matched layers by means of hp-adaptivityNum. Methods Partial Diffe. Equa. : 832858
  14. 14. A.Taflove and S.C.Hagness2005Computational electrodynamics(Artech House, Boston) 3rd ed. Ch7, 273
  15. 15. T.Shimada and K.Hasegawa:2010Perfectly matched layers for elastic waves propagating in anisotropic solids. IEICE Trans. : 215223in Japanese].
  16. 16. ShimadaT.HasegawaK.2010Perfectly matched layers in the cylindrical and spherical coordinates. Jpn. J. Appl. Phys. : 07HB08.
  17. 17. T.Shimada, K.Hasegawa, and S.Sato2009Absorbing characteristics of electromagnetic plane waves in perfectly matched layers discretized by finite element method. Keisan Suri Kogaku Ronbunshu : 04091211in Japanese].
  18. 18. W.C.Chew and J.M.Jin1996Perfectly matched layer in the discretized space: an analysis and optimaization. Electromagnetics : 325340
  19. 19. F.Collino and B.P.Monk:1998Optimizing the perfectly matched layerComput. Methods. Appl. Mech. Eng. : 157171
  20. 20. BerumúdezA.Hervella-NietoL.A.prietoRodrïguezR.2007An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems. J. Compu. Phys. : 469488
  21. 21. A.Berumúdez, L.Hervella-Nieto, A.prieto, and R.Rodrïguez2010Perfectly matched layers for time-haromonic second order elliptic problems. Arch. Comput. Methods Eng. : 77107
  22. 22. ScottW. R.Jr 1994Errors due to spatial discretization and numerical precision in the finite-element methodIEEE Trans. Antennas Propagat. : 15651570

Written By

Koji Hasegawa and Takao Shimada

Published: 10 October 2012