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

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.

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 P MLs,fr omt h eP MLsidea r esu ppr essedb elow-4 5a n d-8 0dBw it h8a n d1 6gridspa c esof 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.

Differential form
A particle displacement vector u, particle velocity vector v, density of momenta P,s t r e s s tensor T and displacement gradient tensor F are given as follows: where ∂ ∂x i and dx i (i = 0, 1, 2) are the contravariant and covariant basis vectors, ⊗ and ∧ represent the tensor product and the cross product, respectively.d is the exterior differential operator.Newton's equation of motion is 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 = 181 Perfectly Matched Layer for Finite Element Analysis of Elastic Waves in Solids Using CCS [2,3,8] given by X i = x i si (τ)dτ = x i siR (τ)+js iI (τ)dτ with the two real functions siR (τ) and siI (τ) , we have the relation Here j is the imaginary unit.

PMLs in the Cartesian, the cylindrical and the spherical coordinates
In the complex coordinate stretching (CCS), we consider that the real coordinate (x 0 , x 1 , x 2 ) 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, (X 0 , X 1 , X 2 ),wehave Here, the superscript c denotes the value in the complex coordinate and the mass density ρ and the stiffness C ijkl (i, j, k, l = X 0 , X 1 , X 2 ) 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) Here i si with h r i and h c i being scale factors of general orthogonal coordinate systems (x 0 , x 1 , x 2 ) and (X 0 , X 1 , X 2 ), respectively.Note that the scale factors h i are given by follows: in the cylindrical coordinate (r, θ, z) h 0 = 1, h 1 = r, h 2 = 1 , and in the spherical coordinate (r, θ, φ) h 0 = 1, h 1 = r, h 2 = r sin θ.In addition, in the Cartesian coordinate, The quotient rule and eqs.( 9)-( 14) yield PML material constants: the mass density ρ PML is and the stiffness is Here, r sin θ in the spherical coordinate system (r, θ, φ) with its complex coordinate (R, Θ, Φ).In addition, s i = si in 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.

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) ,t h e governing equations are 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 s 0 s 1 s 2 and using the assumption of u c i = u i , R = r, Θ = θ and Ẑ = ẑ,wegetthe governing equations in the PML: Here, the mass density of PML is defined as ρ PMLA = s 0 s 1 s 2 ρ.When we rewrite components of stress tensors in the PMLs as T PMLA ij = s 0 s 1 s 2 s j T c ij , we may identify eqs.( 20)∼( 22) to eqs.( 17)∼ (19).
Next we consider the displacement gradient ∇u c =[Γ c kl ] in the complex coordinate (R, Θ, Z).Using definition du = ∇u • ( xi h i dx i ) and the assumption of u c i = u i , R = r, Θ = θ, Ẑ = ẑ, and applying CCS with a complex coordinate (R, Θ, Z), we have the relation: Hence we have Γ c kl = 1 s l Γ kl .Using the quotient rule and the constitutive equation, T c ij = C ijkl Γ c kl , we get the constitutive equation of the PML in the real coordinate (r, θ, z): T PMLA ij = s 0 s 1 s 2 s j s l C ijkl Γ kl .Therefore, we may define the stiffness of the PML derived by the analytic continuation:

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 = s 0 s 1 s 2 ρ and C PMLA ijkl = s 0 s 1 s 2 s j s l C ijkl .The mass density agree with our result, eq. ( 15), because multiplying the stress tensors by the Jacobian of the coordinate transformation, s 0 s 1 s 2 , 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 T PMLA ij = s 0 s 1 s 2 s j T c ij and use the assumption To show a difference between PML material constants, we consider an isotropic solid with following stiffness constants in the Cartesian coordinate Here, λ and μ are the Lamé constants of an isotropic solid, the subscripts i, j, k and l denote the x i -, x j -,x k -a n dx l -axis, respectively, and δ ij is the Kronecker delta.Components of stiffness tensors derived from the differential form and analytic continuation, C PML ijkl and C PMLA ijkl , respectively, are given by Components of a stress tensor in a PML material of an isotropic solid in the Cartesian coordinate.
Table 1 shows all components of the stress tensor computed with C PML ijkl and C PMLA ijkl .C PMLA ijkl gives Tij = Tji (i = j) and we predict that rotational forces may be observed.With C PML ijkl , however, we have a symmetric stress tensor, T ij = T ji (i = j).

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, θ p and θ s are propagation angles of P-waves and SV-or SH-waves, L is thickness of the PML, k i and k r,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 are exp(jωt), wherejistheimaginaryunitandω is the angular frequency.
When the stiffness component of the isotropic solid C ijkl (i, j, k, l = x, y, z) is given by C ijkl = λδ ij δ kl + μ(δ ik δ jl + δ jl δ ik ) where λ and μ are the Lamé constants and δ ij is the Kronecker delta, the stiffness component of its PML C PML ijkl is  Here s i (i = x, y, z) is a coordinate stretching factor of i-direction.[2] The mass density of the PML ρ PML is given by 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 x-y plane, we have a differential equation in one variable x:f r o m Newton's equation of motion and constitutive equation we get the differential equation in the PML where u i 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: Here s xI (x) is the imaginary part of s x 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] u i (−0)=s i u i (+0), (30) At the terminal of PML, x = L, the boundary condition is

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.In the half isotropic solid, the field distribution, h 12 nn + 1 components of the particle displacement and stress, can be expressed by superposition of incident and reflected plane waves: Where R l , k r,l , u l (l = 0, 1, 2) and u i 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.When r = 0,wehave where the superscript T denotes transpose and î(i = x, y, z) is the unit vector of the i-direction.
Derivative of (33) with respect to x is and for r = 0 we have 187 Perfectly Matched Layer for Finite Element Analysis of Elastic Waves in Solids

