Components of a stress tensor in a PML material of an isotropic solid in the Cartesian coordinate.
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. 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. 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 and by using analytic continuation in anisotropic solids in 2002. Recently validity and usefulness of PMLs derived from the analytic continuation in piezoelectric solids was demonstrated. [9, 10, 11] Hastings et. al. 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 about explanations and finite element (FE) implementation of PMLs for time-harmonic elastodynamics, Michler et al. about derivation of material constants of PML in FE method by analytic continuation, and Taflove and Hagness 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. In addition, the rule for determining PML parameters in the Cartesian coordinate holds in the cylindrical and spherical coordinates.
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. 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. For FD-TD method Collino and Monk also carried out such an investigation. 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 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. 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.
2. Derivation of perfectly matched layers for elastic waves by using complex coordinate stretching and differential form
2.1. Differential form
A particle displacement vector, particle velocity vector, density of momenta, stress tensor and displacement gradient tensor are given as follows:
where and are the contravariant and covariant basis vectors, and represent the tensor product and the cross product, respectively. 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, , the relation of tensor components is
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 is, or 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, , we have
Here, the superscript c denotes the value in the complex coordinate and the mass density and the stiffness 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 and of the general orthogonal coordinate system are dual to the unit vectors and, we have
Here with and being scale factors of general orthogonal coordinate systems and, respectively. Note that the scale factors are given by follows: in the cylindrical coordinate, and in the spherical coordinate,,. In addition, in the Cartesian coordinate,.
The quotient rule and eqs. (9)- (14) yield PML material constants: the mass density is
and the stiffness is
Here, , , in the cylindrical coordinate system with its complex coordinate, and, , in the spherical coordinate system with its complex coordinate. In addition, 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.
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 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.
First we consider Newton’s equation of motion. In a cylindrical coordinate, the governing equations are
Here, we use phasor notation. The time dependences of the fields are where is angular frequency. Applying CCS with a complex coordinate, multiplying the CCS equations by and using the assumption of, , and, we get the governing equations in the PML:
Next we consider the displacement gradient in the complex coordinate. Using definition and the assumption of, , , , and applying CCS with a complex coordinate, we have the relation:
Hence we have.
Using the quotient rule and the constitutive equation, , we get the constitutive equation of the PML in the real coordinate. Therefore, we may define the stiffness of the PML derived by the analytic continuation:
2.4. Comparison with PML material constants derived from differential forms and the analytic continuation
By the analytic continuation, Zheng and Huang derived the mass density and stiffness of PML in the cylindrical and spherical coordinates: and. The mass density agree with our result, eq. (15), because multiplying the stress tensors by the Jacobian of the coordinate transformation, , 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 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 and denote the -, -, - and -axis, respectively, and is the Kronecker delta. Components of stiffness tensors derived from the differential form and analytic continuation, and, respectively, are given by
Table 1 shows all components of the stress tensor computed with and. gives and we predict that rotational forces may be observed. With, however, we have a symmetric stress tensor,.
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, and are propagation angles of P-waves and SV- or SH-waves, is thickness of the PML, and 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, where is the imaginary unit and is the angular frequency.
When the stiffness component of the isotropic solid is given by where and are the Lamé constants and is the Kronecker delta, the stiffness component of its PML is
Here is a coordinate stretching factor of -direction. The mass density of the PML is given by
Here is the mass density of the isotropic solid. For examining absorbing performance of PMLs in the direction, taking an assumption of considering fields being consisted by plane waves propagating on the plane, we have a differential equation in one variable: from Newton’s equation of motion and constitutive equation we get the differential equation in the PML
where is the component of the particle displacement in the -direction ().
In this case, we may choose the coordinate stretching factor as follows:
Here is the imaginary part of and therefore a real function, which controls absorbing performance of propagating waves in PMLs.
Boundary conditions at the interface of isotropic solid and PML, , are the nonslip condition and the continuous condition of the normal component of the stress:
At the terminal of PML, , the boundary condition is
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 and a formula for reflection coefficients.
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:
Where, , and 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, we have
where the superscript denotes transpose and is the unit vector of the -direction. Derivative of (33) with respect to is
and for we have
Using eqs.(34) and (38), we have
Equation (31) yields
After we solve the distributions of the particle displacements in the PML by COMSOL MultiPhysics, the reflection coefficients are computed with the following relation derived from eq. (34):
Here and 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, to discretized wave number. Table 2 shows discretized wave numbers for nodal finite elements with the polynomial interpolate function as shown in Fig.2 after Scott. Here and are the -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 -axis propagation constant.
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:
Here, is 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 is the amplitude at x=0 in the isotropic solid (m=0) or PML (m=1). and are shown in Table 3. Here and 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 set.
Using the boundary conditions at x=0 and x=L, and eliminating, we get the relation
where [X] is the square matrix with four columns and rows given by
The reflection coefficients at the boundary, we obtain and with when the incident wave is the P-wave and in the case of SV-wave incidence we obtain and with.
3.2. Computed results
Figure 4 shows the computed results of the reflection coefficient dependence on in the case of the SH-wave incidence with incident angle, the attenuation coefficient and normalized thickness. Here is the intrinsic wave number of the SH-wave in the isotropic solid. Decreasing the interval between adjacent nodal points, 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, dB. 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, , gives a better approximation than. 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.
Next, we consider P-wave or SV-wave scattering problems. Figure 6 shows the computed result of the reflection coefficient dependence on in case of the P-wave perpendicular incidence with the attenuation coefficient and normalized thickness. 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 and normalized thickness. 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 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.
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.
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.