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

## 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.

## 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

where

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,

Using CCS[2, 3, 8] given by

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

Here, the superscript c denotes the value in the complex coordinate and the mass density

Here

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

and the stiffness is

Here,

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

Here, we use phasor notation. The time dependences of the fields are

Here, the mass density of PML is defined as

Next we consider the displacement gradient

Hence we have

Using the quotient rule and the constitutive equation,

### 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:

To show a difference between PML material constants, we consider an isotropic solid with following stiffness constants in the Cartesian coordinate

Table 1 shows all components of the stress tensor computed with

## 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

We use the phasor notation and assume that the time dependences of all fields are

When the stiffness component of the isotropic solid

Here

Here

where

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

Here

Boundary conditions at the interface of isotropic solid and PML,

At the terminal of PML,

### 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

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

where the superscript

and for

Here

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

Equation (31) yields

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

Here

#### 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

#### 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 *k*_{y} 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]

Here, *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 *x*=0 in the isotropic solid (*m*=0) or PML (*m*=1). *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

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

Here,

The reflection coefficients at the boundary

### 3.2. Computed results

Figure 4 shows the computed results of the reflection coefficient dependence on *β* 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

Increasing *s*_{xI}, 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.

## 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.