Applications of the perfectly matched layers in a discontinuous fluid media

(1)

IT 12 023

Examensarbete 30 hp Juni 2012

Applications of the perfectly

matched layers in a discontinuous fluid media

Ali Dorostkar

Institutionen för informationsteknologi

(2)

(3)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:

Box 536 751 21 Uppsala Telefon:

018 – 471 30 03 Telefax:

018 – 471 30 00 Hemsida:

http://www.teknat.uu.se/student

Abstract

Applications of the perfectly matched layers in a discontinuous fluid media

Ali Dorostkar

In this thesis we study the applications of the PML in a multi-layered media. The PML is constructed for the scalar wave equation and the convergence and stability of the continuous problem is studied using the normal mode analysis. A high order accurate semi-discrete problem is constructed by approximating the spatial derivatives with high order finite difference operators satisfying the summation-by-parts properties.

To have a stable semi-discrete approximation of the problem, we impose boundary conditions as well as interface conditions using the simultaneous approximation term technique. In order to gain accuracy, a transformed interface condition is constructed for the PML. The semi-discrete problem is approximated using second order accurate central difference scheme. To achieve higher order accuracy we modify the time marching scheme to eliminate truncation errors. Numerical experiments are presented showing that using the proposed transformed interface conditions, higher order of accuracy and convergence are achieved.

Examinator: Jarmo Rantakokko Ämnesgranskare: Gunilla Kreiss Handledare: Kenneth Duru

(4)

(5)

1 Introduction

In this thesis, we consider the perfectly matched layers (PML) for wave equations in a discontinuous media. Wave propagation problems are mostly formulated in an unbounded media. When numerical methods such as finite di↵erence methods are used to solve these problems, the domain must be truncated to a bounded domain with artificial boundary conditions. One major problem to this approach is that numerical reflections can be introduced into the solution. A suitable approach for domain truncation is the perfectly matched layer (PML) [1, 2, 6, 12]. The PML is an absorbing layer placed around the computational domain which absorbs outgoing waves.

The perfect matching property implies that no reflection is introduced into the computational domain as the waves penetrate the PML.

An efficient way to numerically solve wave propagation problems is to use a high-order accurate spatial scheme with a high-order accurate time marching method [7, 10]. Also to have an accurate numerical method for initial boundary value problems (IBVP), the boundary conditions must be accurately imposed. The SBP-SAT scheme is a stable technique to couple boundary conditions when using finite di↵erence methods [3, 7, 9, 10, 11]. The SBP operators are finite di↵erence operators that satisfy the summation-by-parts (SBP) properties. The simultaneous approximation term (SAT) is a method to impose boundary conditions weakly using penalty terms [3, 11]. The SAT method is also used for numerical interface treatment for wave propagation in discontinuous medias [11].

In theory the PML is a continuous unbounded layer surrounding the computational domain. However, when numerical methods are used, the PML changes to a discrete, finite width layer. In the discrete setting, in order to ensure accuracy and discrete stability, extra care should be taken when using the PML. For instance, if there are physical boundaries extending into the PML, the underlying boundary conditions must also be extended from the physical domain into the PML [4].

The PML and other artificial boundary conditions are often constructed by assuming a homogeneous and infinite media. However, real media are heterogeneous or discontinuous. For example in underwater acoustics it is relevant to consider waveguides consisting of several layers such as air, water, soft and hard sediment and bedrock layers. Applications arising in geophysics and electromagnetic problems can be composed of layers of rock, water and possibly oil.

(8)

The ultimate goal of the project is to investigate the efficiency of the PML in a layered fluid media. We are particularly interested in computational set-ups comprising of multi-block domains. The set-up consists of smaller structured domains that are patched together to a global domain using the SBP-SAT technology.

Here, we consider the scalar wave equation in a discontinuous media. We derive interface conditions and formulate a SBP-SAT approximation. Using the energy method, we derive sharp estimates of the penalty parameters, thus proving strict stability. We derive the PML and the corresponding modified interface conditions. Stability of the continuous PML is proven by a modal ansatz. In the discrete setting our modified interface conditions are crucial in order to ensure accuracy.

