Boundary waves and stability of the perfectly ... - Semantic Scholar
Boundary waves and stability of the perfectly matched layer Kenneth Duru ∗, Gunilla Kreiss
†
April 10, 2012
Abstract We study the stability of the perfectly matched layer (PML) for symmetric second order hyperbolic partial differential equations on the upper half plane, with boundary conditions at y = 0. Using a mode analysis, we develop a technique to analyze the stability of the corresponding initial boundary value problem on the half plane. We apply our technique to the PML for the elastic wave equations subject to free surface and homogenous Dirichlet boundary conditions, and to the PML for the curl–curl Maxwell’s equation subject to insulated walls and perfectly conducting walls boundary conditions. The conclusion is that these half–plane problems do not support temporally modes. Keywords: elastic waves, electromagnetic waves, boundary waves, surface waves, glancing waves, perfectly matched layer, stability, normal mode analysis.
1
Introduction
There are many wave propagation problems where boundary phenomenon such as surface or glancing waves are dominant. Typical examples are in non-destructive testing, seismology, earthquake engineering, ultrasonics, and ground penetrating radar (GPR) technologies. These problems can be described by symmetric time–dependent partial differential equations (PDE) ∗
Division of Scientific Computing, Department of Information Technology, Uppsala University Sweden † Division of Scientific Computing, Department of Information Technology, Uppsala University Sweden
1
in semi-bounded domains. In order to perform numerical simulations, the unbounded extensions of the domains must be truncated, by introducing artificial boundaries. Thus efficient and reliable domain truncation is vital, since it enables more accurate numerical simulations. Perhaps, an efficient approach to truncate unbounded domains is to surround the artificial boundaries with the perfectly matched layer (PML) [2, 17, 26, 28]. The PML consists in extending the computational domain to a layer, where the underlying equations are modified such that waves decay rapidly in the layer. A very important property of the PML is perfect matching. This means that waves propagate from the physical space into the PML without reflections. A PML is usually derived by assuming a homogeneous media and an infinite domain in all directions. In a domain with physical boundaries, the application of the PML introduces boundary corners where physical boundary conditions interact with the PML. In order to ensure accuracy, the underlying boundary conditions must be accurately extended from the interior into the PML. Thus leading to an initial boundary value problem (IBVP) for the PML. These types of problems can be found many application areas, including seismic, optical and gravitational waves. Examples include glancing waves and surface waves in electromagnetic and elastic wave propagation problems. By construction the PML is perfectly matched, but there is no guarantee that all solutions decay with time. The analysis of temporal stability of PMLs is therefore a main topic of research. For Cauchy problems, the stability of the PML can be predicted by the so–called geometric stability condition [7] for several hyperbolic systems [25, 24, 7, 12]. In several studies, numerical experiments have shown that PMLs which are Cauchy stable can exhibit (numerical) growth when some physical boundaries are imposed [9, 27, 1]. For first order hyperbolic PDEs, the theory of IBVPs is rather well–developed. The Friedrichs theory [35, 36, 37] for symmetric systems with maximally dissipative boundary conditions was available, but not sufficient for nonsymmetric systems or more general boundary conditions. A more general theory is the eigenvalue theory based on the normal mode analysis [38]. H–O. Kreiss and his co–workers developed this theory in a more general and systematic manner [33, 34, 39]. Therefore this theory is often called the Kreiss theory/technique. By the Kreiss theory, all boundary stable problems must satisfy the so–called determinant condition [33, 34, 39]. This theory has recently been extended to second order hyperbolic systems [5]. However, the PML is a more complicated system. The PML is non–standard; it is an integro–differential equation. Of course by introducing auxiliary variables the PML can be rewritten as a system of PDEs. For second order 2
systems, the PML is a combination of second order PDEs for the modified wave equation and a first order system of PDEs for the auxiliary differential equations. Unfortunately no stability theory exists for the corresponding IBVP including the PML. In this paper we consider second order symmetric hyperbolic systems on the upper half–plane with boundary conditions at y = 0. First we derive a PML that truncates the boundary in the x–direction, by a complex change of variables ∂/∂x → 1/Sx ∂/∂x. Here, Sx is the PML (complex) metric. We also transform the x–derivatives (∂/∂x) in the boundary condition with the PML metric Sx (∂/∂x → 1/Sx ∂/∂x). Choosing auxiliary variables, we construct new transformed boundary conditions in the PML. Using a mode analysis we develop a stability theory for the corresponding IBVP. A main result is that if the characteristic function for the undamped problem is a homogeneous function1 , then the introduction of the PML will immediately move all roots of the characteristic equation into the negative half complex– plane. We apply our technique to the PML for the elastic wave equation subject to free surface and homogeneous Dirichlet boundary conditions, and to the PML for the curl–curl Maxwell’s equation subject to insulated walls and perfectly conducting walls boundary conditions. The conclusion is that these half–plane problems do not support temporally growing modes. The outline of the paper is as follows. In Section 2, we introduce our model problem and discuss the stability analysis using the normal mode analysis. In Section 3, we derive the PML and the modified boundary conditions. In Section 4, we develop a technique for the stability analysis of the corresponding IBVP for the PML and formulate our main results. In sections 5 and 6, we apply our technique to the elastic wave equations and the Maxwell’s equations, respectively. In the last section we draw conclusions and suggest future work.
2
Symmetric hyperbolic systems
Consider the second order hyperbolic system of PDEs in two space dimensions ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ T ∂u A + B + C + C . (1) = ∂t2 ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x Here, the unknown u = (u1 , u2 , · · · , um )T , is a real valued vector function, (x, y) are the spatial coordinates and t denotes time. The coefficients A, B, C are real matrices and define the media. 1
Let f (v) be a function with the vector argument v. If f (αv) = αk f (v) with α 6= 0 and k ∈ Z, then f (v) is a homogeneous function of degree k.
3
Definition 1 The system (1) is called symmetric hyperbolic if the symbol Ψ (kx , ky ) = kx2 A + ky2 B + kx ky C + C T ,
(2)
is symmetric and positive semi-definite for all kx , ky ∈ R. In general the coefficient matrices A, B are symmetric and positive semi– definite. If (1) is symmetric hyperbolic the eigenvalues of the symbol Ψ(kx , ky ) are real and non–negative, and the eigenvectors are orthogonal. We set the smooth initial data u (x, y, 0) = f (x, y) ,
2.1
∂u (x, y, 0) = g (x, y) . ∂t
Plane waves and slowness diagrams
In order to understand the wave propagation properties of the model (1) it is useful to consider wave-like solutions u = u0 est−ikx x−iky y ,
u0 ∈ Rn ,
kx , ky , x, y ∈ R,
0 ≤ t.
(3)
In equation (3), (kx , ky ) ∈ R2 is the wave vector, and u0 is a vector of constant amplitude called the polarization vector. By inserting (3) in (1) we have the following solvability condition, often called the dispersion relation F0 (s, kx , ky ) ≡ det s2 I + Ψ (kx , ky ) = 0, Ψ (kx , ky ) = kx2 A + ky2 B + kx ky C + C T .
(4)
Note that the symbol Ψ(kx , ky ) is a symmetric matrix. The polarization vector u0 is the eigenvector of the matrix Ψ(kx , ky ), and the eigenvalue is −s2 . It should be noted that if (1) is symmetric hyperbolic, then the roots s of the dispersion relation (4) are distinct and purely imaginary. These roots are called the physical roots and we will refer to the corresponding modes as the physical modes. Since the physical roots are distinct and purely imaginary, and the eigenvectors u0 are orthogonal, the Cauchy problem (1) is well-posed2 . If we write s = iω(kx , ky ), where ω ∈ R is called the temporal frequency, we 2
We also note that if (1) is symmetric hyperbolic with periodic boundary conditions, it is possible to derive a strict energy estimate.
4
can also introduce kx ky , ), the normalized propagation direction, |k| |k| ω ω Vp = ( , ), the phase velocity, kx ky kx ky S = ( , ), the slowness vector, ω ω Vg = ∇k ω(kx , ky ), the group velocity, q and |k| = kx2 + ky2 . K=(
(5)
By inspection the dispersion relation (4) is homogeneous in ω, kx , ky , we can rewrite (4) as F0 (i, S) = 0.
(6)
The wave propagation properties of a certain media can be described by plotting the slowness diagrams defined by points in the S−plane satisfying (6). The group velocity can be expressed as ∂F (iω, k , k ) −1 0 x y Vg = ∇k ω(kx , ky ) = − ∇k F0 (iω, kx , ky ). ∂ω
(7)
Thus, the group velocity is normal to the slowness curve. We are now ready to introduce the geometric stability condition [7]. Definition 2 (Geometric stability condition) Let Vg = ((Vg )x , (Vg )y ) and S = (Sx , Sy ) denote the Cartesian components of the group velocity and the slowness curve respectively. The geometric stability condition in the x−direction is Sx × (Vg )x ≥ 0, for all points on the slowness curve. By Theorem 2 in [7], the geometric stability in the x–direction is necessary for stability of the vertical PML.
2.2
Boundary conditions and the normal mode analysis
Now consider (1) on the upper half–plane Ω = {(x, y) | − ∞ < x < ∞, 0 ≤ y < ∞}. At y = 0 we set the boundary condition ∂ ∂ ∂ L , , u = 0. (8) ∂t ∂x ∂y Here L is a first order partial differential operator acting only at the boundary.
5
For some boundary operators L, it is possible to establish well–posedness by deriving an energy estimate, see [13]. We will use the normal mode analysis [5, 6] to investigate well–posedness and temporal growth properties of the PML. To begin with, we introduce the L2 –norm Z ∞ 2 kv (y) k = v∗ vdy,
(9)
0
where v∗ denotes the conjugate transpose of v. Lemma 1 The problem (1) with the boundary condition (8), is not well– posed if there are non-trivial solutions of the form ˆ (y) , u = est+ikx x u
kˆ u (y) k < ∞,
0 ≤ y < ∞,
(10)
with <s > 0. Proof : If (10) is a solution then ˆ (βy) , uβ = esβt+ikx βx u
kˆ u (βy) k < ∞,
(11)
is also a solution, for any β > 0. Thus, if <s > 0, then we can always construct solutions (11) which can grow arbitrarily fast in time. Thus the problem is ill–posed. We will demonstrate how to use the normal mode analysis to reformulate Lemma 1 as an algebraic condition, the so–called determinant condition. Consider the symmetric hyperbolic system (1) in the upper half plane and the boundary condition (8) at y = 0. We take the Laplace transform in time and the Fourier transform in the tangential direction x. We obtain an eigenvalue problem ˆ = − kx2 Aˆ s2 u u + ikx C + C T
dˆ ˆ u d2 u +B 2, dy dy
d ˆ ˆ = 0, L s, ikx , u dy
y = 0,
0 ≤ y < ∞, kˆ u (y) k < ∞.
<s > 0, (12) (13)
Here (s, kx ) are the dual variables to (t, x). For the ordinary differential equation, ODE, (12) we seek solutions of the form ˆ = Φeκy , u and we have s2 I + kx2 A − ikx κ C + C T − κ2 B Φ = 0. 6
(14)
Note that Φj (s, kx , κ) are the eigenvectors of the matrix P (s, kx , κ) = s2 I + kx2 A − ikx κ C + C T − κ2 B,
(15)
and the corresponding eigenvalues θj (s, kx , κ) must satisfy θj (s, kx , κ) = 0.
(16)
Observe that if we set κ = iky , we get exactly the dispersion relation (4). The matrix P (s, kx , κ) is complex symmetric and the eigenvectors Φj (s, kx , κ) are complete. The eigenvalues θj (s, kx , κ) are quadratic functions of κj , and the complex roots κj of (16) come in pairs (+)
κj (s, kx ),
(+)
0,
(−)
0,
< 0,
where the superscript (+) denotes the root with positive real part and the superscript (−) denotes the root with negative real part. The general solution of (12) is ˆ (y) = u
m X
(+)
(+)
δj Φj eκj
y
+
j
m X
(−)
(−)
δj Φj eκj
y
.
(18)
j (+)
We consider only bounded solutions in L2 (kˆ u (y) k < ∞), then δj and we have ˆ (y) = u
m X
δj Φj eκj y ,
0.
(21)
We have Theorem 1 Consider the symmetric hyperbolic system (1) in the upper half plane and the boundary condition (8) at y = 0. The problem (1) with (8) is not well–posed if there are s with <s > 0 satisfying the determinant condition F0 (s, kx ) ≡ det Υ (s, kx ) = 0. 7
Theorem 1 is a very important result. The whole process of finding IBVPs which are not well–posed reduces to solving the algebraic equation (21). Note that most relevant physical problems are well–posed. Examples are clamped wall and free–surface boundary conditions in elastodynamics, hard wall and soft wall boundary conditions in acoustics, and insulated wall and perfectly conducting wall boundary conditions in electromagnetism. For the above examples there are energy estimates see [13], thus their determinant condition (21) is satisfied. It is important to note that using the normal mode analysis there can be a technical difficulty of establishing sharp estimates in the physical space, especially for roots corresponding to generalized eigenvalues, see [39, 5]. From now on we will assume that the system (1) is symmetric hyperbolic. We are only interested in physical problems which are well–posed. Therefore we will also assume that the problem (1) with the boundary condition (8) is well–posed. Note that if (1) with (8) is well–posed, then for all roots sj of F0 (s, kx ) = 0 we can write sj (kx ) = −aj0 (kx ) + ibj0 (kx ) , where
aj0 (kx ) ≥ 0,
aj0 (kx ) , bj0 (kx ) ∈ R,
j = 1, . . . , 2m.
(22)
Furthermore, if F0 (s, kx ) is homogeneous of degree 2m we can introduce the scaling s kx η0 = , k10 = ≡ ±1. (23) |kx | |kx | Note that all roots ηj0 of F0 (η 0 , k10 ) = 0 are complex numbers, and satisfy ηj0 = −a0j0 + ib0j0 , where
a0j0 ≥ 0,
a0j0 , b0j0 ∈ R, j = 1, . . . , 2m.
(24)
We can distinguish the following cases; • When 0, the problem (1) with (8) is called boundary stable. If there are non–homogeneous boundary data, the solutions on the boundary can be estimated in terms of the data. • When 0, 0 ≤ y < ∞. The PML will be derived using the complex coordinate stretching technique [11]. In the complex metric we will also introduce the complex frequency shift [10]. To begin, we take the Laplace transform in time Z ∞ ˆ (x, y, s) = u e−st u(x, y, t)dt. (25) 0
We then analytically continue the transformed wave equation on the complex contour, Z x σ(z) Υ(x) = x + dz, α +s 0 where oscillating solutions are turned into exponentially decaying solutions. Here, σ(x) ≥ 0 is the damping function and α ≥ 0 is the complex frequency shift (CFS). We denote the complex metric Sx :=
dΥ(x) σ(x) =1+ . dx α+s
If α = 0 we call it the standard PML and if α > 0 we call it the CFS-PML. Changing the coordinate x → Υ(x), we have ˆ ˆ ˆ ˆ 1 ∂ 1 ∂u ∂u ∂ ∂u 1 T ∂u ˆ= s2 u A +C + B + C . (26) Sx ∂x Sx ∂x ∂y ∂y ∂y Sx ∂x Multiplying (26) by Sx leads to ˆ ˆ ˆ ˆ ∂ 1 ∂u ∂u ∂ ∂u 2 T ∂u ˆ= s Sx u A +C + Sx B +C . ∂x Sx ∂x ∂y ∂y ∂y ∂x
(27)
Observe that in (26) the flux in the y-direction changes from B and from B
ˆ ˆ ∂u ∂u + CT ∂y ∂x
ˆ ˆ ∂u 1 T ∂u + C ∂y Sx ∂x
to B
ˆ ˆ ∂u 1 T ∂u + C , ∂y Sx ∂x
to Sx B
ˆ ˆ ∂u ∂u + CT , ∂y ∂x
in (27). Choosing the auxiliary variables ˆ= v
ˆ 1 ∂u , (s + α) Sx ∂x
ˆ = w 9
ˆ 1 ∂u , α + s ∂y
ˆ= q
αˆ u , α+s
(28)
and inverting the Laplace transforms yielding ∂2u ∂u ∂ ∂u ∂u ∂ ∂u T ∂u A + B +σ = +C +C ∂t2 ∂t ∂x ∂x ∂y ∂y ∂y ∂x ∂ (σv) ∂(σw) −A +B + σα (u − q) , ∂x ∂y ∂v ∂u = − (σ + α) v + , ∂t ∂x ∂w ∂u = − αw + , ∂t ∂y ∂q =α (u − q) . ∂t
(29)
Away from the boundary at y = 0, the PML model (29) is perfectly matched to the system (1) if σ (0) = 0. There are no reflections as waves pass the interface x = 0. That is, waves of all frequencies and angles are transmitted from the interior into the layer without reflections. The PML damps waves so that they become insignificant before they reach the outer boundary of the computational domain. We initialize the fields in the PML with homogeneous initial data u(x, y, 0) = 0, v(x, y, 0) = 0,
3.1
∂u (x, y, 0) = 0, ∂t
w(x, y, 0) = 0,
q(x, y, 0) = 0.
Boundary conditions
The accuracy of the PML is ensured by the perfect matching property. This can be shown mathematically by constructing the general solutions to the wave equation and the PML equation in the Laplace-Fourier space, see [14]. For a PML that depends on the x–direction, the analysis in [14] assumes that in the y-direction the domain is unbounded (or periodic) such that a Fourier transform in the y-direction is possible. If the underlying problem is bounded in one or more directions with non-periodic boundary conditions, a different technique is needed. To ensure perfect matching, the underlying boundary conditions must be simultaneously transformed using the PML metric Sx . Now consider the boundary condition (8) and introduce the Laplace transform (25). Still in the Laplace space we introduce the complex metric Sx , having 1 ∂ ∂ ˆ ˆ = 0. L s, u (30) , Sx ∂x ∂y In order to localize the transformed boundary condition (30) in time we introduce the auxiliary variables (28) and invert the Laplace transform, to 10
get L
∂ ∂ ∂ , , , σ, α (u, v, w, q) = 0. ∂t ∂x ∂y
(31)
Note that (31) is a general representation of the transformed boundary condition. Depending on the boundary condition not all the auxiliary variables v, w, q, will be required. Also note that the boundary condition (31) is a lower order modification of the boundary condition (8). If the damping vanishes, σ(x) ≡ 0, the auxiliary variables decouple completely and we recover the physical boundary condition (8).
4
Stability analysis
In order for the PML to be computationally useful, the PML must not support growth. Any growth in the layer can propagate into the computational domain and pollute the solution. The question of stability is not trivial, it is therefore a topic of active research. The perfectly matched layer is indeed a variable coefficient problem. In addition, boundary conditions can introduce many more challenging difficulties. However, we can localize the problem and split the problem into a Cauchy problem and a half-plane problem with constant coefficients. For the Cauchy problem, growth can be predicted by the so called geometric stability condition [7]. The boundary conditions can then be analyzed by using a mode analysis. If the geometric stability condition is not violated and the roots of the characteristics equations corresponding to the boundary conditions lie on the left half–complex plane, we conclude that the PML does not support growing solutions. If the geometric stability condition is violated and the roots of the characteristics equations corresponding to the boundary conditions lie on the left half–complex plane, we say that the PML supports bounded growth. At least the boundary conditions can not cause additional growth, and therefore are not harmful to the PML. In any case if the roots of the characteristics equations corresponding to the boundary conditions lie on the right half–complex plane, the boundary condition supports exponentially growing modes, and can be harmful to the PML. In this case, we say that the boundary condition is unstable.
4.1
Cauchy problem
To keep this paper self–contained we summarize the arguments leading to the geometric stability condition. More elaborate discussions can be found in [7, 12, 29]. To begin with, we freeze all coefficients and consider the Cauchy problem (29). We then assume solutions of the form V = V0 est−ikx x−iky y , 11
V = (u, v, w, q)T .
(32)
Inserting (32) in (29), we have the dispersion relation F (η, kx , ky , σ, α) ≡ (s + α)m F0 (s (s + σ + α) , (s + α) kx , ky (s + σ + α)) = 0. (33) For compactness, we introduce the normalized variables, η=
s , |k|
k1 =
kx , |k|
k2 =
ky , |k|
=
σ1 , |k|
γ=
α , |k|
|k| =
q kx2 + ky2 , (34)
and we have F (η, k1 , k2 , , γ) ≡ (η + γ)m F0 (η (η + + γ) , (η + γ) k1 , k2 (η + + γ)) = 0, (35) where F0 (η0 , k1 , k1 ) is given by (4). The system (29) is stable if for any (k1 , k2 ) satisfying k12 + k22 = 1, the roots η(k1 , k2 , , γ) have negative real parts. The existence of solutions of roots η with positive real parts would correspond to plane wave solutions with exponential growth in time. A stable system does not admit such solutions. We will use standard perturbation arguments. By setting = 0 in (35) we have F (η, k1 , k2 , 0, γ) = (η + γ)3m F0 (η, k1 , k2 ) . (36) Clearly (36) has 5m roots, 2m of which are the purely imaginary roots of (4) and correspond to the physical modes. The remaining 3m roots are real and non-positive η = −γ ≤ 0. (37) We call the corresponding modes non-physical. Since the roots depend continuously on the coefficients, it is apparent that at intermediate frequencies and with small enough damping, the roots η of the non-physical modes will remain in the left half of the complex plane. The following result was proven in [12]. Lemma 2 (Lemma 2 in [12]) At sufficiently high frequencies, the parameter γ > 0 will stabilize the non-physical modes for any ≥ 0. We see that if γ > 0, then the instability in the PML at constant coefficients can not come from the non-physical modes. The physical modes are simple (distinct), and can therefore be expanded in powers of , η = η0 + η + O 2 . (38) Here, η0 is a purely imaginary root of F0 (η, k1 , k2 ). At sufficiently high frequencies, 0. If the geometric stability condition in the x-direction is violated, then there are unstable physical modes at sufficiently high frequencies. We see that away from the boundaries, if the geometric stability condition is not violated the PML does not support growth. In order to understand the effect of the boundary conditions we perform the modal analysis as in Section 2.2.
4.2
Half plane problem
Now consider the constant coefficients time-dependent PML (29) in the upper half plane with the modified boundary condition (31) at y = 0. We take the Laplace transform in time and the Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables and we have dˆ ˆ u d2 u ˆ = − k˜x2 Aˆ s2 u u + ik˜x C + C T + B 2 , 0 < y < ∞, dy dy ∂ ˆ = 0, y = 0, u Lˆ s, ik˜x , ∂y where
<s > 0, (39)
(40)
kx k˜x = . Sx
Here (s, kx ) are the dual variables to (t, x). Introducing the modified wave– number k˜x simplifies the remaining analysis. Note that equations (39) and (40) are completely analogous to equations (12) and (13). Thus the remaining part of the analysis can be adapted from the analysis in Section 2.2. Repeating all the steps leading to theorem 1, we arrive at the determinant condition F0 s, k˜x ≡ det Υ s, k˜x 6= 0, <s > 0. (41) The stability of the boundary condition (31) at y = 0 can be formulated as an algebraic condition, Definition 3 Consider the PML (29) in the upper half plane subject to the modified boundary condition (31) at y = 0. The problem (29) with (31) is not well–posed if for any σ ≥ 0 there is a root s with <s > 0 satisfying the determinant condition s+α ˜ kx = 0. (42) F0 s, kx ≡ F0 s, s+α+σ 13
We can prove the following Theorem 2 Consider the standard PML (29) (with α ≡ 0), and the modified boundary condition (31) at y = 0. Assume that the undamped problem is well–posed, and consequently all roots s0j (kx ) of F0 (s0 , kx ) = 0 satisfy <s0j (kx ) ≤ 0 If the characteristic polynomial F0 (s0 , kx ) is homogeneous of degree 2m in both s0 and kx , then the roots of (42) are s = 0, sj (kx ) = −σ + s0j (kx ) ,
j = 1, . . . , 2m.
Proof : Introduce α = 0 in (42), ince F0 s, k˜x is homogeneous we write F0 s,
s kx s+σ
≡
s s+σ
2m F (s + σ, kx ) = 0.
Therefore, we have a multiple root s = 0 and F0 (s + σ, kx ) = 0. The corresponding roots of F0 (s + σ, kx ) = 0 are sj (kx ) = −σ + s0j (kx ) ,
j = 1, . . . , 2m.
Consider now the CFS–PML with α > 0. Note that if F0 (s, kx ) is homogeneous of degree 2m in both s and kx , we have 2m s+α s+α s+α+σ F0 s, kx ≡ F0 s , kx = 0. (43) s+α+σ s+α+σ s+α In order to ensure compactness, we introduce the normalized variables, η0 =
s , |kx |
we obtain
k10 =
η0 + γ 0 0 η + γ 0 + 0
kx ≡ ±1, |kx |
2m
F0 η
0η
0 = 0
σ1 , |kx |
+ γ 0 + 0 0 , k1 η0 + γ 0
γ0 =
α , |kx |
(44)
Note that since η 0 = −γ 0 does not satisfy (42) we have 2m η0 + γ 0 6= 0. η 0 + γ 0 + 0
= 0.
(45)
(46)
Therefore we consider 0 0 0 0η + γ + 0 F0 η , k1 = 0. η0 + γ 0
We have the result 14
(47)
Theorem 3 Consider the CFS–PML (29) (with α > 0) subject to the modified boundary condition (31) at y = 0, and the characteristics polynomial (47). Assume that the undamped problem is well–posed, that is all roots s0j of F0 (s0 , kx ) satisfy s0j = −aj0 + ibj0 ,
aj0 , bj0 ∈ R,
where
aj0 ≥ 0,
j = 1, . . . , 2m.
If the characteristic polynomial F0 (s0 , kx ) is homogeneous of degree 2m in 0 both s0 and kx , then for each root s0j of F0 (s0 , kx ), there are two roots, η1j 0 0 0 and η2j , of (47). The two roots, η1j and η2j satisfy 1) When a00j = b00j = 0, 0 = 0, η 0 = − (γ 0 + 0 ) . η1j 2j
2) When b00j = 0, 0 η1j
0
0
=− γ + +
a00j
/2 −
0 η2j = − γ 0 + 0 + a00j /2 +
r r
γ 0 + 0 + a00j
2
/4 − a00j γ 0 ,
γ 0 + 0 + a00j
2
/4 − a00j γ 0 .
3) When a00j = 0, b00j 6= 0, 0 < 0, 0 the corresponding roots have non–positive real parts. Since the roots are continuous functions of 0 > 0 there will be roots with positive real parts only if there is a purely imaginary root of (50) for some 0 > 0. We will show that this is impossible. Assume that ηj0 = iξj0 with ξj0 ∈ R, is a root for some b00j 6= 0, γ 0 , 0 > 0. Inserting ηj0 = iξj0 in (50) and considering the real and the imaginary parts separately we have −ξj02 + a00j γ 0 + b00j ξj0 = 0, γ 0 + 0 + a00j ξj0 − b00j γ 0 = 0. The two relations above are simultaneously satisfied only if γ0 = 0
or a00j = b00j = 0
or γ 0 = −(a00j
or γ 0 = 0, 0 = −a00j s a00j + 0 0 0 + ) ± ib0j . a00j
or γ 0 = −0 − a00j , b00j = 0,
Since b00j 6= 0, and γ 0 , 0 > 0, are purely real and positive these are contradictions. The proof of the theorem is complete.
This is a very useful result. If we can show that the characteristic polynomial (21) of the undamped problem is a homogeneous function we can immediately conclude that the corresponding boundary condition is stable. The introduction of the PML will move all non–zero roots, s, into the negative half complex–plane. If there are zero roots they will not be perturbed by the PML. Thus the boundary condition (31) can not be harmful to the 17
PML (29). We conclude with the following remark. Consider the half–plane PML problem (29) with (31). If the undamped problem and sj denotes is well–posed the roots of the characteristics equation F0 s, k˜x = 0, we can distinguish the following cases • no growth: The geometric stability condition is satisfied and <sj ≤ 0. • bounded growth: The geometric stability condition is violated and <sj ≤ 0. • ill–posedness: <sj > 0.
5
Linear elastodynamics
In this section, we will apply the results of Section 4 to the PML for the elastic wave equation, subject to a free–surface and a homogeneous Dirichlet boundary condition. Consider the elastic wave equation for a source free heterogeneous media in two space dimensions ρ
∂2u =∇·T, ∂t2
x ∈ Ω,
t > 0.
(57)
Here, x = (x, y) are the spatial variables, t denotes time, Ω is the upper half plane and ρ is the density of the media. The unknown u = (u1 , u2 )T , is the displacement vector and T is the stress tensor. The x and y components of the stress tensor T are ∂u ∂u ∂u ∂u Tx = A +C , Ty = B + CT . (58) ∂x ∂y ∂y ∂x The matrices A, B, C are given by the elastic coefficients. If time is scaled appropriately, then the system (57) can be rewritten as ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ T ∂u A + B + C + C . (59) = ∂t2 ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x Note that (59) is completely analogous to the model (1). In an orthotropic media whose principal axes coincide with the (x, y) axes, the fourth-order symmetric and positive definite tensor C of elastic coefficients [32] is
c11 c12 0 C = c12 c22 0 , 0 0 c33
where
c11 , c12 , c22 , c33 > 0,
18
c11 c22 − c212 > 0.
Thus, the coefficient matrices in (59) are c11 0 c33 0 A= , B= , 0 c33 0 c22
C=
0 c12 . c33 0
(60)
For the particular case of isotropic media the elastic coefficients depend on the Lam´e parameters λ, µ, c11 = c22 = 2µ + λ,
c33 = µ,
c12 = λ.
At the boundary y = 0, the boundary condition is • a free–surface: Ty := B
∂u ∂u + CT = 0, ∂y ∂x
y = 0,
(61)
or • a clamped wall: u = 0,
y = 0.
(62)
The perfectly matched layer for the elastic wave equation (59) is the same as (29). Away from the physical boundary y = 0, perfect matching is achieved if the damping function vanishes at the interface σ (0) = 0. In order to achieve perfect matching on the boundary at y = 0, we must apply the complex change of variables ∂ 1 ∂ → ∂x Sx ∂x
(63)
to the boundary conditions. Consider the free-surface boundary condition (61) in the y-direction, Ty = 0. By taking Laplace transform in time and applying the complex change of variables (63) in (61) we have B
ˆ ˆ ∂u 1 T ∂u + C = 0. ∂y Sx ∂x
(64)
In order to localize the modified boundary condition (64) in time we introduce the auxiliary variables (28). Depending on our particular choice of auxiliary variables (28) we have two possible modified boundary conditions in the PML, ∂u ∂u B + CT − σC T v = 0, (65) ∂y ∂x B
∂u ∂u + CT + σBw = 0. ∂y ∂x 19
(66)
Note that in the continuous setting the modified boundary conditions, (65) and (66), are equivalent and correspond to the modified boundary condition (64) in the Laplace space. Observe also that if the damping vanishes, σ = 0, we recover the traction–free boundary condition (61). On a clamped boundary, we set the displacement field to zero, u = 0.
(67)
Theorem 4 Consider the constant coefficient PML (29) for isotropic elastic media in the upper half-plane Ω = (−∞, ∞) × (0, ∞), and subject to the modified boundary condition (65) or (66) at y = 0, or the homogeneous Dirichlet condition (67) at y = 0. The corresponding initial boundary value problem does not support growing modes. Proof: We know that in isotropic media the geometric stability condition is always satisfied, thus by lemma 3 the corresponding constant coefficients Cauchy problem is stable. Now consider the corresponding constant coefficient half plane problem. As before, we take Laplace transform in time and Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables, leading to the system of second order ordinary differential equations dˆ ˆ d2 u u ˆ = − k˜x2 Aˆ s2 u u + B 2 + ik˜x C + C T , dy dy
(68)
where
kx k˜x = . Sx The corresponding boundary condition is • a free–surface: B
dˆ u(ik˜x , 0, s) ˆ (ik˜x , 0, s) = 0, + ik˜x C T u dy
(69)
or • a clamped wall: ˆ (ik˜x , 0, s) = 0. u
(70)
The ordinary differential equations (68) has the general solution ˆ (ik˜x , y, s) = δ1 Φ1 (ik˜x , s)e−κ1 y + δ2 Φ2 (ik˜x , s)e−κ2 y u + δ3 Φ1 (ik˜x , s)eκ1 y + δ4 Φ2 (ik˜x , s)eκ2 y . 20
(71)
For the particular case of an isotropic media κ1 , κ2 are given by s s 2 2 ˜ kx s + µkx s2 + (λ + 2µ)k˜x2 , , κ2 = ± , k˜x = κ1 = ± µ 2µ + λ Sx √ with the square root z → z being defined as the one with non–negative real part. In the appendix we have proven that if <s > 0, then 0, α ≥ 0. Next we consider the homogeneous Dirichlet boundary condition (70). We insert the solution (72) in (70) to determine, δ1 , δ2 . ˆ (ik˜x , 0, s) = δ1 Φ1 (ik˜x , s) + δ2 Φ2 (ik˜x , s), u
(78)
The corresponding linear system is 1 ˜
− iκk1x
˜x ik κ2
! δ1 = 0. δ2 1
In order to ensure a non-trivial solution we must have ! ˜x ik 1 κ2 det = 0. ˜ − iκk1x 1
(79)
(80)
Computing the determinant (80) we have s2 s2 (s + α + σ)2 + (3µ + λ) (s + α)2 kx2 F2 (s, k˜x ) ≡ k˜x2 − κ2 κ1 ≡ = 0. µ (2µ + λ) (s + α + σ)2 (81) 22
Again the characteristic equation F2 is homogeneous of degree 4 in both s, kx . At σ ≡ 0 the roots of the numerator of F2 are p s = 0, s = −α, s = ±i (3µ + λ)kx . By theorems 2 and 3 we have <s ≤ 0 when σ > 0, α ≥ 0. This completes the proof. Note that for anisotropic elastic media the expressions we obtain are in general more complicated. However, the characteristic polynomial F(s, kx ) is homogeneous in both s and kx . Thus by theorems 2 and 3, when the geometric stability condition is satisfied, the corresponding PML subject to the free–surface boundary condition (66) or the homogenous Dirichlet boundary condition (67) does not support growth.
6
Maxwell’s equations
Consider the curl–curl Maxwell’s equations in a linear, source free, homogeneous isotropic media 1 ∂2u ∂u + µσ ∗ = −∇ × ∇ × u, 2 2 c ∂t ∂t
1 c= √ . µ
(82)
Here the vector u = (u1 , u2 )T is the electric field, c is the wave speed, µ, and σ ∗ are the relative magnetic permeability, the relative electric permittivity and the conductivity of the media respectively. If time is scaled appropriately we can write r ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ µ ∂u ∗ T ∂u +σ = A + B + C + C . ∂t2 ∂t ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x (83) where the coefficient matrices are given by 1 0 0 0 A= , B= , 0 0 0 1
C=
0 0 . −1 0
(84)
At the boundary y = 0, the boundary condition is • an insulated wall: ∂u1 ∂u2 − = 0, ∂y ∂x
y = 0,
(85)
or • a perfectly conducting (PEC) wall: u1 = 0, 23
y = 0.
(86)
The perfectly matched layer for the Maxwell’s wave equation (83) is r ∂u ∂ ∂u ∂ ∂u ∂ ∂2u µ ∂u ∂ T ∂u ∗ A + B + C + C + σ +σ = ∂t2 ∂t ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x r ∂ (σv) µ ∂(σw) −A +B + σα (u − q) − σ ∗ σ (u − q) , ∂x ∂y ∂u ∂v = − (σ + α) v + , ∂t ∂x ∂w ∂u = −αw + , ∂t ∂y ∂q = α (u − q) . ∂t (87) As before σ (x) ≥ 0 denotes the damping function and α ≥ 0 is the complex frequency shift. Away from the physical boundary y = 0, perfect matching is achieved if the damping function vanishes at the interface σ (0) = 0. Again, in order to achieve perfect matching on the boundary y = 0, we ∂ ∂ must apply the complex change of variables ∂x → S1x ∂x to the boundary conditions. Consider the insulated wall boundary condition (85). We apply the complex change of variables and we have ∂u ˆ1 1 ∂u ˆ2 − = 0. ∂y Sx ∂x
(88)
In order to localize the modified boundary condition (88) in time we introduce the auxiliary variables (28). Choosing auxiliary variables we have two possible modified boundary conditions in the PML, ∂u1 ∂u2 − + σv2 = 0, ∂y ∂x
(89)
∂u1 ∂u2 − + σw1 = 0. (90) ∂y ∂x Also note that in the continuous setting the modified boundary conditions, (89) and (90), are equivalent and correspond to the modified boundary condition (88) in the Laplace space. On a PEC boundary, we set the tangential electric field to zero, u1 = 0.
(91)
In the analyses below we will assume a lossless media σ ∗ = 0. In a conducting media, the conductivity σ ∗ > 0 will always yield a stabilizing effect. Theorem 5 Consider the constant coefficient PML (87) in the upper halfplane Ω = (−∞, ∞) × (0, ∞) with the coefficient matrices given by (84), and subject to the boundary condition (89), (90) or (91) at y = 0. The corresponding initial boundary value problem does not support growing modes. 24
Proof: Without loss of generality we set σ ∗ = 0. We know that the dispersion relation of the Maxwell equation (83) is a circle centered at the origin. The corresponding slowness diagram is the unit circle centered at the origin. Thus the geometric stability condition is always satisfied and the corresponding constant coefficient Cauchy problem is stable. Consider now the corresponding constant coefficients half plane problem. As before, we take Laplace transform in time and Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables, leading to the system of second order ordinary differential equations dˆ ˆ d2 u u ˆ = − k˜x2 Aˆ s2 u u + B 2 + ik˜x C + C T , dy dy
(92)
where
kx . k˜x = Sx The corresponding boundary condition is • an insulated wall: dˆ u1 (ik˜x , 0, s) − ik˜x u ˆ2 (ik˜x , 0, s) = 0, dy
(93)
u ˆ1 (ik˜x , 0, s) = 0.
(94)
or • a PEC wall:
For the system (92), we seek a solution of the form ˆ = Φeκy , u and we have s2 I + kx2 A − ikx κ C + C T − κ2 B Φ = 0.
(95)
Note that Φj (j = 1, 2), are the eigenvectors of the matrix P (s, kx , κ) = s2 I + kx2 A − ikx κ C + C T − κ2 B,
(96)
and the corresponding eigenvalues are θj . The first eigenvalue θ1 corresponds to a non-propagating degenerate mode and satisfies θ1 (s, kx , κ) := (s + σ)2 = 0. 25
Thus s = −σ, and the corresponding mode is stable. The second eigenvalue θ2 (s, kx , κ) corresponds to a physical propagating mode and leads to the general solution ˆ (ikx , y, s) = δ1 Φ(ikx , s)e−κy + δ2 Φ(ikx , s)eκy , u
(97)
where κ is given by kx , k˜x = Sx
q κ = ± s2 + k˜x2 , and the eigenfunction is " Φ=
1 ˜x −ik κ
# .
To ensure a non-trivial solution the unknown constants δ1 , δ2 must not be zero at the same time. The condition kˆ u(y)k < ∞ which is equivalent to ˆ (ik˜x , y, s) = 0, leads to δ2 = 0 and limy→∞ u ˆ (ikx , y, s) = δ1 Φ(ik˜x , s)e−κy . u
(98)
To determine the constant δ1 we use the boundary conditions (93) or (94). Consider now the insulated wall boundary condition (93) and insert (98), we have κ2 + k˜x2 1 s2 (s + α + σ)2 + 2 (s + α)2 kx2 ≡ = 0. κ κ (s + α + σ)2
(99)
All roots of (99) satisfy F1 ≡
s2 (s + α + σ)2 + 2 (s + α)2 kx2 = 0. (s + α + σ)2
(100)
Note that F1 is homogeneous of degree 2 in both s, kx . At σ ≡ 0 the roots of the numerator of F1 are √ s = −α, s = 0, s = ±i 2kx . By theorems 2 and 3 we have <s ≤ 0 when σ > 0, α ≥ 0. Next consider the PEC boundary condition (94). We insert the solution (98) in (94) to determine δ1 , u ˆ1 (ik˜x , 0, s) = δ1 = 0.
(101)
Thus the PEC wall supports only trivial solutions. This completes the proof. 26
7
Conclusions and outlook
In this work we have developed a method to analyze the PML when physical boundary conditions are important. Our method is based on a mode analysis on a half plane. We apply our method to equations of linear elastodynamics subject to free–surface boundary condition and homogeneous Dirichlet boundary condition, and to curl–curl Maxwell’s equations subject to an insulated wall and a PEC wall. For the free–surface boundary condition of elastodynamics and insulated wall of electromagnetism, we propose new modified boundary conditions. Finally, we proved that these half plane problems do not support growth when the geometric stability is not violated. In a coming paper we will evaluate numerically the modified boundary conditions proposed in this paper. We will also consider waveguides. We will point out some numerical and practical issues.
Appendix In [6] the following result was proven Lemma 4 (Corollary 1 in [6]) Consider the complex variable s s2 κ= + kx2 , kx ∈ R, µ
(102)
where s = a + ib, and a, b ∈ R. If a > 0, then 0, c0 µ
c0 =κ = √ b, µ
with
0 < c0 ≤ 1.
Here we will prove the theorem Theorem 6 Consider the complex variable s s2 ˜ 2 kx + kx , k˜x = , kx ∈ R, κ= µ Sx
µ>0
(103)
where
s+α+σ , α, σ ≥ 0. s+α If a > 0, there is a constant δ > 0 such that 0. s = a + ib,
Sx =
In order to ensure compactness, we introduce the normalization (44), we have s η 02 k0 2 κ0 = + k˜0 1 , k˜10 = 10 , (104) µ Sx 27
where η 0 = a0 + ib0 , Consider Sx0 =
Sx0 =
η 0 + γ 0 + 0 , η0 + γ 0
Note that
η 0 + γ 0 + 0 . η0 + γ 0
with γ 0 , 0 ≥ 0.
1 = a0s + i0 b0s , Sx0
where a0s =
(a0 + γ 0 )2 + b02 ((a0
+
γ 0 ) (a0 b0s
+
=
γ0
+
0 )
+
b02 )2
+
a0 + γ 0
(0 b0 )2
a0 + γ 0 + 0 + b02 > 0,
(a0 + γ 0 )2 + b02 ((a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 )2 + (0 b0 )2
b0 .
Let η˜0 = η 0 Sx0 . If η 0 = a0 + ib0 , we have
0
η˜ =
a0 (a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 + 0 b02 (a0 + γ 0 )2 + b02
+ib
0
(a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 − 0 a0 (a0 + γ 0 )2 + b02 (105)
If a0 > 0, then 0
0.
(106)
Also note that when 0 ≡ 0, we have 0. Corollary 1 Consider the complex variables κ ˜ 0 = Sx0 κ0 , η˜0 = Sx0 η 0 with s κ ˜0 =
η˜02 + k 0 21 , µ
µ>0
(107)
where η 0 = a0 + ib0 ,
Sx0 =
η 0 + γ 0 + 0 , η0 + γ 0
with
|η 0 |, |k10 | ≤ 1,
and
γ 0 , 0 ≥ 0.
If a0 > 0, then 0, √
†
April 10, 2012
Abstract We study the stability of the perfectly matched layer (PML) for symmetric second order hyperbolic partial differential equations on the upper half plane, with boundary conditions at y = 0. Using a mode analysis, we develop a technique to analyze the stability of the corresponding initial boundary value problem on the half plane. We apply our technique to the PML for the elastic wave equations subject to free surface and homogenous Dirichlet boundary conditions, and to the PML for the curl–curl Maxwell’s equation subject to insulated walls and perfectly conducting walls boundary conditions. The conclusion is that these half–plane problems do not support temporally modes. Keywords: elastic waves, electromagnetic waves, boundary waves, surface waves, glancing waves, perfectly matched layer, stability, normal mode analysis.
1
Introduction
There are many wave propagation problems where boundary phenomenon such as surface or glancing waves are dominant. Typical examples are in non-destructive testing, seismology, earthquake engineering, ultrasonics, and ground penetrating radar (GPR) technologies. These problems can be described by symmetric time–dependent partial differential equations (PDE) ∗
Division of Scientific Computing, Department of Information Technology, Uppsala University Sweden † Division of Scientific Computing, Department of Information Technology, Uppsala University Sweden
1
in semi-bounded domains. In order to perform numerical simulations, the unbounded extensions of the domains must be truncated, by introducing artificial boundaries. Thus efficient and reliable domain truncation is vital, since it enables more accurate numerical simulations. Perhaps, an efficient approach to truncate unbounded domains is to surround the artificial boundaries with the perfectly matched layer (PML) [2, 17, 26, 28]. The PML consists in extending the computational domain to a layer, where the underlying equations are modified such that waves decay rapidly in the layer. A very important property of the PML is perfect matching. This means that waves propagate from the physical space into the PML without reflections. A PML is usually derived by assuming a homogeneous media and an infinite domain in all directions. In a domain with physical boundaries, the application of the PML introduces boundary corners where physical boundary conditions interact with the PML. In order to ensure accuracy, the underlying boundary conditions must be accurately extended from the interior into the PML. Thus leading to an initial boundary value problem (IBVP) for the PML. These types of problems can be found many application areas, including seismic, optical and gravitational waves. Examples include glancing waves and surface waves in electromagnetic and elastic wave propagation problems. By construction the PML is perfectly matched, but there is no guarantee that all solutions decay with time. The analysis of temporal stability of PMLs is therefore a main topic of research. For Cauchy problems, the stability of the PML can be predicted by the so–called geometric stability condition [7] for several hyperbolic systems [25, 24, 7, 12]. In several studies, numerical experiments have shown that PMLs which are Cauchy stable can exhibit (numerical) growth when some physical boundaries are imposed [9, 27, 1]. For first order hyperbolic PDEs, the theory of IBVPs is rather well–developed. The Friedrichs theory [35, 36, 37] for symmetric systems with maximally dissipative boundary conditions was available, but not sufficient for nonsymmetric systems or more general boundary conditions. A more general theory is the eigenvalue theory based on the normal mode analysis [38]. H–O. Kreiss and his co–workers developed this theory in a more general and systematic manner [33, 34, 39]. Therefore this theory is often called the Kreiss theory/technique. By the Kreiss theory, all boundary stable problems must satisfy the so–called determinant condition [33, 34, 39]. This theory has recently been extended to second order hyperbolic systems [5]. However, the PML is a more complicated system. The PML is non–standard; it is an integro–differential equation. Of course by introducing auxiliary variables the PML can be rewritten as a system of PDEs. For second order 2
systems, the PML is a combination of second order PDEs for the modified wave equation and a first order system of PDEs for the auxiliary differential equations. Unfortunately no stability theory exists for the corresponding IBVP including the PML. In this paper we consider second order symmetric hyperbolic systems on the upper half–plane with boundary conditions at y = 0. First we derive a PML that truncates the boundary in the x–direction, by a complex change of variables ∂/∂x → 1/Sx ∂/∂x. Here, Sx is the PML (complex) metric. We also transform the x–derivatives (∂/∂x) in the boundary condition with the PML metric Sx (∂/∂x → 1/Sx ∂/∂x). Choosing auxiliary variables, we construct new transformed boundary conditions in the PML. Using a mode analysis we develop a stability theory for the corresponding IBVP. A main result is that if the characteristic function for the undamped problem is a homogeneous function1 , then the introduction of the PML will immediately move all roots of the characteristic equation into the negative half complex– plane. We apply our technique to the PML for the elastic wave equation subject to free surface and homogeneous Dirichlet boundary conditions, and to the PML for the curl–curl Maxwell’s equation subject to insulated walls and perfectly conducting walls boundary conditions. The conclusion is that these half–plane problems do not support temporally growing modes. The outline of the paper is as follows. In Section 2, we introduce our model problem and discuss the stability analysis using the normal mode analysis. In Section 3, we derive the PML and the modified boundary conditions. In Section 4, we develop a technique for the stability analysis of the corresponding IBVP for the PML and formulate our main results. In sections 5 and 6, we apply our technique to the elastic wave equations and the Maxwell’s equations, respectively. In the last section we draw conclusions and suggest future work.
2
Symmetric hyperbolic systems
Consider the second order hyperbolic system of PDEs in two space dimensions ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ T ∂u A + B + C + C . (1) = ∂t2 ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x Here, the unknown u = (u1 , u2 , · · · , um )T , is a real valued vector function, (x, y) are the spatial coordinates and t denotes time. The coefficients A, B, C are real matrices and define the media. 1
Let f (v) be a function with the vector argument v. If f (αv) = αk f (v) with α 6= 0 and k ∈ Z, then f (v) is a homogeneous function of degree k.
3
Definition 1 The system (1) is called symmetric hyperbolic if the symbol Ψ (kx , ky ) = kx2 A + ky2 B + kx ky C + C T ,
(2)
is symmetric and positive semi-definite for all kx , ky ∈ R. In general the coefficient matrices A, B are symmetric and positive semi– definite. If (1) is symmetric hyperbolic the eigenvalues of the symbol Ψ(kx , ky ) are real and non–negative, and the eigenvectors are orthogonal. We set the smooth initial data u (x, y, 0) = f (x, y) ,
2.1
∂u (x, y, 0) = g (x, y) . ∂t
Plane waves and slowness diagrams
In order to understand the wave propagation properties of the model (1) it is useful to consider wave-like solutions u = u0 est−ikx x−iky y ,
u0 ∈ Rn ,
kx , ky , x, y ∈ R,
0 ≤ t.
(3)
In equation (3), (kx , ky ) ∈ R2 is the wave vector, and u0 is a vector of constant amplitude called the polarization vector. By inserting (3) in (1) we have the following solvability condition, often called the dispersion relation F0 (s, kx , ky ) ≡ det s2 I + Ψ (kx , ky ) = 0, Ψ (kx , ky ) = kx2 A + ky2 B + kx ky C + C T .
(4)
Note that the symbol Ψ(kx , ky ) is a symmetric matrix. The polarization vector u0 is the eigenvector of the matrix Ψ(kx , ky ), and the eigenvalue is −s2 . It should be noted that if (1) is symmetric hyperbolic, then the roots s of the dispersion relation (4) are distinct and purely imaginary. These roots are called the physical roots and we will refer to the corresponding modes as the physical modes. Since the physical roots are distinct and purely imaginary, and the eigenvectors u0 are orthogonal, the Cauchy problem (1) is well-posed2 . If we write s = iω(kx , ky ), where ω ∈ R is called the temporal frequency, we 2
We also note that if (1) is symmetric hyperbolic with periodic boundary conditions, it is possible to derive a strict energy estimate.
4
can also introduce kx ky , ), the normalized propagation direction, |k| |k| ω ω Vp = ( , ), the phase velocity, kx ky kx ky S = ( , ), the slowness vector, ω ω Vg = ∇k ω(kx , ky ), the group velocity, q and |k| = kx2 + ky2 . K=(
(5)
By inspection the dispersion relation (4) is homogeneous in ω, kx , ky , we can rewrite (4) as F0 (i, S) = 0.
(6)
The wave propagation properties of a certain media can be described by plotting the slowness diagrams defined by points in the S−plane satisfying (6). The group velocity can be expressed as ∂F (iω, k , k ) −1 0 x y Vg = ∇k ω(kx , ky ) = − ∇k F0 (iω, kx , ky ). ∂ω
(7)
Thus, the group velocity is normal to the slowness curve. We are now ready to introduce the geometric stability condition [7]. Definition 2 (Geometric stability condition) Let Vg = ((Vg )x , (Vg )y ) and S = (Sx , Sy ) denote the Cartesian components of the group velocity and the slowness curve respectively. The geometric stability condition in the x−direction is Sx × (Vg )x ≥ 0, for all points on the slowness curve. By Theorem 2 in [7], the geometric stability in the x–direction is necessary for stability of the vertical PML.
2.2
Boundary conditions and the normal mode analysis
Now consider (1) on the upper half–plane Ω = {(x, y) | − ∞ < x < ∞, 0 ≤ y < ∞}. At y = 0 we set the boundary condition ∂ ∂ ∂ L , , u = 0. (8) ∂t ∂x ∂y Here L is a first order partial differential operator acting only at the boundary.
5
For some boundary operators L, it is possible to establish well–posedness by deriving an energy estimate, see [13]. We will use the normal mode analysis [5, 6] to investigate well–posedness and temporal growth properties of the PML. To begin with, we introduce the L2 –norm Z ∞ 2 kv (y) k = v∗ vdy,
(9)
0
where v∗ denotes the conjugate transpose of v. Lemma 1 The problem (1) with the boundary condition (8), is not well– posed if there are non-trivial solutions of the form ˆ (y) , u = est+ikx x u
kˆ u (y) k < ∞,
0 ≤ y < ∞,
(10)
with <s > 0. Proof : If (10) is a solution then ˆ (βy) , uβ = esβt+ikx βx u
kˆ u (βy) k < ∞,
(11)
is also a solution, for any β > 0. Thus, if <s > 0, then we can always construct solutions (11) which can grow arbitrarily fast in time. Thus the problem is ill–posed. We will demonstrate how to use the normal mode analysis to reformulate Lemma 1 as an algebraic condition, the so–called determinant condition. Consider the symmetric hyperbolic system (1) in the upper half plane and the boundary condition (8) at y = 0. We take the Laplace transform in time and the Fourier transform in the tangential direction x. We obtain an eigenvalue problem ˆ = − kx2 Aˆ s2 u u + ikx C + C T
dˆ ˆ u d2 u +B 2, dy dy
d ˆ ˆ = 0, L s, ikx , u dy
y = 0,
0 ≤ y < ∞, kˆ u (y) k < ∞.
<s > 0, (12) (13)
Here (s, kx ) are the dual variables to (t, x). For the ordinary differential equation, ODE, (12) we seek solutions of the form ˆ = Φeκy , u and we have s2 I + kx2 A − ikx κ C + C T − κ2 B Φ = 0. 6
(14)
Note that Φj (s, kx , κ) are the eigenvectors of the matrix P (s, kx , κ) = s2 I + kx2 A − ikx κ C + C T − κ2 B,
(15)
and the corresponding eigenvalues θj (s, kx , κ) must satisfy θj (s, kx , κ) = 0.
(16)
Observe that if we set κ = iky , we get exactly the dispersion relation (4). The matrix P (s, kx , κ) is complex symmetric and the eigenvectors Φj (s, kx , κ) are complete. The eigenvalues θj (s, kx , κ) are quadratic functions of κj , and the complex roots κj of (16) come in pairs (+)
κj (s, kx ),
(+)
0,
(−)
0,
< 0,
where the superscript (+) denotes the root with positive real part and the superscript (−) denotes the root with negative real part. The general solution of (12) is ˆ (y) = u
m X
(+)
(+)
δj Φj eκj
y
+
j
m X
(−)
(−)
δj Φj eκj
y
.
(18)
j (+)
We consider only bounded solutions in L2 (kˆ u (y) k < ∞), then δj and we have ˆ (y) = u
m X
δj Φj eκj y ,
0.
(21)
We have Theorem 1 Consider the symmetric hyperbolic system (1) in the upper half plane and the boundary condition (8) at y = 0. The problem (1) with (8) is not well–posed if there are s with <s > 0 satisfying the determinant condition F0 (s, kx ) ≡ det Υ (s, kx ) = 0. 7
Theorem 1 is a very important result. The whole process of finding IBVPs which are not well–posed reduces to solving the algebraic equation (21). Note that most relevant physical problems are well–posed. Examples are clamped wall and free–surface boundary conditions in elastodynamics, hard wall and soft wall boundary conditions in acoustics, and insulated wall and perfectly conducting wall boundary conditions in electromagnetism. For the above examples there are energy estimates see [13], thus their determinant condition (21) is satisfied. It is important to note that using the normal mode analysis there can be a technical difficulty of establishing sharp estimates in the physical space, especially for roots corresponding to generalized eigenvalues, see [39, 5]. From now on we will assume that the system (1) is symmetric hyperbolic. We are only interested in physical problems which are well–posed. Therefore we will also assume that the problem (1) with the boundary condition (8) is well–posed. Note that if (1) with (8) is well–posed, then for all roots sj of F0 (s, kx ) = 0 we can write sj (kx ) = −aj0 (kx ) + ibj0 (kx ) , where
aj0 (kx ) ≥ 0,
aj0 (kx ) , bj0 (kx ) ∈ R,
j = 1, . . . , 2m.
(22)
Furthermore, if F0 (s, kx ) is homogeneous of degree 2m we can introduce the scaling s kx η0 = , k10 = ≡ ±1. (23) |kx | |kx | Note that all roots ηj0 of F0 (η 0 , k10 ) = 0 are complex numbers, and satisfy ηj0 = −a0j0 + ib0j0 , where
a0j0 ≥ 0,
a0j0 , b0j0 ∈ R, j = 1, . . . , 2m.
(24)
We can distinguish the following cases; • When 0, the problem (1) with (8) is called boundary stable. If there are non–homogeneous boundary data, the solutions on the boundary can be estimated in terms of the data. • When 0, 0 ≤ y < ∞. The PML will be derived using the complex coordinate stretching technique [11]. In the complex metric we will also introduce the complex frequency shift [10]. To begin, we take the Laplace transform in time Z ∞ ˆ (x, y, s) = u e−st u(x, y, t)dt. (25) 0
We then analytically continue the transformed wave equation on the complex contour, Z x σ(z) Υ(x) = x + dz, α +s 0 where oscillating solutions are turned into exponentially decaying solutions. Here, σ(x) ≥ 0 is the damping function and α ≥ 0 is the complex frequency shift (CFS). We denote the complex metric Sx :=
dΥ(x) σ(x) =1+ . dx α+s
If α = 0 we call it the standard PML and if α > 0 we call it the CFS-PML. Changing the coordinate x → Υ(x), we have ˆ ˆ ˆ ˆ 1 ∂ 1 ∂u ∂u ∂ ∂u 1 T ∂u ˆ= s2 u A +C + B + C . (26) Sx ∂x Sx ∂x ∂y ∂y ∂y Sx ∂x Multiplying (26) by Sx leads to ˆ ˆ ˆ ˆ ∂ 1 ∂u ∂u ∂ ∂u 2 T ∂u ˆ= s Sx u A +C + Sx B +C . ∂x Sx ∂x ∂y ∂y ∂y ∂x
(27)
Observe that in (26) the flux in the y-direction changes from B and from B
ˆ ˆ ∂u ∂u + CT ∂y ∂x
ˆ ˆ ∂u 1 T ∂u + C ∂y Sx ∂x
to B
ˆ ˆ ∂u 1 T ∂u + C , ∂y Sx ∂x
to Sx B
ˆ ˆ ∂u ∂u + CT , ∂y ∂x
in (27). Choosing the auxiliary variables ˆ= v
ˆ 1 ∂u , (s + α) Sx ∂x
ˆ = w 9
ˆ 1 ∂u , α + s ∂y
ˆ= q
αˆ u , α+s
(28)
and inverting the Laplace transforms yielding ∂2u ∂u ∂ ∂u ∂u ∂ ∂u T ∂u A + B +σ = +C +C ∂t2 ∂t ∂x ∂x ∂y ∂y ∂y ∂x ∂ (σv) ∂(σw) −A +B + σα (u − q) , ∂x ∂y ∂v ∂u = − (σ + α) v + , ∂t ∂x ∂w ∂u = − αw + , ∂t ∂y ∂q =α (u − q) . ∂t
(29)
Away from the boundary at y = 0, the PML model (29) is perfectly matched to the system (1) if σ (0) = 0. There are no reflections as waves pass the interface x = 0. That is, waves of all frequencies and angles are transmitted from the interior into the layer without reflections. The PML damps waves so that they become insignificant before they reach the outer boundary of the computational domain. We initialize the fields in the PML with homogeneous initial data u(x, y, 0) = 0, v(x, y, 0) = 0,
3.1
∂u (x, y, 0) = 0, ∂t
w(x, y, 0) = 0,
q(x, y, 0) = 0.
Boundary conditions
The accuracy of the PML is ensured by the perfect matching property. This can be shown mathematically by constructing the general solutions to the wave equation and the PML equation in the Laplace-Fourier space, see [14]. For a PML that depends on the x–direction, the analysis in [14] assumes that in the y-direction the domain is unbounded (or periodic) such that a Fourier transform in the y-direction is possible. If the underlying problem is bounded in one or more directions with non-periodic boundary conditions, a different technique is needed. To ensure perfect matching, the underlying boundary conditions must be simultaneously transformed using the PML metric Sx . Now consider the boundary condition (8) and introduce the Laplace transform (25). Still in the Laplace space we introduce the complex metric Sx , having 1 ∂ ∂ ˆ ˆ = 0. L s, u (30) , Sx ∂x ∂y In order to localize the transformed boundary condition (30) in time we introduce the auxiliary variables (28) and invert the Laplace transform, to 10
get L
∂ ∂ ∂ , , , σ, α (u, v, w, q) = 0. ∂t ∂x ∂y
(31)
Note that (31) is a general representation of the transformed boundary condition. Depending on the boundary condition not all the auxiliary variables v, w, q, will be required. Also note that the boundary condition (31) is a lower order modification of the boundary condition (8). If the damping vanishes, σ(x) ≡ 0, the auxiliary variables decouple completely and we recover the physical boundary condition (8).
4
Stability analysis
In order for the PML to be computationally useful, the PML must not support growth. Any growth in the layer can propagate into the computational domain and pollute the solution. The question of stability is not trivial, it is therefore a topic of active research. The perfectly matched layer is indeed a variable coefficient problem. In addition, boundary conditions can introduce many more challenging difficulties. However, we can localize the problem and split the problem into a Cauchy problem and a half-plane problem with constant coefficients. For the Cauchy problem, growth can be predicted by the so called geometric stability condition [7]. The boundary conditions can then be analyzed by using a mode analysis. If the geometric stability condition is not violated and the roots of the characteristics equations corresponding to the boundary conditions lie on the left half–complex plane, we conclude that the PML does not support growing solutions. If the geometric stability condition is violated and the roots of the characteristics equations corresponding to the boundary conditions lie on the left half–complex plane, we say that the PML supports bounded growth. At least the boundary conditions can not cause additional growth, and therefore are not harmful to the PML. In any case if the roots of the characteristics equations corresponding to the boundary conditions lie on the right half–complex plane, the boundary condition supports exponentially growing modes, and can be harmful to the PML. In this case, we say that the boundary condition is unstable.
4.1
Cauchy problem
To keep this paper self–contained we summarize the arguments leading to the geometric stability condition. More elaborate discussions can be found in [7, 12, 29]. To begin with, we freeze all coefficients and consider the Cauchy problem (29). We then assume solutions of the form V = V0 est−ikx x−iky y , 11
V = (u, v, w, q)T .
(32)
Inserting (32) in (29), we have the dispersion relation F (η, kx , ky , σ, α) ≡ (s + α)m F0 (s (s + σ + α) , (s + α) kx , ky (s + σ + α)) = 0. (33) For compactness, we introduce the normalized variables, η=
s , |k|
k1 =
kx , |k|
k2 =
ky , |k|
=
σ1 , |k|
γ=
α , |k|
|k| =
q kx2 + ky2 , (34)
and we have F (η, k1 , k2 , , γ) ≡ (η + γ)m F0 (η (η + + γ) , (η + γ) k1 , k2 (η + + γ)) = 0, (35) where F0 (η0 , k1 , k1 ) is given by (4). The system (29) is stable if for any (k1 , k2 ) satisfying k12 + k22 = 1, the roots η(k1 , k2 , , γ) have negative real parts. The existence of solutions of roots η with positive real parts would correspond to plane wave solutions with exponential growth in time. A stable system does not admit such solutions. We will use standard perturbation arguments. By setting = 0 in (35) we have F (η, k1 , k2 , 0, γ) = (η + γ)3m F0 (η, k1 , k2 ) . (36) Clearly (36) has 5m roots, 2m of which are the purely imaginary roots of (4) and correspond to the physical modes. The remaining 3m roots are real and non-positive η = −γ ≤ 0. (37) We call the corresponding modes non-physical. Since the roots depend continuously on the coefficients, it is apparent that at intermediate frequencies and with small enough damping, the roots η of the non-physical modes will remain in the left half of the complex plane. The following result was proven in [12]. Lemma 2 (Lemma 2 in [12]) At sufficiently high frequencies, the parameter γ > 0 will stabilize the non-physical modes for any ≥ 0. We see that if γ > 0, then the instability in the PML at constant coefficients can not come from the non-physical modes. The physical modes are simple (distinct), and can therefore be expanded in powers of , η = η0 + η + O 2 . (38) Here, η0 is a purely imaginary root of F0 (η, k1 , k2 ). At sufficiently high frequencies, 0. If the geometric stability condition in the x-direction is violated, then there are unstable physical modes at sufficiently high frequencies. We see that away from the boundaries, if the geometric stability condition is not violated the PML does not support growth. In order to understand the effect of the boundary conditions we perform the modal analysis as in Section 2.2.
4.2
Half plane problem
Now consider the constant coefficients time-dependent PML (29) in the upper half plane with the modified boundary condition (31) at y = 0. We take the Laplace transform in time and the Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables and we have dˆ ˆ u d2 u ˆ = − k˜x2 Aˆ s2 u u + ik˜x C + C T + B 2 , 0 < y < ∞, dy dy ∂ ˆ = 0, y = 0, u Lˆ s, ik˜x , ∂y where
<s > 0, (39)
(40)
kx k˜x = . Sx
Here (s, kx ) are the dual variables to (t, x). Introducing the modified wave– number k˜x simplifies the remaining analysis. Note that equations (39) and (40) are completely analogous to equations (12) and (13). Thus the remaining part of the analysis can be adapted from the analysis in Section 2.2. Repeating all the steps leading to theorem 1, we arrive at the determinant condition F0 s, k˜x ≡ det Υ s, k˜x 6= 0, <s > 0. (41) The stability of the boundary condition (31) at y = 0 can be formulated as an algebraic condition, Definition 3 Consider the PML (29) in the upper half plane subject to the modified boundary condition (31) at y = 0. The problem (29) with (31) is not well–posed if for any σ ≥ 0 there is a root s with <s > 0 satisfying the determinant condition s+α ˜ kx = 0. (42) F0 s, kx ≡ F0 s, s+α+σ 13
We can prove the following Theorem 2 Consider the standard PML (29) (with α ≡ 0), and the modified boundary condition (31) at y = 0. Assume that the undamped problem is well–posed, and consequently all roots s0j (kx ) of F0 (s0 , kx ) = 0 satisfy <s0j (kx ) ≤ 0 If the characteristic polynomial F0 (s0 , kx ) is homogeneous of degree 2m in both s0 and kx , then the roots of (42) are s = 0, sj (kx ) = −σ + s0j (kx ) ,
j = 1, . . . , 2m.
Proof : Introduce α = 0 in (42), ince F0 s, k˜x is homogeneous we write F0 s,
s kx s+σ
≡
s s+σ
2m F (s + σ, kx ) = 0.
Therefore, we have a multiple root s = 0 and F0 (s + σ, kx ) = 0. The corresponding roots of F0 (s + σ, kx ) = 0 are sj (kx ) = −σ + s0j (kx ) ,
j = 1, . . . , 2m.
Consider now the CFS–PML with α > 0. Note that if F0 (s, kx ) is homogeneous of degree 2m in both s and kx , we have 2m s+α s+α s+α+σ F0 s, kx ≡ F0 s , kx = 0. (43) s+α+σ s+α+σ s+α In order to ensure compactness, we introduce the normalized variables, η0 =
s , |kx |
we obtain
k10 =
η0 + γ 0 0 η + γ 0 + 0
kx ≡ ±1, |kx |
2m
F0 η
0η
0 = 0
σ1 , |kx |
+ γ 0 + 0 0 , k1 η0 + γ 0
γ0 =
α , |kx |
(44)
Note that since η 0 = −γ 0 does not satisfy (42) we have 2m η0 + γ 0 6= 0. η 0 + γ 0 + 0
= 0.
(45)
(46)
Therefore we consider 0 0 0 0η + γ + 0 F0 η , k1 = 0. η0 + γ 0
We have the result 14
(47)
Theorem 3 Consider the CFS–PML (29) (with α > 0) subject to the modified boundary condition (31) at y = 0, and the characteristics polynomial (47). Assume that the undamped problem is well–posed, that is all roots s0j of F0 (s0 , kx ) satisfy s0j = −aj0 + ibj0 ,
aj0 , bj0 ∈ R,
where
aj0 ≥ 0,
j = 1, . . . , 2m.
If the characteristic polynomial F0 (s0 , kx ) is homogeneous of degree 2m in 0 both s0 and kx , then for each root s0j of F0 (s0 , kx ), there are two roots, η1j 0 0 0 and η2j , of (47). The two roots, η1j and η2j satisfy 1) When a00j = b00j = 0, 0 = 0, η 0 = − (γ 0 + 0 ) . η1j 2j
2) When b00j = 0, 0 η1j
0
0
=− γ + +
a00j
/2 −
0 η2j = − γ 0 + 0 + a00j /2 +
r r
γ 0 + 0 + a00j
2
/4 − a00j γ 0 ,
γ 0 + 0 + a00j
2
/4 − a00j γ 0 .
3) When a00j = 0, b00j 6= 0, 0 < 0, 0 the corresponding roots have non–positive real parts. Since the roots are continuous functions of 0 > 0 there will be roots with positive real parts only if there is a purely imaginary root of (50) for some 0 > 0. We will show that this is impossible. Assume that ηj0 = iξj0 with ξj0 ∈ R, is a root for some b00j 6= 0, γ 0 , 0 > 0. Inserting ηj0 = iξj0 in (50) and considering the real and the imaginary parts separately we have −ξj02 + a00j γ 0 + b00j ξj0 = 0, γ 0 + 0 + a00j ξj0 − b00j γ 0 = 0. The two relations above are simultaneously satisfied only if γ0 = 0
or a00j = b00j = 0
or γ 0 = −(a00j
or γ 0 = 0, 0 = −a00j s a00j + 0 0 0 + ) ± ib0j . a00j
or γ 0 = −0 − a00j , b00j = 0,
Since b00j 6= 0, and γ 0 , 0 > 0, are purely real and positive these are contradictions. The proof of the theorem is complete.
This is a very useful result. If we can show that the characteristic polynomial (21) of the undamped problem is a homogeneous function we can immediately conclude that the corresponding boundary condition is stable. The introduction of the PML will move all non–zero roots, s, into the negative half complex–plane. If there are zero roots they will not be perturbed by the PML. Thus the boundary condition (31) can not be harmful to the 17
PML (29). We conclude with the following remark. Consider the half–plane PML problem (29) with (31). If the undamped problem and sj denotes is well–posed the roots of the characteristics equation F0 s, k˜x = 0, we can distinguish the following cases • no growth: The geometric stability condition is satisfied and <sj ≤ 0. • bounded growth: The geometric stability condition is violated and <sj ≤ 0. • ill–posedness: <sj > 0.
5
Linear elastodynamics
In this section, we will apply the results of Section 4 to the PML for the elastic wave equation, subject to a free–surface and a homogeneous Dirichlet boundary condition. Consider the elastic wave equation for a source free heterogeneous media in two space dimensions ρ
∂2u =∇·T, ∂t2
x ∈ Ω,
t > 0.
(57)
Here, x = (x, y) are the spatial variables, t denotes time, Ω is the upper half plane and ρ is the density of the media. The unknown u = (u1 , u2 )T , is the displacement vector and T is the stress tensor. The x and y components of the stress tensor T are ∂u ∂u ∂u ∂u Tx = A +C , Ty = B + CT . (58) ∂x ∂y ∂y ∂x The matrices A, B, C are given by the elastic coefficients. If time is scaled appropriately, then the system (57) can be rewritten as ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ T ∂u A + B + C + C . (59) = ∂t2 ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x Note that (59) is completely analogous to the model (1). In an orthotropic media whose principal axes coincide with the (x, y) axes, the fourth-order symmetric and positive definite tensor C of elastic coefficients [32] is
c11 c12 0 C = c12 c22 0 , 0 0 c33
where
c11 , c12 , c22 , c33 > 0,
18
c11 c22 − c212 > 0.
Thus, the coefficient matrices in (59) are c11 0 c33 0 A= , B= , 0 c33 0 c22
C=
0 c12 . c33 0
(60)
For the particular case of isotropic media the elastic coefficients depend on the Lam´e parameters λ, µ, c11 = c22 = 2µ + λ,
c33 = µ,
c12 = λ.
At the boundary y = 0, the boundary condition is • a free–surface: Ty := B
∂u ∂u + CT = 0, ∂y ∂x
y = 0,
(61)
or • a clamped wall: u = 0,
y = 0.
(62)
The perfectly matched layer for the elastic wave equation (59) is the same as (29). Away from the physical boundary y = 0, perfect matching is achieved if the damping function vanishes at the interface σ (0) = 0. In order to achieve perfect matching on the boundary at y = 0, we must apply the complex change of variables ∂ 1 ∂ → ∂x Sx ∂x
(63)
to the boundary conditions. Consider the free-surface boundary condition (61) in the y-direction, Ty = 0. By taking Laplace transform in time and applying the complex change of variables (63) in (61) we have B
ˆ ˆ ∂u 1 T ∂u + C = 0. ∂y Sx ∂x
(64)
In order to localize the modified boundary condition (64) in time we introduce the auxiliary variables (28). Depending on our particular choice of auxiliary variables (28) we have two possible modified boundary conditions in the PML, ∂u ∂u B + CT − σC T v = 0, (65) ∂y ∂x B
∂u ∂u + CT + σBw = 0. ∂y ∂x 19
(66)
Note that in the continuous setting the modified boundary conditions, (65) and (66), are equivalent and correspond to the modified boundary condition (64) in the Laplace space. Observe also that if the damping vanishes, σ = 0, we recover the traction–free boundary condition (61). On a clamped boundary, we set the displacement field to zero, u = 0.
(67)
Theorem 4 Consider the constant coefficient PML (29) for isotropic elastic media in the upper half-plane Ω = (−∞, ∞) × (0, ∞), and subject to the modified boundary condition (65) or (66) at y = 0, or the homogeneous Dirichlet condition (67) at y = 0. The corresponding initial boundary value problem does not support growing modes. Proof: We know that in isotropic media the geometric stability condition is always satisfied, thus by lemma 3 the corresponding constant coefficients Cauchy problem is stable. Now consider the corresponding constant coefficient half plane problem. As before, we take Laplace transform in time and Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables, leading to the system of second order ordinary differential equations dˆ ˆ d2 u u ˆ = − k˜x2 Aˆ s2 u u + B 2 + ik˜x C + C T , dy dy
(68)
where
kx k˜x = . Sx The corresponding boundary condition is • a free–surface: B
dˆ u(ik˜x , 0, s) ˆ (ik˜x , 0, s) = 0, + ik˜x C T u dy
(69)
or • a clamped wall: ˆ (ik˜x , 0, s) = 0. u
(70)
The ordinary differential equations (68) has the general solution ˆ (ik˜x , y, s) = δ1 Φ1 (ik˜x , s)e−κ1 y + δ2 Φ2 (ik˜x , s)e−κ2 y u + δ3 Φ1 (ik˜x , s)eκ1 y + δ4 Φ2 (ik˜x , s)eκ2 y . 20
(71)
For the particular case of an isotropic media κ1 , κ2 are given by s s 2 2 ˜ kx s + µkx s2 + (λ + 2µ)k˜x2 , , κ2 = ± , k˜x = κ1 = ± µ 2µ + λ Sx √ with the square root z → z being defined as the one with non–negative real part. In the appendix we have proven that if <s > 0, then 0, α ≥ 0. Next we consider the homogeneous Dirichlet boundary condition (70). We insert the solution (72) in (70) to determine, δ1 , δ2 . ˆ (ik˜x , 0, s) = δ1 Φ1 (ik˜x , s) + δ2 Φ2 (ik˜x , s), u
(78)
The corresponding linear system is 1 ˜
− iκk1x
˜x ik κ2
! δ1 = 0. δ2 1
In order to ensure a non-trivial solution we must have ! ˜x ik 1 κ2 det = 0. ˜ − iκk1x 1
(79)
(80)
Computing the determinant (80) we have s2 s2 (s + α + σ)2 + (3µ + λ) (s + α)2 kx2 F2 (s, k˜x ) ≡ k˜x2 − κ2 κ1 ≡ = 0. µ (2µ + λ) (s + α + σ)2 (81) 22
Again the characteristic equation F2 is homogeneous of degree 4 in both s, kx . At σ ≡ 0 the roots of the numerator of F2 are p s = 0, s = −α, s = ±i (3µ + λ)kx . By theorems 2 and 3 we have <s ≤ 0 when σ > 0, α ≥ 0. This completes the proof. Note that for anisotropic elastic media the expressions we obtain are in general more complicated. However, the characteristic polynomial F(s, kx ) is homogeneous in both s and kx . Thus by theorems 2 and 3, when the geometric stability condition is satisfied, the corresponding PML subject to the free–surface boundary condition (66) or the homogenous Dirichlet boundary condition (67) does not support growth.
6
Maxwell’s equations
Consider the curl–curl Maxwell’s equations in a linear, source free, homogeneous isotropic media 1 ∂2u ∂u + µσ ∗ = −∇ × ∇ × u, 2 2 c ∂t ∂t
1 c= √ . µ
(82)
Here the vector u = (u1 , u2 )T is the electric field, c is the wave speed, µ, and σ ∗ are the relative magnetic permeability, the relative electric permittivity and the conductivity of the media respectively. If time is scaled appropriately we can write r ∂2u ∂ ∂u ∂ ∂u ∂ ∂u ∂ µ ∂u ∗ T ∂u +σ = A + B + C + C . ∂t2 ∂t ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x (83) where the coefficient matrices are given by 1 0 0 0 A= , B= , 0 0 0 1
C=
0 0 . −1 0
(84)
At the boundary y = 0, the boundary condition is • an insulated wall: ∂u1 ∂u2 − = 0, ∂y ∂x
y = 0,
(85)
or • a perfectly conducting (PEC) wall: u1 = 0, 23
y = 0.
(86)
The perfectly matched layer for the Maxwell’s wave equation (83) is r ∂u ∂ ∂u ∂ ∂u ∂ ∂2u µ ∂u ∂ T ∂u ∗ A + B + C + C + σ +σ = ∂t2 ∂t ∂x ∂x ∂y ∂y ∂x ∂y ∂y ∂x r ∂ (σv) µ ∂(σw) −A +B + σα (u − q) − σ ∗ σ (u − q) , ∂x ∂y ∂u ∂v = − (σ + α) v + , ∂t ∂x ∂w ∂u = −αw + , ∂t ∂y ∂q = α (u − q) . ∂t (87) As before σ (x) ≥ 0 denotes the damping function and α ≥ 0 is the complex frequency shift. Away from the physical boundary y = 0, perfect matching is achieved if the damping function vanishes at the interface σ (0) = 0. Again, in order to achieve perfect matching on the boundary y = 0, we ∂ ∂ must apply the complex change of variables ∂x → S1x ∂x to the boundary conditions. Consider the insulated wall boundary condition (85). We apply the complex change of variables and we have ∂u ˆ1 1 ∂u ˆ2 − = 0. ∂y Sx ∂x
(88)
In order to localize the modified boundary condition (88) in time we introduce the auxiliary variables (28). Choosing auxiliary variables we have two possible modified boundary conditions in the PML, ∂u1 ∂u2 − + σv2 = 0, ∂y ∂x
(89)
∂u1 ∂u2 − + σw1 = 0. (90) ∂y ∂x Also note that in the continuous setting the modified boundary conditions, (89) and (90), are equivalent and correspond to the modified boundary condition (88) in the Laplace space. On a PEC boundary, we set the tangential electric field to zero, u1 = 0.
(91)
In the analyses below we will assume a lossless media σ ∗ = 0. In a conducting media, the conductivity σ ∗ > 0 will always yield a stabilizing effect. Theorem 5 Consider the constant coefficient PML (87) in the upper halfplane Ω = (−∞, ∞) × (0, ∞) with the coefficient matrices given by (84), and subject to the boundary condition (89), (90) or (91) at y = 0. The corresponding initial boundary value problem does not support growing modes. 24
Proof: Without loss of generality we set σ ∗ = 0. We know that the dispersion relation of the Maxwell equation (83) is a circle centered at the origin. The corresponding slowness diagram is the unit circle centered at the origin. Thus the geometric stability condition is always satisfied and the corresponding constant coefficient Cauchy problem is stable. Consider now the corresponding constant coefficients half plane problem. As before, we take Laplace transform in time and Fourier transform in the tangential direction x. We can eliminate all the auxiliary variables, leading to the system of second order ordinary differential equations dˆ ˆ d2 u u ˆ = − k˜x2 Aˆ s2 u u + B 2 + ik˜x C + C T , dy dy
(92)
where
kx . k˜x = Sx The corresponding boundary condition is • an insulated wall: dˆ u1 (ik˜x , 0, s) − ik˜x u ˆ2 (ik˜x , 0, s) = 0, dy
(93)
u ˆ1 (ik˜x , 0, s) = 0.
(94)
or • a PEC wall:
For the system (92), we seek a solution of the form ˆ = Φeκy , u and we have s2 I + kx2 A − ikx κ C + C T − κ2 B Φ = 0.
(95)
Note that Φj (j = 1, 2), are the eigenvectors of the matrix P (s, kx , κ) = s2 I + kx2 A − ikx κ C + C T − κ2 B,
(96)
and the corresponding eigenvalues are θj . The first eigenvalue θ1 corresponds to a non-propagating degenerate mode and satisfies θ1 (s, kx , κ) := (s + σ)2 = 0. 25
Thus s = −σ, and the corresponding mode is stable. The second eigenvalue θ2 (s, kx , κ) corresponds to a physical propagating mode and leads to the general solution ˆ (ikx , y, s) = δ1 Φ(ikx , s)e−κy + δ2 Φ(ikx , s)eκy , u
(97)
where κ is given by kx , k˜x = Sx
q κ = ± s2 + k˜x2 , and the eigenfunction is " Φ=
1 ˜x −ik κ
# .
To ensure a non-trivial solution the unknown constants δ1 , δ2 must not be zero at the same time. The condition kˆ u(y)k < ∞ which is equivalent to ˆ (ik˜x , y, s) = 0, leads to δ2 = 0 and limy→∞ u ˆ (ikx , y, s) = δ1 Φ(ik˜x , s)e−κy . u
(98)
To determine the constant δ1 we use the boundary conditions (93) or (94). Consider now the insulated wall boundary condition (93) and insert (98), we have κ2 + k˜x2 1 s2 (s + α + σ)2 + 2 (s + α)2 kx2 ≡ = 0. κ κ (s + α + σ)2
(99)
All roots of (99) satisfy F1 ≡
s2 (s + α + σ)2 + 2 (s + α)2 kx2 = 0. (s + α + σ)2
(100)
Note that F1 is homogeneous of degree 2 in both s, kx . At σ ≡ 0 the roots of the numerator of F1 are √ s = −α, s = 0, s = ±i 2kx . By theorems 2 and 3 we have <s ≤ 0 when σ > 0, α ≥ 0. Next consider the PEC boundary condition (94). We insert the solution (98) in (94) to determine δ1 , u ˆ1 (ik˜x , 0, s) = δ1 = 0.
(101)
Thus the PEC wall supports only trivial solutions. This completes the proof. 26
7
Conclusions and outlook
In this work we have developed a method to analyze the PML when physical boundary conditions are important. Our method is based on a mode analysis on a half plane. We apply our method to equations of linear elastodynamics subject to free–surface boundary condition and homogeneous Dirichlet boundary condition, and to curl–curl Maxwell’s equations subject to an insulated wall and a PEC wall. For the free–surface boundary condition of elastodynamics and insulated wall of electromagnetism, we propose new modified boundary conditions. Finally, we proved that these half plane problems do not support growth when the geometric stability is not violated. In a coming paper we will evaluate numerically the modified boundary conditions proposed in this paper. We will also consider waveguides. We will point out some numerical and practical issues.
Appendix In [6] the following result was proven Lemma 4 (Corollary 1 in [6]) Consider the complex variable s s2 κ= + kx2 , kx ∈ R, µ
(102)
where s = a + ib, and a, b ∈ R. If a > 0, then 0, c0 µ
c0 =κ = √ b, µ
with
0 < c0 ≤ 1.
Here we will prove the theorem Theorem 6 Consider the complex variable s s2 ˜ 2 kx + kx , k˜x = , kx ∈ R, κ= µ Sx
µ>0
(103)
where
s+α+σ , α, σ ≥ 0. s+α If a > 0, there is a constant δ > 0 such that 0. s = a + ib,
Sx =
In order to ensure compactness, we introduce the normalization (44), we have s η 02 k0 2 κ0 = + k˜0 1 , k˜10 = 10 , (104) µ Sx 27
where η 0 = a0 + ib0 , Consider Sx0 =
Sx0 =
η 0 + γ 0 + 0 , η0 + γ 0
Note that
η 0 + γ 0 + 0 . η0 + γ 0
with γ 0 , 0 ≥ 0.
1 = a0s + i0 b0s , Sx0
where a0s =
(a0 + γ 0 )2 + b02 ((a0
+
γ 0 ) (a0 b0s
+
=
γ0
+
0 )
+
b02 )2
+
a0 + γ 0
(0 b0 )2
a0 + γ 0 + 0 + b02 > 0,
(a0 + γ 0 )2 + b02 ((a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 )2 + (0 b0 )2
b0 .
Let η˜0 = η 0 Sx0 . If η 0 = a0 + ib0 , we have
0
η˜ =
a0 (a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 + 0 b02 (a0 + γ 0 )2 + b02
+ib
0
(a0 + γ 0 ) (a0 + γ 0 + 0 ) + b02 − 0 a0 (a0 + γ 0 )2 + b02 (105)
If a0 > 0, then 0
0.
(106)
Also note that when 0 ≡ 0, we have 0. Corollary 1 Consider the complex variables κ ˜ 0 = Sx0 κ0 , η˜0 = Sx0 η 0 with s κ ˜0 =
η˜02 + k 0 21 , µ
µ>0
(107)
where η 0 = a0 + ib0 ,
Sx0 =
η 0 + γ 0 + 0 , η0 + γ 0
with
|η 0 |, |k10 | ≤ 1,
and
γ 0 , 0 ≥ 0.
If a0 > 0, then 0, √
