Perfectly matched layer for second-order time-domain elastic wave equation: formulation and stability
PPerfectly matched layer for second-order time-domainelastic wave equation: formulation and stability
Hisham Assi ∗ , Richard S. C. Cobbold Institute of Biomaterials and Biomedical Engineering, University of Toronto,164 College Street, Toronto, M5S 3G9, Canada
Abstract
A time domain system of equations is proposed to model elastic wave propaga-tion in an unbounded two-dimensional anisotropic solid using perfectly matchedlayer (PML). Starting from a system of first-order frequency domain stress-velocity equations and using complex coordinate stretching approach with a two-parameter stretch function, a second-order formulation is obtained. The finalsystem, which consists of just two second order equations along with four aux-iliary equations, is smaller than existing formulations, thereby simplifying theproblem and reducing the computational cost. The discrete stability of the so-lutions for a given mesh size is examined with the help of a plane-wave analysisof the corresponding continuous problem. It is shown that increasing the scalingparameter of the stretch function leads to significant stability improvements forcertain anisotropic media that have known issues. Numerical computations fordifferent isotropic and anisotropic media are used to illustrate the results.
Keyword : Perfectly matched layers; Elastic waves; Discrete stability; Secondorder time-domain
Contents ∗ Email address: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] D ec Formulation of PML for elastic wave propagation 14
Numerical simulations of wave propagation in an unbounded media need special trunca-tion methods to avoid spurious wave reflections from the computational domain bound-aries. Absorbing boundary conditions (ABCs) [1] were first used. Such conditions workwell when the waves are normally incident as in the case for 1D simulations, but thisapproach has limitations for higher dimensions. A more effective technique, as firstdescribed by B´erenger in 1994 [2], is to terminate the computational domain with aperfectly matched layer (PML). Figure 1 illustrates the use of such a layer consisting ofa hypothetical absorbing material that terminates the computational domain in sucha way that the waves decay exponentially with negligible reflections from the outerboundaries, regardless of the incident angle. This is true for the case of an infinitelyfine mesh i.e, for the continuous limit. In practice, a non-zero mesh element size causessome numerical reflections from the inner boundary of the PML, but these can be madevery small, making PML an efficient means for modeling a variety of wave phenomenasuch as electromagnetic waves, acoustic waves in fluids, and elastic waves in solids.For electromagnetic wave simulations, B´erenger [2] showed that by adding specific2igure 1: Illustrating the use of a perfectly matched layer (PML) for achieving near-perfect modeling of the solution to the unbounded wave radiation problem.conductivity parameters to Maxwell’s equations perfect matching and decaying of thepropagating waves in the PML could be achieved. An alternative method is to assumethat the material contained within the PML is a uniaxial anisotropic media [3–5],generally referred to as the uniaxial PML approach. In this method the original formof the wave equation is retained but with frequency-dependent tensors as the materialproperties which makes it suitable for frequency domain simulations. A third methodwith greater generality and flexibility is the complex coordinate stretching approach36]. In fact, the conductivity parameter introduced by B´erenger [2] can be thought ofas a parameter in a stretch function that extends the spatial coordinate in the layer tothe complex plane. The addition of more parameters was subsequently proposed withthe aim of making the method causal [7]. Although the original PML was subsequentlyfound to be causal [8, 9], other benefits accrued from this new multi-parameter stretchfunction. Specifically, it was found that strong absorption occurred for the evanescentwaves, improved absorption occurred at grazing angles [9–11], and improved stabilitywas achieved in the PML for certain anisotropic elastic media [11–13].Many PML formulations have been introduced for elastic wave propagation [9, 13–18] as well as for general hyperbolic equations [19]. Amongst these the split-field for-mulations usually make use of a single parameter stretch function and are typicallydescribed by systems of first order equations with double the number of physical equa-tions such as those used by B´erenger [2]. Unsplit field formulations use the physicalfields variables along with extra auxiliary variables that are typically needed to obtainthe time-domain equations from the frequency-domain equations. The use of multi-parameter stretch function usually requires a convolution to obtain a time-domain for-mulation, leading to the name convolutional PML [20] for many of the unsplit fieldmodels. The majority of these formulation uses a large number of equations (10 ormore) to describe elastic wave propagation in the PML which affects the computationaltime and resources. Stability is a known issue in PMLs [9, 11–13, 21, 22], especiallyfor some anisotropic solids. Some methods for addressing this problem have been pro-posed [9, 12, 13]. In particular, by controlling the stretch function parameters and themesh size the discrete stability was improved for certain cases where the correspondingcontinuous problems were unstable [12, 13, 22].The purpose of this paper is to introduce second order time domain formulationfor elastic wave propagation in isotropic and anisotropic solids in two space dimen-4ions. Second order equations emerge directly from Newton’s second law which makethem more robust as compared to the first order velocity-stress system of equations[22]. Moreover, the second order equations are more readily implemented in commonnumerical schemes [23], such as those used in PDE software packages like the finiteelement method-based (FEM) COMSOL Multiphysics (COMSOL, Inc., Burlington,Mass., U.S.A.) as used in this work. Other advantages accrue from using this formula-tion. First, it has a smaller number of equations than the classical and convolutionalmodels, thereby simplifying the numerical implementation. Second, it has greater long-time stability for certain anisotropic media that are typically unstable in classical PMLsimulations. A simple method to further improve the discrete stability is proposed. Inthe next section we describe the background needed for obtaining the PML equations.This is followed by the derivation of our second order formulation. Then, with the helpof a plane-wave analysis, the stability analysis is formulated. Numerical results arepresented and discussed for both isotropic and anisotropic media.
The propagation of waves an elastic medium can be described using Newton’s secondlaw, Hook’s law, and the linear approximation of the strain. These lead to the followingthree equations respectively: ρ ∂ u i ∂t = d (cid:88) j =1 ∂σ ij ∂x j , (1) σ ij = d (cid:88) k,l =1 C ijkl ε kl , (2)5 kl = 12 (cid:18) ∂u k ∂x l + ∂u l ∂x k (cid:19) , (3)where u i are the components of particle displacement vector, σ ij , and ε kl are the compo-nents of the symmetric stress and strain tensors respectively, C ijkl are the componentsof the fourth order elasticity tensor with the following symmetries: C ijkl = C ijlk = C jikl ,and C ijkl = C klij , and d is the number of space dimensions which is 2 for this work.The source of energy that excites the elastic medium can either be embedded in theboundary conditions or added as a load victor to (1). The above three equations to-gether with the symmetry properties of the elasticity tensor enable the problem to beexpressed as two second order equations in terms of the displacement vector: ρ ∂ u i ∂t = (cid:88) j =1 ∂∂x j (cid:32) (cid:88) k,l =1 C ijkl ∂u k ∂x l (cid:33) . (4)Another way to formulate the problem is through a system of first order equationsin term of stress and velocity. These can be obtained using the same equations as usedto obtain (4), leading to ρ ∂v i ∂t = (cid:88) j =1 ∂σ ij ∂x j ∂σ ij ∂t = (cid:88) k,l =1 C ijkl ∂v k ∂x l , (5)where v i = ∂u i /∂t is the velocity vector component. In such a formulation five firstorder equations are needed to describe the problem. Namely, two velocity vector com-ponents, v i and four stress tensor components, σ ij , which are reduced to three due tothe symmetry in the stress tensor ( σ ij = σ ji ).6 .2 Materials properties All media considered in this work are orthotropic, which is a special case of an anisotropicmedia whose axes of symmetry coincide with x and x . For such a medium the elastic-ity tensor has only four independent components. For simplicity and consistency withthe notation commonly used [24], we replace indices 11 →
1, 22 →
2, 12 →
3, and21 →
3, so that the Hooks law for orthotropic media becomes σ σ σ σ = C C C C C C C C ∂u ∂x ∂u ∂x ∂u ∂x ∂u ∂x . (6)The elasticity coefficients are displayed in this notation in table Table 1.For the purpose of validation, testing, and stability analysis, we chose five mediawhose characteristics are shown in Table 1. Material I is isotropic ( C = C = C + 2 C ) while the others are the anisotropic materials. In particular, media II, III,IV are identical to media II, III, IV as specified by B´ecache et al. [25], and media V,which was also studied in [9, 12], corresponds to zinc crystal. The isotropic mediumwas used to test our PML and, by comparison with theoretical predictions, to validatethe results of our numerical simulations. The anisotropic media was mainly used tostudy the stability.Table 1: Elasticity coefficients for the materials examined.Material C C C C I 7.8 7.8 2 3.8II 20 20 2 3.8III 4 20 2 7.5IV 10 20 6 2.5V 16.5 6.2 3.96 57 .3 Plane waves and slowness curves
To better understanding the wave propagation properties for equation like (4), it isuseful to consider plane wave solutions of the form u = u e i ( k · x − ωt ) , (7)where u ∈ C is the polarization vector, or the amplitude of the wave with wavevector k ∈ R and angular frequency ω ∈ C , and i = √−
1. The dispersion relation between k and ω , can be obtained by substituting (7) into (4). Assuming ρ = 1 and that C ijkl are constants this results in a fourth order polynomial given by F ( ω, k ) = det (cid:32) ω δ ik − (cid:88) j,l =1 C ijkl k j k l (cid:33) = 0 , (8)where δ ik is the Kronecker delta function. For an orthotropic medium this can bewritten as F ( ω, k , k ) = ω − ω (cid:2) ( C + C ) k + ( C + C ) k (cid:3) + C C k + C C k + (cid:0) C C − c − C C (cid:1) k k = 0 , (9)which is the characteristic polynomial of (4) for the orthotropic case. We will refer forthe four roots of (9), ω n ( k , k ) where ( n = 1 ...
4) as the physical modes.Consider the following two conditions on the elasticity tensor C > , C > , C > , and C C > C C (cid:54) = C and C (cid:54) = C . (10)If the first condition is satisfied then the four roots of (10) are all real. Moreover, if thesecond condition is also satisfied then the four roots will be distinct enabling the group8elocity to be defined by V g = ∇ k ω = − ∇ k F ( ω, k ) ∂F ( ω, k ) /∂ω (11)which specifies the direction of energy transport. The slowness vector defined by S = k /ω provides a convenient means for understanding the dispersion relations. Since (10)is homogeneous in k and ω , it can be expressed as F (1 , S , S ) = 0 . (12)For the materials in Table 1 slowness curves, which are the plot of (12), are shownin Figure 2. The inner curve corresponds to the fast wave (the longitudinal or quasi-longitudinal) and the outer curve corresponds to slow waves (shear or quasi-shear). Thephase velocity, V = ω/ | k | = ± / | S | , in each propagation direction can be obtainedfrom the slowness curves. In this work, the maximum and minimum phase velocity fora given material will be referred to as c max , and c min respectively. In addition, followingfrom (11), the direction of the group velocity is normal to these curves. B´ecache et al. [25] found that the stability of the split-field classical PML depends on the shape ofslowness curves and they called this the geometrical stability condition. A perfectly matched layer can be constructed by the analytic continuation of the spa-tial coordinate to the complex domain inside the PML region [6, 26, 27]. Assumingthat the region sufficiently far from that containing the sources and inhomogenities(see Figure 1) is linear and homogeneous, the radiation solution can be written as asuperposition of harmonic plane waves [26]. Because these waves are analytic functions9 s I − II − s III − s s IV Group velocity Phase velocity Slow wave Fast wave 0.5 0 0.5 s V Figure 2: Slowness curves for all the materials whose properties are given in Table 1.Phase and group velocity are indicated for selected points.of the space coordinate, the radiation solutions are also analytic and are subject toanalytic continuation [18, 26, 27].A coordinate transformation x j → ˜ x j ( x j ) : R → C is performed where ˜ x j ( x j ) hasthe value of x j inside the physical domain and is continuous everywhere. Since homo-geneity was assumed close to and inside the PML region, x j appears in the differentialequations only as a partial derivative. Thus, the original wave equation in x j can betransformed into a one in ˜ x j merely by replacing 1 / ∂x j by 1 / ∂ ˜ x j . This transformedequation has the same solution in the physical domain as the original equation, butwithin the PML, it can be made an exponentially decaying solution with no reflectionsat the interface. Unfortunately, solving this differential equation along contours in thecomplex plane can be challenging. This can be avoided by transforming the complex10oordinate back to the real coordinate x j [26].Within the PML the spatial coordinate in the PDEs only appears in the form ofspatial partial derivatives. As a result, instead of defining the transformation x → ˜ x ,the relation between ∂ ˜ x j and ∂x j suffices for the transformation. If the complex stretchfunction is defined as their ratio, i.e., s j ( x j ) = ∂ ˜ x j ( x j ) / ∂x j , then ∂∂ ˜ x j = 1 s j ( x j ) ∂∂x j . (13)Since the stretch function is a complex function in x j , it can be expressed in the two-parameters form: s j ( x j , ω ) = α j ( x j ) (cid:20) i β j ( x j ) ω (cid:21) , (14)where the damping coefficient, β j ≥
0, is responsible for damping the propagatingwave inside the PML. Moreover, the scaling coefficient, α j >
0, is responsible for eitherstretching ( α j >
1) or compressing (0 < α < ω , was added to make the damping wavevector independent. In the physical domain(see Figure 1) ˜ x j ( x j ) = x j , so that β j = 0 and α j = 1, whereas in the PML, β j > α j can differ from 1.To illustrate the effect of the complex coordinate stretching, consider the simple caseof the 1D oscillatory solution shown in Figure 3 (A). Figure 3(B) shows the wave for β ( x ) > α ( x ) = 1 throughout. It can be seen that an exponentiallydamped wave given is present in the PML. Figure 3(C) shows the cases for α ( x ) > β ( x ) > α ( x ) resultingin an apparent increase in the number of cycles, which is equivalent to increasing thespatial frequency, k , in the original coordinate. As subsequently shown, this conceptcan be used to improve the discrete stability. The damping also increased in (C), sincethe coordinate stretching makes the wave travels more and hence, decays more. If the11 A m p li t ude A −
101 x=x B −
101 Spatial variable, x A m p li t ude C
01 Spatial variable, x D Figure 3: Illustrating the effect of complex coordinate stretching for a 1D plane wave,shown in (A), propagating into a PML. The point x = x marks the beginningof the PML; the shaded region in (B), (C), and (D). (B) Shows the case wherethe damping coefficient β (cid:62)
0, and the scaling coefficient α = 1. In (C) thesame value of β as in (B) was used while α >
1. (D) Illustrating the case ofan evanescent wave with α > α ( x ) and β ( x ) are reversed. Thus, if α ( x ) > α j ( x j ) and β j ( x j ). Despite the absence of a rigorous methodology for their choice [17, 18], polyno-12ial functions are often used. For the scaling coefficient, this can be expressed as α j ( x j ) = | x j | < x α j − (cid:16) | x j |− x d (cid:17) m if x ≤ | x j | ≤ x + d, (15)and for the damping coefficient β j ( x j ) = | x j | < x ˜ β j (cid:16) | x j |− x d (cid:17) n if x ≤ | x j | ≤ x + d, (16)where d is the thickness of the PML, 2 x is the dimension of the physical domain, whichis a square centered at the origin as shown in Figure 1, m and n are the polynomialorders, and ˜ α j and ˜ β j are constants that represent the maximum values of α and β respectively. The value of ˜ β j can be expressed in terms of the desired amplitudereflection coefficient ( R j ) due to the reflection from the outer boundary of the PML.For normal incidence, and assuming α j = 1, it can be shown that˜ β j = c max ( n + 1)2 d ln (cid:18) R j (cid:19) , (17)where c max is the highest wave speed which in the case of an isotropic solid, is the longi-tudinal wave speed. The choice of ˜ α j in (15) depends on the desired scaling (stretchingor compression) of the original coordinate. The scaling of the original coordinate issimply the derivative of the real part of ˜ x j with respect to x j , which is equal to α j ( x j ).Hence, the value of ˜ α j is simply the maximum scaling of the original coordinate in thej th direction. The orders of the polynomial functions, m and n , in (15) and (16) cantheoretically be any integer, or even zero. Linear and quadratic polynomials are usuallyused, and will be used in this work unless mentioned otherwise.13hen α j is set equal to unity in (14), the stretch function simplifies s j ( x j , ω ) =1+[ iβ j ( x j ) / ω ], which is the classical stretch function. Another form of the stretch func-tion was introduced by Kuzuoglu and Mittra [7] who added a frequency–shift parameter γ ( x ), such that s j ( x j , ω ) = α j ( x j ) + [ β j ( x j ) / ( γ j ( x j ) − iω )], leading to a PML formula-tions that are usually called convolutional frequency shift (CFS-PML). We chose to usea two-parameter stretch function as described in (14). Besides terminating the evanes-cent waves, other advantages accrue from making α ( x ) (cid:54) = 1. As will be shown, it can beused to improve the stability in the PML. Moreover, the choice of α ( x ) > s j ( x j , ω ) , α j ( x j ) , β j ( x j )will not be used in the remainder of this work. All other coefficients of the PDEs areassumed to be space-dependent only. With the help of the above background, our time-domain PML formulation can beintroduced for the wave propagation in unbounded solids. The derivation starts fromthe first order velocity-stress equations in the frequency domain and concludes with asecond order PML time domain equations in term of the velocity field.
Because the stretch function s j ( x j , ω ) is a function of frequency, the PML formulationwhich uses complex coordinate stretching starts in the frequency domain, and then, ifneeded, the time domain formulation can be obtained by using the inverse Fourier trans-form. The frequency-domain PML equations can be obtained from Fourier transformsof (4) by replacing x by ˜ x , followed by the use of (13) to transform the coordinates-14tretched equations back to the original coordinates yielding: − ω ˆ u i s s ρ = (cid:88) j =1 ∂∂x j (cid:32) (cid:88) k,l =1 s s C ijkl s j s l ∂ ˆ u k ∂x l (cid:33) . (18)In this expression it should be noted that inside the physical domain where s = s = 1,(18) reduces to the frequency domain form of(4). In the PML region, (18) can belooked at as the original equation but with a fictitious medium whose density is s s ρ and whose elasticity tensor is s s C ijkl /s j s l . Both of these coefficients are now complexand frequency dependent.To obtain the velocity-stress formulation for the PML, we proceed in a similarmanner to that used to obtain (18) leading to: − iω ˆ v i ρs s = (cid:88) j =1 s s s j ∂ ˆ σ ij ∂x j − iω ˆ σ ij = (cid:88) k,l =1 C ijkl s l ∂ ˆ v k ∂x l . (19)Either (18) or (19) can be be used for frequency domain simulation. In addition, (19)will be used in the next section to obtain the time domain equations. We proceed by first splitting each of the stress field components in (19) into two non-physical components, σ ij and σ ij , while keeping the velocity field components unsplit.Since s s /s j in (19) does not depend on x j it can be placed inside the x j derivative,15eading to − iω ˆ v i ρ s s = (cid:88) j =1 ∂∂x j (cid:32) s s s j (cid:88) l =1 ˆ σ lij (cid:33) − iω ˆ σ lij = (cid:88) k =1 C ijkl s l ∂ ˆ v k ∂x l . (20)Multiplying the first by − iω , the second by s l , and expanding s and s using (14),results in ρ (cid:2) ( − iω ) + ( − iω ) ( β + β ) + β β (cid:3) ˆ v i = (cid:88) j =1 α j ∂∂x j (cid:34)(cid:18) − iω + β β β j (cid:19) (cid:88) l =1 ˆ σ lij (cid:35) ( − iω ) ˆ σ lij + β l ˆ σ lij = (cid:88) k =1 C ijkl α l ∂ ˆ v k ∂x l . (21)The time domain form of (21) can now be obtained by taking its inverse Fouriertransform without a need for convolution ( − iω ⇒ ∂/∂t ), leading to ρ (cid:20) ∂ v i ∂t + ( β + β ) ∂v i ∂t + β β v i (cid:21) = (cid:88) j =1 α j ∂∂x j (cid:34) (cid:88) l =1 (cid:32) ∂σ lij ∂t + β β β j σ lij (cid:33)(cid:35) ∂σ lij ∂t + β l σ lij = (cid:88) k =1 C ijkl α l ∂v k ∂x l . (22)By substituting ∂σ lij /∂t from the second to the first and simplifying, yields ρ (cid:20) ∂ v i ∂t + ( β + β ) ∂v i ∂t + β β v i (cid:21) = (cid:88) j =1 α j ∂∂x j (cid:34) (cid:88) k,l =1 C ijkl α l ∂v k ∂x l + (cid:88) l =1 (cid:18) β β β j − β l (cid:19) σ lij (cid:35) ∂σ lij ∂t + β l σ lij = (cid:88) k =1 C ijkl α l ∂v k ∂x l . (23)Noting that if j (cid:54) = l , then ( β β / β j ) − β l = β l − β l = 0 so that only four of theeight split stress components (cid:0) σ lij (cid:1) , namely σ jij remain in the first equation. Thesefour non-physical split stress components are needed to solve for the velocity field and16ill be considered as auxiliary variables denoted by σ jij ≡ A ij . Thus, our time domainPML formulation consists of two second-order velocity field equations and four auxiliaryequations that can be expressed as˜ ρ (cid:18) ∂ v i ∂t + b ∂v i ∂t + c v i (cid:19) = (cid:88) j =1 ∂∂x j (cid:34)(cid:32) (cid:88) k,l =1 ˜ C ijkl ∂v k ∂x l (cid:33) + a j A ij (cid:35) ∂A ij ∂t + β j A ij = (cid:88) k =1 C ijkj α j ∂v k ∂x j , (24)where ˜ ρ = α α ρ , ˜ C ijkl = α α C ijkl / α j α l , a j = ( α α / α j ) [( β β / β j ) − β j ], b = β + β , and c = β β . It should be noted that the number of equations in (24) is lessthan that present in the classical form and the convolutional form (typically 10 and13 equations respectively [23]). Other time domain PML formulations follow a similarpattern.If preferred, a set of displacement time domain PML equations can readily be ob-tained by integrating (24) with respect to time. Since the coefficients of (24) are timeindependent and A ij is only an auxiliary variable, this results in equations of the sameform as the above equations but with the velocity field, v i replaced by the displacementfield, u i . It should be noted that, in the physical domain, the two equations in (24)are decoupled, and the displacement form of the first one is identical to the originalequation, (4), which should be the case for any valid PML formulation. In general, when m, n (cid:54) = 0 in (15) and(16), (24) is a variable coefficient PDE in thePML. However, to study the stability of the variable coefficient problem, it is helpfulto assume constant coefficients, which allows use of the plane wave analysis approach[12, 13, 25]. In the physical domain, we know that the roots of the characteristic17olynomial are real and there is no stability issue, but in the PML complex roots canbe present leading to potential instability. When ω is complex, the plane wave solution,as given by (7), becomes u = u e (cid:61){ ω } t e i ( k · x −(cid:60){ ω } t ) . Thus, the sign of the imaginarypart of ω determines the stability of (24). Specifically, if (cid:61) { ω } >
0, the solution growsexponentially with time, alternatively if (cid:61) { ω ( k ) } (cid:54) , ∀ k ∈ R , (25)then (24) is stable.Numerical results and studies [13, 25, 31] have shown that instability starts in in oneor both directions of the PML, but not in the corner region where the full PML equationis involved. Just one direction for the stability analysis will be considered, namely, the x direction, where β = 0 and α = 1. For this case the 8 th order characteristicpolynomial of (24) is F ( ω, k , k , β , α ) ≡ F (cid:20) ( ω + iβ ) ω, k α ω, k ( ω + iβ ) (cid:21) = 0 . (26)where F is defined by (10). Assuming ω = iη , which makes (26) a real-coefficient 8 th order polynomial in η and, according to the complex conjugate root theorem its roots, η ( k , k , β , α ), come in complex conjugate pairs. Hence Lemma 1. , The roots of (26) , ω ( k , k , β , α ) , come in pairs: each pair has the sameimaginary part and the real parts differ only in sign. If none of the four pair of roots of (26) has a positive imaginary part, stability inthe x direction of the PML is assured.First, consider the case in which α = 1 . For this case (26) is identical to theequation for ˜ F pml as given by B´ecache et al [25] as part of the dispersion relation of the18lassical split-field PML (see their equation (64)). Using the perturbation techniques,they studied the stability of F ( ω, k , k , β ,
1) and found, among other results, thefollowing:1. All the necessary and sufficient stability conditions could be expressed in termsof the elasticity coefficients.2. High frequency stability geometric condition (Theorem 2 of their work):It is necessary that all points on the slowness curve satisfy S j × ( V g ) j (cid:62) , (27)for the PML in the x j direction to be stable. This means that the j th componentof the group velocity is in the same direction as the j th component of the slownessvector, which can be readily identified on the slowness curves shown in Figure 2.Violating this condition usually causes the most severe instability. The geometricstability was also found to be necessary condition for other PML formulations[12, 13].3. Because of the symmetries in the orthotropic media, it is enough to consider thefirst quarter of the k − space ( k > k > α (cid:54) = 1. By inspection, it is evident that the roots of F ( ω, α k , k , β , α ) = 0 are the same as the roots of F ( ω, k , k , β Corollary 1.
Changing the scaling parameter, α from unity will cause any root of F =0 to be moved in k − space. For the continuous case, such a movement can never causeany unstable roots to become stable. Therefore, the necessary and sufficient condition or the stability of our constant coefficient continuous problem, as defined by (24) , areexactly the same as the ones reported by B´ecache et al [25] for their split-field system. k × (m ) k × ( m ) A × (m ) B I m ag i na r y pa r t o f × ( s ) Figure 4: Illustrating the effect of incorporation of the scaling coefficient α on the rootsof (26), for the continuous, constant coefficient problem. The color maps showthe imaginary part of the unstable pair of roots, (cid:61) { ω } , for material III. (A)For the classical case of α = 1. (B) For a case in which α = 2. The rootsmerely shifted to higher wavevectors.Figure 4 shows the effect of increasing the scaling parameter on the roots of (26)for material III in Table 1. In (A), the imaginary part of the unstable pair of roots areshown for a range of wavevectors in the first quarter of the k − space, for the case of α = 1. In (B), the same pair of roots is plotted over the same range of wavevectors butfor the case of α = 2. Indeed, as suggested in Corollary 1, the roots were just shifted.Corollary 1 shows that incorporating the scaling parameter will not improve ourcontinuous constant coefficient problem in (24). Though, since PML is meant to beused for numerical simulation, the more relevant question is whether the stability of thediscrete problem that corresponds to an unstable continuous problem can be improved?In fact, this was shown to be case if the unstable continuous modes are not well resoledby the discrete mesh [12, 13, 22], specially for second order formulations [22].If the unstable modes of the PML formulation shown by (24) were in higher wavevec-20or range than can be resolved by the mesh, then our discrete model could be expectedto stable. On the other hand, if the unstable continuous modes are resolvable, in-creasing α shifts the modes to higher wavevectors which might improve the discretestability. This will be the case if the modes of the lower wavevector, which now coverthe resolvable range, have a smaller imaginary part. To investigate this, we return tothe dispersion relation given by (26) and let ξ = k / α (remember α > β , the roots of the dispersion relation are continuous functions in term of ξ ,and k thanks to the implicit function theorem . Noting that decreasing the value of ξ isequivalent to increasing α or decreasing k , as ξ → ω + iβ ) (cid:0) ω − C k (cid:1) (cid:0) ω − C k (cid:1) = 0 , (28)which admits no solution, ω ( k ), with a positive imaginary part. in fact, two of thefour pairs of roots of (28) have the imaginary parts equal to − β , while the imaginaryparts of the other two pairs are equal to zero. Since the root of the dispersion relationare continuous functions in term of ξ = k /α , it follows that: Theorem 1.
By increasing α beyond a certain threshold, the discrete stability of (24) starts to improve. In fact, the results given in subsection 5.2 provide evidence that supports Theorem 1
In all our discrete studies, the source of excitation was a 1 mm diameter infinite cylinderembedded in an infinite 2D medium. To model the infinite medium we assumed aphysical domain of 1.0 cm surrounded by a 1.0 mm PML. The boundary of the cylinderwas assumed to vibrate normally (unless mentioned otherwise) with a velocity, whose21ormalized time-dependence is given by the first derivative of a Gaussian, i.e., v ( t ) = −√ e πf ( t − t ) e − π f ( t − t ) (29)where f is the dominant frequency and t is a source delay time. For all numericalexperiments f = 1500 Hz and t = 1 ms. 90% of the energy of the signal is containedbelow the frequency f c = 1900 Hz.COMSOL Multiphysics was used in combination with MATLAB to numericallysolve (24) using the finite element method. Dirichlet boundary conditions were usedthroughout: specifically, v = v ( t ) ˆn on the surface of the cylinder and v = 0 onthe outer boundary of the computational domain, where ˆn is the normal unit vectorto cylinder surface. A square mesh was used for the PML region, but we retaineda triangular shape in the physical domain. The choice of an appropriate mesh sizeis governed by the shortest wavelength of significance for the propagating pulse, i.e., c min /f c . Since a second order shape function was used in our finite element methodthe mesh size was taken to be h = ( c min /f c ), which corresponds to ten degrees offreedom per wavelength. For time discretization we used an implicit method, specificallythe generalized alpha method. Compared to explicit methods the stability of implicitmethods is not as sensitive to the choice of the time step, time step size of just lessthan h /c max was used, which is sufficient to make optimal use of the mesh. Simulation of wave propagation in unbounded isotropic solid is presented in Figure 5where snapshots of the propagation pulse described by (29) are shown for three instantsof time. To test the accuracy with which these simulations describe the propagatingpulse, we made use of the exact solution for a monochromatic compressional wave22 −
101 Time (ms) N o r m li z ed v e l o c i t y A ↓ B ↓ C ↓ D x (mm) x ( mm ) B − −
55 00.51 x (mm) x ( mm ) C − −
55 00.51 x (mm) D − − − ① ② ③ Figure 5: (B), (C) and (D) are snapshot images showing the amplitude of the particlevelocity for a transient longitudinal wave propagating in the isotropic solidmedium listed in Table 1. The radiation originates from a surface of a 1 mmdiameter cylinder that radial with the velocity profile shown in (A). Markedon the time axis of (A) are the times at which the snapshots in (B), (C) and(D) are taken. Note that (B) and (C) have linear scales, while (D) is in dB’s.The points (cid:192) , (cid:193) , and (cid:194) in (B) are in the physical domain where the solutionsare compared to the analytical solutions in Figure 6.caused by an infinitely long vibrating cylinder in an unbounded isotropic solid [24]. Bymultiplying this with the Fourier transform of (29), then taking the inverse Fouriertransform the time-domain analytical solution was obtained and compared to the FEMresults. As shown in Figure 6 the agreement is excellent, thereby providing good evi-dence for the effectiveness of our PML formulation in simulating unbounded media and23he correctness of the FEM model. − v N o r m li z ed − v N o r m li z ed TheoryFEM − Time (ms) v , v N o r m li z ed − Time (ms) v N o r m li z ed TheoryFEM TheoryFEMTheoryFEM ① ①② ③
A BC D
Figure 6: Validation results: the three points (cid:192) , (cid:193) , and (cid:194) marked on Figure 5 (B)are the locations in the physical domain where the particle velocities wereboth simulated and analytically calculated. The solid line is the theoreticaland the dashed line is from the FEM simulation. (A) and (B) show the twocomponents of the velocity field at point (cid:192) . (C) Velocity field at point (cid:193) . (D)Showing both components of the velocity field at point (cid:194) .Another measure of the effectiveness of the PML can be obtained by looking atthe manner in which the energy in the physical domain evolves in time to ensure thatno energy is reflected back into the physical domain. There are several ways of doingthis [12, 13, 32], one of which is to calculate the maximum magnitude of the particlevelocity in the physical domain (cid:13)(cid:13)(cid:13)(cid:112) v + v (cid:13)(cid:13)(cid:13) ∞ , and to see how this evolves in time. Thisis shown in Figure 7 for the isotropic material as well as for material II, both of whichhave no stability issues. The discrepancies in the energy curve is due to the fact that (cid:13)(cid:13)(cid:13)(cid:112) v + v (cid:13)(cid:13)(cid:13) ∞ is a local measure at the maximum-valued point, and not an averagedmeasure over the whole physical domain like other norms, which on the other hand24
10 2010 − − − Isotropic time (ms) E ne r g y II Figure 7: Showing the evolution of energy in the physical domain, as represented by (cid:107) (cid:112) v + v (cid:107) ∞ for the isotropic material and material II.makes it more sensitive measure to any reflection. The last three materials in Table 1 violate the stability conditions as described by byB´ecache et al [25]. For these, the plane wave analysis was used in order to study thestability. This approach assumed that all the coefficients of the PDE, including α j and β j are constant throughout the PML. In spite of these assumptions, the planewave analysis provides a valuable guide for achieving stability in the discrete variable-coefficients problem [12, 13, 25]. The imaginary parts of the roots, (cid:61) { ω ( k ) } , of (26) were numerically obtained, usingMATLAB, over a range of wavevectors appropriate to our analysis. Since the materialsbeing considered are orthotropic, it is sufficient to study the first quarter of the k − space[25]. As discussed earlier, the stretch function parameters were assumed to be constants.For all cases, β that corresponds to a reflection coefficient R = 1 × − was used.25aterial III is the most challenging in terms of stability [12, 13, 25, 32] since itseverely violates the geometric stability as expressed in (27). This is evident fromthe slowness curve of Figure 2. For this material the effect of coordinate stretching,making α >
1, was examined in detail and reported in Figure 4 which was discussedin section 4, and Figure 8 which will be discussed below. k × (m ) k × ( m ) A × (m ) B I m ag i na r y pa r t o f × ( s ) Figure 8: Illustrating the effect of incorporation of the scaling coefficient α on thediscrete stability. The color maps show the continuous imaginary part ofthe unstable pair of roots, (cid:61) { ω } , for material III. (A) For α = 1. (B)For α = 10. The dashed lines indicate the highest wavevectors that canbe resolved for the mesh size used in the discrete simulations (see text fordetails).Figure 8 contain two panels each of which shows the imaginary part of the unstablepair of roots of (26). Panel (A) corresponds to using the classical stretch function, α = 1, while in (B) α = 10 was used. As one would expect, in (B) the roots wereshifted to even higher wavevectors than in the case of α = 2 in Figure 4 (B). Though,the continuous problem still unstable because the positive imaginary part only shifted.But our interest is in discrete solutions so that the question now arises as to what wouldbe the effect of this shift on the discrete problem.To answer this question, we note that the highest spatial frequency that can benumerically resolved in each direction is π / h . Dashed lines are included in both26raphs of Figure 8 to represent this threshold. It is clear from (A) that unstable rootswith positive imaginary part are present in the wavevectors range that can be resolvedby discrete models, i.e., below the dashed lines. Hence, we expect the FEM simulationsto be unstable for this case. On the other hand in (B), the unstable roots are shiftedbeyond the wavevectors range that can be discretely resolved. Therefore substantialincrease in the stability of the FEM simulations is expected. Similar results were alsoobtained for the x direction but, because the violation in the x direction for thismaterial is very severe a higher value for α was needed to ensure stability over thesame range of wavevectors.Similar plane-wave analyses were performed for materials IV, and V. For material IV,even with α = 1, the unstable pair of roots were found to occur at higher wavevectorsthan those that can be numerically resolved and hence, these should be stable in theFEM simulations. For material V, the unstable pair were below the dashed line overfor α = 1, suggesting the possibility of a numerical instability. For the discrete FEM simulation, α j and β j are not constants, rather they are functionsof x j as shown in (15) and (16). Since the unstable modes are usually the quasi-shearmodes [25], the media was excited by tangential vibrations of the cylinder surfacein order to have most of the wave energy in that mode. Figure 9 shows the FEMresult for the three unstable materials using the classical stretch function, i.e., withoutintroducing any scaling coefficients. This was achieved by setting ˜ α j = 1 in (15). In(17) the reflection coefficients were chosen to be R j = 1 . − , and in (15) and (16) m = n = 2 were used. Each row in this figure shows three snapshots for the wavepropagating in materials III, IV, and V, respectively. In the last column, to bettershow the amount of energy that remains in the computational domain, a dB scale has27een used. As expected from the plane wave analysis Figure 9 (F) shows that evenafter a long time ( 20 ms), material IV is stable. On the other hand, for material V, asshown in (I), some instabilities have emerged in PML region. Material III shows seriousinstabilities that appear to start after the arrival of the slow wave to the PML region( ∼ x ( mm ) time = 2.0 (ms) A M a t e r i a l III −
505 01 4.3 (ms) B C dB 600700 x ( mm ) D M a t e r i a l I V −
505 01 3.5 (ms) E F dB − x (mm) x ( mm ) G M a t e r i a l V − −
505 01 x (mm)3.9 (ms) H − x (mm)20.0 (ms) I − − Figure 9: Snapshot images showing the waveforms, originating from same cylinder asshown in Figure 5, but propagating in three different anisotropic solid media,namely III, IV, and V as specified in Table 1. The middle column snapshottimes were chosen to approximately correspond to the quasi-shear wave beingabsorbed by the PML. The color maps on the third column are in decibelscale.Figure 10 shows propagation snapshots for materials III and V at the same times28s in Figure 9, but with the value ˜ α = 1, ˜ α = 10 for V, and ˜ α = 20, ˜ α = 90, and m = n = 8 for III. Note that α j changes from 1 to ˜ α j though the PML, hence, higherorder polynomial were used for high ˜ α j in order to get smoother change in the PDEcoefficients at the interface between the physical domain and the PML. The comparisonof these two figures shows the effect of increasing the scaling parameter of the stretchfunction on the stability. While the instabilities disappeared for all directions in materialV and in the x direction for material III, some instability remained in the x directioncausing some energy to be reflected back to the physical domain. This is likely dueto the severity of the violation of the geometric stability in the x direction for thismaterial. Nevertheless, comparing Figure 10 (C) and Figure 9 (C) (noting the use ofdB scales), the use of a higher value for the scaling coefficient, α j results in a majorimprovement in stability for material III. This conclusion is also evident in Figure 11that shows the manner in which the energy in the physical domain evolves in time asrepresented by (cid:13)(cid:13)(cid:13)(cid:112) v + v (cid:13)(cid:13)(cid:13) ∞ in materials III and V in both cases. Using PML approach we have addressed the problem of wave propagation in an un-bounded, linear anisotropic solid in two dimensions. A time-domain second order PDEhas been derived using complex coordinate stretching. An important advantage of ourformulation is the small number of equations. Specifically, two second order equationsalong with four auxiliary equations which, to the best of knowledge, is the smallestnumber so far reported to describe wave propagation in solids using a time-domainPML formulation. This simplifies the problem and reduces the computational resourcesneeded. Moreover, by reducing the formulation to a second order, use can be made ofa wider variety of second order numerical schemes.29 x ( mm ) time=2.0 (ms) A M a t e r i a l III −
505 01 4.3 (ms) B C dB − x (mm) x ( mm ) D M a t e r i a l V − −
505 01 x (mm)3.9 (ms) E − x (mm)20.0 (ms) F − − Figure 10: Propagation snapshots as in Figure 9, but just for materials III and V, afterintroducing the scaling coefficient, ˜ α j . For material III ˜ α = 20 and ˜ α = 90.For material V ˜ α = 10 and ˜ α = 1. Comparison with Figure 9 shows thestability improvement for both materials.With help of the plane-wave analysis, we were able to stabilize the discrete PMLproblem for a wide range of otherwise unstable anisotropic media. This was achievedby increasing the value of the scaling parameter ˜ α j sufficiently to move the unstableroots out of the discretely resolved range of spatial frequencies. Only two parametersstretch function was used in our formulation, while more parameters are usually usedin formulations that were reported with methods to stabilize the problems. Whileachieving one the best reported results in stabilizing the PML problem, our method hasthe advantage of being simple. Discrete stability can be simply improved by increasingthe value of the scaling parameter. 30
10 15 2010 − A E ne r g y − − − B − − − time (ms) E ne r g y C − − − time (ms) D Figure 11: Showing the evolution of energy in the physical domain, as represented by (cid:107) (cid:112) v + v (cid:107) ∞ , for the same numerical experiments in figures Figure 9 andFigure 10. (A) Material III, for the case of ˜ α j = 1 and (B) Material III, forthe case of ˜ α = 20 , ˜ α = 90.(C) Material V, for the case of ˜ α j = 1 and (D)Material V, for the case of ˜ α = 10 , ˜ α = 1. Acknowledgements
The authors wish to thank Prof. Adrian Nachman and Prof. Mary Pugh of the Univer-sity of Toronto Department of Mathematics for their helpful advice. RSCC is gratefulto the Natural Sciences and Engineering Council (NSERC) for support under grant eferences [1] B. Engquist and A. Majda, “Absorbing Boundary Conditions for the NumericalSimulation of Waves,”
Math. Comput. , vol. 31, no. 139, pp. 629–651, 1977.[2] J.-P. B´erenger, “A perfectly matched layer for the absorption of electromagneticwaves,”
J. Comput. Phys. , vol. 114, no. 2, pp. 185–200, 1994.[3] Z. S. Sacks, D. M. Kingsland, and R. Lee, “A perfectly matched anisotropicabsorber for use as an absorbing boundary condition,”
IEEE Trans. AntennasPropag. , vol. 43, no. 12, pp. 1460–1463, 1995.[4] J. A. Roden and S. D. Gedney, “Efficient implementation of the uniaxial-basedPML media in three-dimensional nonorthogonal coordinates with the use of theFDTD technique,”
Microwave Opt. Technol. Lett. , vol. 14, no. 2, pp. 71–75, 1997.[5] S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for thetruncation of FDTD lattices,”
IEEE Trans. Antennas Propag. , vol. 44, no. 12, pp.1630–1639, 1996.[6] W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modifiedmaxwell’s equations with stretched coordinates,”
Microwave Opt. Technol. Lett. ,vol. 7, no. 13, pp. 599–604, 1994.[7] M. Kuzuoglu and R. Mittra, “Frequency dependence of the constitutive parametersof causal perfectly matched anisotropic absorbers,”
IEEE Microw. Guided WaveLett. , vol. 6, no. 12, pp. 447–449, 1996.[8] F. L. Teixeira and W. C. Chew, “On causality and dynamic stability of perfectlymatched layers for FDTD simulations,”
IEEE Trans. Microwave Theory Tech. ,vol. 47, no. 6, pp. 775–785, 1999. 329] K. C. Meza-Fajardo and A. S. Papageorgiou, “A Nonconvolutional, Split-Field,Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic ElasticMedia: Stability Analysis,”
Bull. Seismol. Soc. Am. , vol. 98, no. 4, pp. 1811–1836,2008.[10] J.-P. B´erenger, “Application of the CFS PML to the absorption of evanescentwaves in waveguides,”
IEEE J. Sel. Areas Commun. , vol. 12, no. 6, pp. 218–220,2002.[11] E. B´ecache, P. G. Petropoulos, and S. D. Gedney, “On the long-time behavior ofunsplit perfectly matched layers,”
IEEE Trans. Antennas Propag. , vol. 52, no. 5,pp. 1335–1342, 2004.[12] K. Duru and G. Kreiss, “A well-posed and discretely stable perfectly matched layerfor elastic wave equations in second order formulation,”
Commun. Comput. Phys. ,vol. 11, no. 5, pp. 1643–1672, 2012.[13] D. Appel¨o and G. Kreiss, “A new absorbing layer for elastic waves,”
J. Comput.Phys. , vol. 215, no. 2, pp. 642–660, 2006.[14] F. Collino and C. Tsogka, “Application of the perfectly matched absorbing layermodel to the linear elastodynamic problem in anisotropic heterogeneous media,”
Geophysics , vol. 66, no. 1, pp. 294–307, 2001.[15] F. Hastings, J. B. Schneider, and S. L. Broschat, “Application of the perfectlymatched layer (PML) absorbing boundary condition to elastic wave propagation,”
J. Acoust. Soc. Am. , vol. 100, no. 5, pp. 3061–3069, 1996.[16] F. H. Drossaert and A. Giannopoulos, “Complex frequency shifted convolutionPML for FDTD modelling of elastic waves,”
Wave Motion , vol. 44, no. 7-8, pp.593–604, 2007. 3317] W. C. Chew and Q.-H. Liu, “Perfectly matched layers for elastodynamics: A newabsorbing boundary condition,”
J. Comput. Acoust. , vol. 4, no. 4, pp. 341–359,1996.[18] S. Kucukcoban and L. F. Kallivokas, “Mixed perfectly-matched-layers for directtransient analysis in 2D elastic heterogeneous media,”
Comput. Meth. Appl. Mech.Eng. , vol. 200, no. 1-4, pp. 57–76, 2011.[19] D. Appel¨o, T. Hagstrom, and G. Kreiss, “Perfectly Matched Layers for HyperbolicSystems: General Formulation, Well-posedness, and Stability,”
SIAM J. Appl.Math. , vol. 67, no. 1, pp. 1–23, 2006.[20] J. A. Roden and S. D. Gedney, “Convolution PML (CPML): An efficient FDTDimplementation of the CFS-PML for arbitrary media,”
Microwave Opt. Technol.Lett. , vol. 27, no. 5, pp. 334–339, 2000.[21] P. R. Loh, A. F. Oskooi, M. Ibanescu, M. Skorobogatiy, and S. G. Johnson, “Fun-damental relation between phase and group velocity, and application to the failureof perfectly matched layers in backward-wave structures,”
Phys. Rev. E , vol. 79,no. 6, 2009.[22] G. Kreiss and K. Duru, “Discrete stability of perfectly matched layers foranisotropic wave equations in first and second order formulation,”
BIT Numer.Math. , vol. 53, no. 3, pp. 641–663, Mar. 2013.[23] D. Komatitsch and R. Martin, “An unsplit convolutional perfectly matched layerimproved at grazing incidence for the seismic wave equation,”
Geophysics , vol. 72,no. 5, p. SM155, 2007.[24] A. I. Beltzer,
Acoustics of solids . New York: Springer-Verlag, 1988.3425] E. B´ecache, S. Fauqueux, and P. Joly, “Stability of perfectly matched layers, groupvelocities and anisotropic waves,”
J. Comput. Phys. , vol. 188, no. 2, pp. 399–433,2003.[26] S. G. Johnson. (2008) Notes on Perfectly Matched Layers ( PMLs), MIT OpenCourse Ware, .[27] F. L. Teixeira and W. C. Chew, “Complex space approach to perfectly matchedlayers: a review and some new developments,”
Int. J. Numer. Modell. Electron.Networks Devices Fields , vol. 13, no. 5, pp. 441–455, 2000.[28] W. Zhang and Y. Shen, “Unsplit complex frequency-shifted PML implementa-tion using auxiliary differential equations for seismic wave modeling,”
Geophysics ,vol. 75, no. 4, pp. T141–T154, 2010.[29] F. H. Drossaert and A. Giannopoulos, “A nonsplit complex frequency-shifted PMLbased on recursive integration for FDTD modeling of elastic waves,”
Geophysics ,vol. 72, no. 2, p. T9, 2007.[30] P. G. Petropoulos, “Reflectionless Sponge Layers as Absorbing Boundary Condi-tions for the Numerical Solution of Maxwell Equations in Rectangular, Cylindrical,and Spherical Coordinates,”
SIAM J. Appl. Math. , vol. 60, no. 3, pp. 1037–1058,2000.[31] C. H. Daros, “Material Stability Conditions for a Class of InhomogeneousAnisotropic Media,”
Math. Mech. Solids , vol. 14, no. 4, pp. 377–389, 2007.[32] Y. Li and O. B. Matar, “Convolutional perfectly matched layer for elastic second-order wave equation,”