The outline of the paper is as follows. In section two, we introduce a 1D model problem and the needed background. In section three, we consider a wave propagation problem in a discontinuous media. We derive the interface conditions such that the continuous problem is stable. We construct the semi-discrete approximation of the problem using the SBP-SAT scheme.

We prove the discrete stability of the method by deriving sharp estimates for the penalty parameters. In section four, the PML is constructed for 2D wave equation in a discontinuous media. In the continuous setting, the stability of the PML is proven for both homogeneous and discontinuous media.

We also modify the interface conditions and extend them to the auxiliary variables. In section five, we construct stable semi-discrete approximation of the PML. We also construct the fully discrete problem for the PML. The numerical experiments are presented in section six. We use the SBP operators used in [4] for our experiments. We start by verifying the accuracy of the interface treatment. We continue to study the accuracy of the PML. We summarize the thesis in the last section.

2 Background

In this section we consider the scalar wave equation in a homogenous fluid media with homogeneous Neumann boundary conditions. We will show well- posedness and stability using the energy method. We will introduce SBP operators and derive strictly stable numerical boundary treatments using the simultaneous approximation term (SAT) method.

The scalar wave equation in a homogenous fluid media is 1

c²

@²

@t²u(x, t) = @²

@x²u(x, t), x2 (0, 1), (2.1)

(9)

subject to initial conditions

u(x, 0) = f₁(x), @

@tu(x, 0) = f₂(x). (2.2) The boundary conditions are

@

@xu(0, t) = 0, @

@xu(1, t) = 0. (2.3) Note that f₁(x) and f₂(x) are smooth functions compactly supported in (0, 1). The wave equation (2.1) can describe acoustic wave propagation in fluids. Here u is the acoustic pressure and c is the speed of sound.

2.1 Well-posedness and stability To begin, we define the energy

E(t) = 1 c

@u

@t

2

+ @u

@x

2

. (2.4)

In RHS of (2.4), the first term is the kinetic energy and the second term is the potential energy. We multiply (2.1) by @u/@t and integrate over (0, 1), we have

@

@t Z ₁

0

✓1 c

@u

@t

◆2

dx = Z ₁

0

@u

@t

@²u

@x²dx. (2.5)

Using integration by parts on the RHS of (2.5), we have

@

@t Z ₁

0

✓1 c

@u

@t

◆2

dx = @

@t Z ₁

0

✓@u

@x

◆2

dx + @u

@t

@u

@x

1 0

(2.6)

() @

@t 1 c

@u

@t

2

+ @u

@x

2!

= @u

@t

@u

@x

1 0

= 0, (2.7)

() E(t) = E(0). (2.8)

The energy is conserved. That is the solution u is uniformly bounded by the initial data for all t > 0. Thus the problem (2.1) is well-posed and stable.

2.2 Summation By Parts operators

Consider the computational domain x2 [0, 1]. Let the grid in the x coordi- nate be defined by

xj = jh, j2 0, 1, · · · N 1, N, h = 1

N 1, (2.9)

where N 2 N and the grid function is a vector v of length N.

(10)

Let D₁, D₂ be discrete operators approximating the first and the second derivatives respectively, that is

D1 ⇡ @

@x, D2 ⇡ @²

@x². (2.10)

The discrete operators D₁, D₂ are summation-by-parts (SBP) operators if they satisfy the following properties

D₁ = H ¹Q, Q^T + Q = B, B = diag( 1, 0, 0, ..., 1), (2.11)

D₂ = H ¹( M + BS), M = M^T, v^TM v 0, (2.12)

H = H^T, v^THv > 0. (2.13)

Here v is a vector of length N , Q is almost a skew-symmetric operator, M is called the symmetric part of the operator D2 and the operator S is a one sided approximation of the first derivative @/@x on the boundaries.

H defines a discrete norm, kvkH = v^THv. Note that if H = hI (where I is the identity matrix) we get the standard discrete norm. Also note that there exists SBP operators if the computational domain is unbounded but the structure of these operators defer from the description above.

Discrete operators satisfying the summation-by-parts properties (2.11), (2.12), (2.13) are usually derived using 2r order (r = 1, 2,· · · ) centered di↵erence schemes away from the boundaries and lower order schemes (usually of order r) close to the boundaries. The operator H can be a diagonal norm or a block norm. In our experiments we will use the diagonal norm SBP operators.