Here
(39) Using eqs.( 34) and (38), we have Substituting eqs.( 30) and (41) into eq.( 40), we get the Robin condition: After we solve the distributions of the particle displacements in the PML by COMSOL MultiPhysics, the reflection coefficients R l (l = 0, 1, 2) are computed with the following Order of finite element Discretized wave number βPML Table 2. Discretized wave number in PML.
relation derived from eq. ( 34): Here u(+0) and u i (−0) are particle displacements at PML's incident side given by FEA solution and known incident field vector of displacements.

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.

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 k y o ft h eP -a n dS V -w a v en u m b e r so n l ys i n c eS H -w a v e s 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] u x,m = e −jk y y 4 Here, k x,i is the x-component of the intrinsic wave number for the isotropic solid and the discretized wave number for the PML, the subscripts i = 1a n di = 3 denote P-waves propagating to +x-a n d−x-direction respectively, i = 2a n di = 4d e n o t eS V -w a v e s propagating to +x-a n d−x-direction respectively, and A i,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).f x,i and f y,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 -wavesasshowninFig.1.Intheisotropicregion,wesets x = 1.i 12 34 f x,i cos θ p − sin θ s − cos θ p sin θ s f y,i sin θ p cos θ s sin θ p cos θ s Table 3. Displacement directions of P-and SV-waves, f x,i and f y,i .
Using the boundary conditions at x = 0andx = L, and eliminating A i,1 , we get the relation where [X] is the square matrix with four columns and rows given by Here, e −jk x,1 s x x 000 0e −jk x,2 s x x 00 00 e −jk x,3 s x x 0 000 e −jk The reflection coefficients at the boundary x = 0, we obtain (A 3,0 /A 1,0 and A 4,0 /A 1,0 ) with A 2,0 = 0 when the incident wave is the P-wave and in the case of SV-wave incidence we obtain (A 3,0 /A 2,0 and A 4,0 /A 2,0 ) with A 1,0 = 0.

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 s xI (x)=0.1 and normalized thickness k s L = 24π.H e r ek s is the intrinsic wave number of the SH-wave in the isotropic solid.Decreasing the interval between adjacent nodal points h,t h e 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, 20 log 10 (exp(2k s Ls 2I )) = 20 × 4.8π log 10 e = 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, k s h = 0.1π, gives a better approximation than k s h = 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.
We consider an isotropic solid and its PML with the Poisson ratio σ=0.3 in this section.

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.

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

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

Figure 5 .
Figure 5. Dependence of reflection coefficients of the SH-wave perpendicular incidence with the attenuation coefficient s xI (x)=0.1 and normalized thickness k s L = 24π on 2π/(βh).HereN is the number of FEs.

Figure 7 .
Figure 7. Dependence on P-wave incident angle θ i for k s h = 0.2π and k s L = 24π.

Figure 8 .
Figure 8. Dependence on SV-wave incident angle θ i for k s h = 0.2π and k s L = 24π.