2.2.1 Semi-discrete approximation(SBP+SAT)

Here we will discretize the wave equation (2.1) using SBP operators and numerically impose the boundary condition (2.3) using the SAT method.

The semi-discrete approximation is strictly stable if the discrete equations satisfy discrete energy estimates similar to energy estimate (2.8) of the continuous problem. Strict stability is very important if longtime calculations are required.

Consider the IBVP problem (2.1) subject to the initial condition (2.2) and boundary conditions (2.3). The semi-discrete approximation using the SBP- SAT scheme is

1 c²

d²v

dt² = D₂v ⌧ H ¹BSv, (2.14)

(11)

where D₂ is the SBP operator (2.12), (2.13) and ⌧ is the penalty term. We define discrete energy

E(t) = 1 c

@v

@t

2 H

+ v^TM v. (2.15)

Substituting the SBP operator D₂ in (2.14) and setting ⌧ = 1 we have 1

c² d²v

dt² = H ¹M v. (2.16)

Multiplying (2.16) by (dv/dt)^TH and adding the transpose of the products we have

d

dtE = 0, (2.17)

() E(t) = E(0). (2.18)

The discrete approximation (2.14) is strictly stable.

3 Discontinuous media

In this section we will consider the wave equation in a discontinuous media.

We will derive interface conditions . We will derive a SBP-SAT semi-discrete approximation and prove stability by deriving sharp estimates of the penalty terms for SBP operators with di↵erent order of accuracy.

The discontinuity in the media is reflected on the underlying system of equations by changes (jumps) in the velocity coefficient. We consider 1D scalar wave equation in a domain consisting of two materials coupled at x = 0, the system of equations is

1 c²₁

@²u

@t² = @²u

@x² x2 ⌦1= ( 1, 0) (3.1)

1 c²₂

@²v

@t² = @²v

@x² x2 ⌦2 = (0,1) (3.2)

u(0, x) = f₁(x), @

@tu(0, x) = g₁(x) (3.3)

v(0, x) = f2(x), @

@tv(0, x) = g2(x), (3.4)

subject to the initial conditions (3.3) and (3.4). Here c₁ and c₂ are the wave speeds and f₁(x), f₂(x), g₁(x) and g₂(x) are smooth functions.

At x = 0 we need a relation between u and v in order to compute the

(12)

solution. This conditions are called physical interface conditions. We will

use ⇢

u = v,

c²₁^@u_@x = c²₂^@v_@x, (3.5) as our interface conditions. Physically, the first condition corresponds to the equality of the acoustic pressures while the second condition express that the flux of the waves should be the same on the interface.

We show that using the proposed interface conditions (3.5), the energy estimate can be derived thus proving that the problem is well-posed.

To begin, We define the energy as E = @u

@t

2

⌦1

+ c₁@u

@x

2

⌦1

+ @v

@t

2

⌦2

+ c₂@v

@x

2

⌦2

.

Multiplying (3.1) by u_t and (3.2) by v_t and integrating over the domain we have

d dt

@u

@t

2

⌦1

= Z

⌦1

2@u

@t

@²u

@t²dx = Z

⌦1

2@u

@t

@²u

@x²dx

= c²₁@u

@t

@u

@x|⁰ Z

c²₁ @²u

@t@x

@u

@xdx (integration by parts) ) d

dt

@u

@t

2

⌦1

+ c₁@u

@x

2

⌦1

!

= 2c²₁@u

@t

@u

@x|⁰. (3.6)

In the same manner for v we have ) d

dt

@v

@t

2

⌦2

+ c₂@v

@x

2

⌦2

!

= 2c²₂@v

@t

@v

@x|0. (3.7) Adding (3.6) and (3.7) together we have

d

dtE = 2c²₁@u

@t

@u

@x 2c²₂@v

@t

@v

@x, Using the interface conditions (3.5) we have

d

dtE = 0 () E(t) = E(0).

The energy is conserved and the problem is well-posed.

(13)

3.1 Semi-discrete approximation

We use the SBP-SAT scheme introduced in section 2.2 and impose interface conditions weakly. To simplify, we consider a homogeneous medium with, c₁= c₂and an artificial interface at x = 0. The semi-discrete approximation using the SBP-SAT scheme is

wtt= Dxxw ⌧NHe ¹BDe xw NHe ¹D_x^TBb^Tw+ (3.8)

⌧0He ¹Bwe ⌧1He ¹Bwe t, where w, D_xx, eB, bB, eH and D_x are defined as

w =

✓u v

◆

, B =e

0

@

0 0

Pe

0 0

1

A , P =e

✓ 1 1

1 1

◆

, (3.9)

D_xx=

✓D₂ 0 0 D₂

◆

, B =b 0

@

0 0

Pb

0 0

1

A , P =b

✓ 1 1

1 1

◆

, (3.10)

D_x=

✓D1 0 0 D₁

◆

, H =e

✓H 0

0 H

◆

. (3.11)

In (3.8) the first term in the right hand side is the discrete approximation of

@²w/@x², the second term is discrete imposition of the interface condition c²₁(@u/@x) = c²₂(@v/@x), the third and fourth term impose u = v and the last term imposes @u/@t = @v/@t.

Note that eB is zero everywhere except on the interface, that is ( eBw)i = 0, i /2 {N, N + 1}, ( eBw)_N = u_N v₀ and ( eBw)_{N +1}= v₀ u_N. Also note that 0 is a zero matrix.

3.2 Stability

We analyze the stability of the numerical interface treatment (3.8) for 2nd, 4th and 6th order accurate SBP operators. The aim is to derive sharp estimates of the penalty terms ⌧₀, ⌧_N, ⌧₁, _N such that we have an energy estimate. Note that these estimates are sharp, meaning that if values smaller than these estimates are used our semi-discrete approximation becomes un- stable. To begin with, we need some definitions and theorems.

Definition 1 (Gershgorin disc). Let A be a complex n⇥ n matrix, with entries a_ij . For 1 i  n let Ri=Pn

j=1 j6=i|aij| be the sum of the absolute values of the non-diagonal entries in the i_throw. Let D(aii, Ri) be the closed disc centered at a_ii with radius R_i. Such a disc is called a Gershgorin disc.

Theorem 1 (Gershgorin theorem). Every eigenvalue of A lies within at least one of the Gershgorin discs D(a_ii, R_i).

(14)

Proof. See [8]

Theorem 2. Let A and B be symmetric matrices, let _i(A) denote the i_th eigenvalue, where we order the eigenvalues in increasing order. If _i(B) is nonnegative and A is positive definite, then

i(BA) min(B) i(A) Proof. See [5].

We start by replacing D_xxwith the definition of the SBP operator for second derivative. We also put ⌧_N = ¹₂ and _N = ¹₂ and simplify, then (3.8) can be rewritten as

w_tt= H ¹

✓ M 0

0 M

◆ w +1

2H ¹BDb _xw + 1

2H ¹D^T_xBb^Tw

⌧₀H ¹Bwe ⌧₁H ¹Bwe _t,

Here w and D_x are the same as in (3.9) and (3.11). We introduce two matrices

Q = bBD_x ⌧₀Be (3.12)

M =f

✓ M 0

0 M

◆ 1

2 Q + Q^T . (3.13)

The semi-discrete problem becomes

w_tt = H ¹M wf ⌧₁H ¹Bwe _t. (3.14) Lemma 1. Consider the matrix fM and Q defined in (3.13) and (3.12). If Q can be written as Q = BA where the matrix ee B is defined in (3.9) and the matrix A is symmetric positive definite then fM is symmetric positive semi-definite.

Proof. We know that eB is symmetric positive semi-definite and A is symmetric positive definite, hence by Theorem 2 we have _i( eBA) 0.

We know that eBA + ( eBA)^T is symmetric and since it has positive eigenvalues, it is positive semi-definite. We conclude that Q + Q^T is negative semi-definite. From the definition of SBP operators we know that

✓M 0

0 M

◆

is symmetric positive semi-definite, therefore fM is symmetric positive semi- definite.

(15)

Lemma 2. Consider the system of ODEs (3.14). If fM is symmetric positive semi-definite and ⌧1 0 then

Ew(t) =kwtk²_H_e+ w^TM wf (3.15) is a semi-discrete energy and

E_w(t) = E_w(0) Z

2⌧₁w^T_tBwe _tdt. (3.16) Therefore the semi-discrete approximation (3.14) is stable.

Proof. We start by multiplying (3.14) by w_t^TH from left. We havee

w^T_tHwe _tt= w^T_tM wf ⌧₁w^T_tBwe _t. (3.17) Taking the transpose of (3.17) we have

w^T_ttHwe _t= w^TM wf _t ⌧₁w_t^TBwe _t. (3.18) Adding (3.17) and (3.18) we get

d E

dt = 2⌧1w^T_tBwe t. (3.19) choosing ⌧₁ 0 we see that the energy is bounded. Note that if ⌧₁ = 0 we have the equality E_w(t) = E_w(0) and the energy is conserved and when

⌧₀ > 0, the energy decays.

Next we will show that Q = BA for 2nd, 4th and 6th order accuratee SBP operators. We will also determine sharp estimates of ⌧0 such that A is symmetric positive definite. This immediately leads to stability of (3.14) by Lemma 2.

Second order accurate SBP operator We consider the second order SBP operator. Note that for this we define D_x as

S = 1 h

0 BB

@

1 1 0

0 1 1

1 CC

A , D^x =

✓S 0 0 S

◆

. (3.20)

We replace D_x in (3.12). We have

Q = bBD_x ⌧₀B =e 1 h

0 BB

@

0 0

1 1 1 1

0 0

1 CC A ⌧₀

0 BB

@

0 0

1 1

0 0

1 CC A , (3.21)

(16)

where 0 is a zero matrix. We have (Qw)N = 1

h( uN 1+ uN v0+ v1) ⌧0(uN v0) (Qw)_{N +1}= (Qw)_N

Setting ⌧₀ = ^⌧_h⁰⁰, we have (Qw)N = 1

h uN 1+ uN v0+ v1 ⌧₀⁰uN + ⌧₀⁰v0

= 1

h u_{N 1}+ ⌧₀⁰u_N+ v₀

| {z }

(Aw)N

+1

h u_N + ⌧₀⁰v₀+ v₁

| {z }

(Aw)N +1

.

and

(Qw)_{N +1}= 1

h u_{N 1}+ ⌧₀⁰u_N + v₀

| {z }

(Aw)N

1

h u_N + ⌧₀⁰v₀+ v₁

| {z }

(Aw)N +1

.

We can write Q as

Q = BA,e A = 0 BB B@

⌧₀⁰ 1 0

1 ⌧₀⁰ 1 . .. ... ...

0

1 CC CA.

The matrix A is symmetric. In order to ensure positive eigenvalues we apply the Gershgorin theorem. We have

⌧₀⁰ 2 (3.22)

)⌧0

2

h. (3.23)

From the above we conclude that using the second order accurate SBP operator, the numerical interface treatment (3.8) is stable for

⌧_N = 1

2, _N = 1

2, ⌧₀ 2

h, ⌧₁ 0. (3.24)

Fourth order accurate SBP operator The Fourth order operator Dx

is

S = 1 h

0 BB

@

d₀ d₁ d₂ d₃ d₄ · · ·

· · · d₄ d₃ d₂ d₁ d₀ 1 CC

A , (3.25)

d_i 0, d₀ d₁ d₂+ d₃ d₄= 0. (3.26)

(17)

Substituting D_x in (3.12), we have (Qw)_N = 1

h[ d₄u_{N 4}+ d₃u_{N 3} d₂u_{N 2} d₁u_{N 1}+ d₀u_N+ (3.27) d₀v₀+ d₁v₁+ d₂v₂ d₃v₃+ d₄v₄] ⌧₀(u_N v₀), and

(Qw)N +1= (Qw)N. We set ⌧₀= ^⌧⁰⁰_h^d⁴ and rearrange (3.27) to

(Qw)N = 1

h(Aw)N+ 1

h(Aw)N +1 (3.28)

(Qw)_{N +1}= 1

h(Aw)_N 1

h(Aw)_{N +1}. (3.29) where

(Aw)_N = d₄u_{N 4} (d₃ d₄)u_{N 3} (d₁ d₀)u_{N 2}+ d₀u_{N 1}+ (3.30) + (⌧₀⁰ d₄)u_N + d₀v₀+

(d₁ d₀)v₁ (d₃ d₄)v₂+ d₄v₃, and

(Aw)_{N +1}= d₄u_{N 3} (d₃ d₄)u_{N 2} (d₁ d₀)u_{N 1}+ d₀u_N+ (3.31) + (⌧₀⁰ d₄)v₀+ d₀v₁+

(d₁ d₀)v₂ (d₃ d₄)v₃+ d₄v₄. Note that one can easily reach (3.27) by substituting (3.30) and (3.31) in (3.28). The matrix A is then

A = (a_ij) = 8>

>>

><

>>

:

⌧₀⁰ d₄ i = j d₀ |i j| = 1

d1+ d0 |i j| = 2 d₃+ d₄ |i j| = 3 d₄ |i j| = 4

0 otherwise

(3.32)

Applying the Gershgorin theorem to ensure that A is positive definite we have

a_ii= ⌧₀⁰ d₄ 2 Xn j=1 j6=i

|aij|

)⌧0

2 h

Xn j=1 j6=i

|aij|

)⌧0

2

h(d₀+|d0 d₁| + |d4 d₃| + d4). (3.33)

(18)

We have a stable method if (3.33) holds. Replacing d_i with the values from SBP operators used in this thesis we have

⌧₀ 5.3332

h . (3.34)

Sixth order accurate SBP operator The sixth order operator D_x is

S = 1 h

0 BB

@

d0 d1 d2 d3 d4 d5 d6 · · ·

· · · d6 d5 d4 d3 d2 d1 d0

1 CC

A (3.35)

d_i > 0, d₀+ d₁ d₂ d₃ d₄+ d₅ d₆= 0 Substituting D_x in (3.12), we have

(Qw)_N = 1

h[d₆u_{N 6} d₅u_{N 5}+ d₄u_{N 4}+ d₃u_{N 3}+ d₂u_{N 2}+ (3.36) d₁u_{N 1}+ d₀u_N+

d0v0+ d1v1 d2v2 d3v3 d4v4+ d5v5 d6v6]+

⌧₀(u_N v₀), and

(Qw)_{N +1}= (Qw)_N. We set ⌧₀= ^⌧⁰⁰_h^d⁶ and rearrange (3.36) to (3.28) where

(Aw)_N = d₆u_{N 6}+ ( d₆+ d₅)u_{N 5}+ ( d₆+ d₅ d₄)u_{N 4}+ (3.37) +(d2 d1+ d0)uN 3+ ( d1+ d0)uN 2+ d0uN 1+ (⌧₀⁰ d6)uN+ +d₀v₀+ ( d₁+ d₀)v₁+ (d₂ d₁+ d₀)v₂+ ( d₆+ d₅ d₄)v₃+ +( d₆+ d₅)v₄ d₆v₅,

and

(Aw)_{N +1}= d6u_{N 5}+ ( d6+ d5)u_{N 4}+ ( d6+ d5 d4)u_{N 3}+ (3.38) +(d2 d1+ d0)uN 2+ ( d1+ d0)uN 1+ d0uN+ (⌧₀⁰ d6)v0+ +d₀v₁+ ( d₁+ d₀)v₂+ (d₂ d₁+ d₀)v₃+ ( d₆+ d₅ d₄)v₄+ +( d₆+ d₅)v₅ d₆v₆,

Note that one can easily reach (3.36) by substituting (3.37) and (3.38) in

(19)

(3.28) The matrix A is then

A = (a_ij) = 8>

>>

><

>>

:

⌧₀⁰ d6 i = j

d₀ |i j| = 1

d₁+ d₀ |i j| = 2 d2 d1+ d0 |i j| = 3 d₄+ d₅ d₆ |i j| = 4 d₅ d₆ |i j| = 5

d6 |i j| = 6

0 otherwise

(3.39)

Using the Gershgorin theorem we have

⌧₀⁰ d₆ 2(d₀+|d0 d₁| + |d2 d₁+ d₀| + | d₄+ d₅ d₆| + |d5 d₆| + d6) )⌧0 2

h(d0+|d0 d1| + |d2 d1+ d0| + | d4+ d5 d6| + |d5 d6| + d6).

(3.40) We have a stable 6th order accurate method if (3.40) holds. Replacing d_i with the values from the SBP operators used in this thesis we have

⌧₀ 6.9363

h . (3.41)

4 The Perfectly Matched Layer

In this section we will derive the PML for the 2D wave equation in a discontinuous media. We will derive corresponding modified material interface condition for the PML. We will also prove stability of the PML in homogeneous and discontinuous media.

Consider the wave equation in a discontinuous media 1

c²₁u_1tt= u_1xx+ u_1yy, y 2 ⌦1 = ( 1, 0), x 2 R, (4.1) 1

c²₂u_2tt= u_2xx+ u_2yy, y2 ⌦2= (0,1), x 2 R, (4.2) with c₁ 6= c2 where the physical interface is placed at y = 0 and the condi-

tions are ⇢

u₁ = u₂

c²₁u_1y = c²₂u_2y. (4.3) Assume that we are interested in computing the solution in the right half- plane 0 x < 1. To absorb out-going waves, we introduce the PML in the left half-plane. The PML is derived by analytically continuing the equations

(20)

to the complex contour. We start by taking Fourier transformation in time, we have

!²

c²₁uˆ₁ = ˆu_1xx+ ˆu_1yy, y2 ⌦1, x2 R, (4.4)

!²

c²₂uˆ₂ = ˆu_2xx+ ˆu_2yy, y2 ⌦2, x2 R. (4.5) We analytically continue (4.4) and (4.5) on the complex contour,

˜

x = x + (x)

i ! . (4.6)

Here (x) is a smooth, non-negative and increasing function. Changing The coordinates to ˜x we have

!²

c²₁uˆ₁ = 1 Sx

✓ 1 Sx

ˆ u_1x

◆

x

+ ˆu_1yy, y2 ⌦1, (4.7)

!²

c²₂uˆ₂ = 1 Sx

✓ 1 Sx

ˆ u_2x

◆

x

+ ˆu_2yy, y2 ⌦2, (4.8) where

S_x = 1 + (x)

i! , (x) = d

dx (x). (4.9)

The smooth function (x) which is zero in the interior and positive in the PML is called the damping function. Note that

S_x= 1 + (x)

i ! () 1 = 1 S_x + 1

S_xi !

() 1

S_x = 1 1

S_xi !. (4.10)

Multiplying equations (4.7) and (4.8) by Sx and substituting (4.10) we have S_x!²

c²₁uˆ₁ =

✓✓

1 1

S_xi !

◆ ˆ u_1x

◆

x

+ S_xuˆ_1yy, y2 ⌦1, (4.11) S_x!²

c²₂uˆ₂ =

✓✓

1 1

S_xi !

◆ ˆ u_2x

◆

x

+ S_xuˆ_2yy, y2 ⌦2. (4.12) Note that in (4.11) and (4.12) the flux in the y direction changes to

ˆ

u_1y ! Sxuˆ_1y, uˆ_2y! Sxuˆ_2y. (4.13) To localize in time we introduce two auxiliary variables ˆv = _S¹_x^u_{i !}^ˆ^x and ˆw =

ˆ uy

i ! and apply inverse Fourier transform to (4.11) and (4.12). We have 8>

<

>:

1

c²₁ (u_1tt+ u_1t) = u_1xx+ u_1yy ( v₁)_x+ w_1y v_1t = u_1x v₁, w_1t= u_1y

, y2 ⌦1, (4.14)

Applications of the perfectly matched layers in a discontinuous fluid media

Examensarbete 30 hp Juni 2012

Applications of the perfectly

matched layers in a discontinuous fluid media

Ali Dorostkar

Institutionen för informationsteknologi

Abstract

Applications of the perfectly matched layers in a discontinuous fluid media

Contents

1 Introduction

2 Background

3 Discontinuous media

4 The Perfectly Matched Layer