We propose a simple method for constructing non-reflecting boundary conditions via Perfectly Matched Layer approach. The basic idea of the method is to build a layer with high rate damping properties with are provided by adding the stiff relaxation source terms to all equations of the system. No complicated modification of the system to be solved is then required

Romenski, E.

Titarev, V.A.

Toro, E.F.

English

University of Trento : UNITN-eprints

PERFECTLY MATCHED LAYERS WITH HIGH RATE
DAMPING FOR HYPERBOLIC SYSTEMS
E. ROMENSKI 1, V.A. TITAREV and E. F. TORO
Laboratory of Applied Mathematics, Faculty of Engineering,
University of Trento, Italy
romenski@ing.unitn.it; titarev@mail.ru; toro@ing.unitn.it
Abstract
We propose a simple method for constructing non-reflecting boundary conditions
via Perfectly Matched Layer approach. The basic idea of the method is to build
a layer with high rate damping properties with are provided by adding the stiff
relaxation source terms to all equations of the system. No complicated modification
of the system to be solved is then required.
1 Introduction
The computational problem of wave propagation in infinite domains arises in acoustics,
seismology and electromagnetic waves phenomena. The numerical study is usually refor-
mulated in a bounded artificial computational domain in which infinity is modelled by
some properties of the numerical boundary or its neighborhood. The purpose of such a
reformulation of the problem is to avoid wave reflection from the numerical boundary.
The problem of constructing numerical boundary conditions has been intensively stud-
ied in recent decades. As a result, two main approaches have been developed. The first
approach is the so called ”Absorbing Boundary Conditions” (ABC) approach and con-
sists of formulating direct boundary conditions for the computational region, which would
eliminate wave reflection from the boundary. We refer the reader to some recent papers
[8, 9] and bibliography therein.
Another approach is the ”Perfectly Matched Layer” (PML) approach in which the
computational domain is surrounded by an additional boundary layer. The main idea
is to let the waves propagate out of the computational region. Therefore, it is crucial
that the boundary layer does not generate waves propagating back to the computational
region. Various versions of the PML approach have been proposed in the past, which are
based on introducing new artificial variables and differential equations for these variables.
See [5, 4, 2] and references therein.
In this paper we propose a simple, reliable and efficient PML method suitable for
solving wave propagation problems described by hyperbolic systems, either linear or non-
linear. The basic idea of the method is to formulate a perfectly matched layer as a layer
1On leave from Sobolev Institute of Mathematics, Russian Academy of Sciences, Novosibirsk, Russia
1
with high rate damping properties. These properties are provided by adding the stiff
relaxation source terms to all equations of the system. An asymptotic analysis for small
value of relaxation time τ allow us to make a conclusion that the amplitude of the wave
reflected from the PML is of the second order on τ . Therefore the influence of the angle
of incidence and wavelength of the wave incoming to the PML from the computational
domain on the reflected wave can be made negligible. The main advantages of the method
proposed here are its simplicity and generality. No complicated modification of the system
to be solved is required. The structure of the source terms allows one to use practically
any advection scheme without significant changes. No stability problem exists either.
A similar idea to use source terms with soft damping rate of acoustic waves in the cer-
tain and rather large buffer layer (comparable with computational domain) was proposed
in [3]. In the quoted paper the mesh stretching and filtering together with non-reflecting
boundary conditions have been used in order to provide outflow disturbance without re-
flection for aerodynamic sound generation problems. The absorbing layers with relaxation
damping have been proposed also for electromagnetic equations [1] and for linearized Euler
equations [5]. In these papers absorbing technique requires additional reformulation the
governing equations (splitting in the coordinate direction).
The rest of the paper organized as follows. In Section 2 we formulate the idea of the
new perfectly matched layer for the general case of nonlinear systems of conservation laws
and the method of its numerical solving. In Section 3 we analyze the idea as applied to
linear acoustics. In Section 4 we discuss the application of the method to two-dimensional
equations and present some numerical examples. Conclusions are drawn in Section 5.
2 Equations for perfectly matched layer
We study the wave propagation phenomena described by a hyperbolic system of conserva-
tion laws in an infinite spacial domain. Generally speaking, the system can be nonlinear.
In the three dimensional Cartesian coordinate system xi it can be written as
∂
∂t
Q+
∂
∂x
F (Q) +
∂
∂y
G(Q) +
∂
∂z
H(Q) = 0, (1)
where Q is the vector of conserved variables, F (Q), G(Q), H(Q), are the flux vectors in
the coordinate directions.
To solve the problem numerically we apply the PML strategy, which consists of defin-
ing the finite computational domain D in which we intend to obtain a solution, and then
constructing a surrounding absorbing boundary layer of a prescribed width. This bound-
ary layer must not affect the basic computational domain, meaning that there are no
waves coming back to the basic domain from the boundary layer.
Suppose that the values Q∞ of the conserved variables at infinity are known: Q→ Q∞
if xi → ∞. We formulate the equations which describe the wave propagation in the
2
perfectly matched layer by adding a relaxation source term to (1):
∂
∂t
Q+
∂
∂x
F (Q) +
∂
∂y
G(Q) +
∂
∂z
H(Q) = −1
τ
(Q−Q∞), (2)
where τ is the relaxation time which can be a function of variables Q. We take τ = ∞
inside our computational domain whereas in the PML layer τ is taken to be small.
The idea of adding such relaxation term is to provide rapid (exponential) damping
of all variables and waves inside the layer. It is clear intuitively that the decrease of the
wave amplitude will be faster if the relaxation time tends to zero. Note that the structure
and the type of equations does not change making it possible to use the same numerical
method as for the original system (1).
We now proceed to describe the numerical procedure to be used in the new PML.
Suppose that we have a one-step numerical scheme for solving (1), which gives us a result
in the form
Qn+1 −Qn
∆t
= L(Qn).
Here ∆t is the time step. The simplest numerical procedure is the one step implicit
approximation of the source term in the system (2). That is we use the following formula:
Qn+1 −Qn
∆t
= L(Qn)− 1
τ
(Qn+1 −Q∞). (3)
It gives us the result in the form(
1 +
∆t
τ
)
Qn+1 =
∆t
τ
Q∞ +Qn +∆tL(Qn)
¿From the latter formula we see that if the relaxation time τ is very small (in particular
τ much less than ∆t) then the value of Qn+1 is very close to its value at infinity Q∞.
3 Acoustic wave propagation in the PML
In this section we study the method as applied to acoustic waves leaving the computational
domain. We suppose that the wave of small amplitude can be obtained as a solution of
a linearized isentropic Euler equations which are simply the acoustic equations. Without
loss of generality we can assume that the mass density and speed of sound are equal to
unity. Then the acoustic equations with the added relaxation source terms read as follows:
∂p
∂t
+
∂u
∂x
+
∂v
∂y
= −p
τ
,
∂u
∂t
+
∂p
∂x
= −u
τ
, (4)
∂v
∂t
+
∂p
∂y
= −v
τ
.
3
Here p, u, v are non-dimensional pressure, and components of velocities in x and y direc-
tions, respectively. Inside the computational domain τ is supposed to be infinity (+∞),
hence source terms vanish, whereas in the PML τ is finite and small enough to provide a
high rate of decay of waves amplitude.
We shall study the solution to system (4) in the PML (x, y) ∈ [0, h] × (−∞,+∞)
supposing that the computational region is located from the left side of PML (x ≤ 0).
Following the analysis in [4] we suppose that a harmonic plane wave of the form
(p, u, v)T = (1, α, β)T exp(iω(t− αx− βy))
propagates from the main computational domain to the PML. Here α, β (α2 + β2 = 1)
represent the angle of wave incidence and ω is the normalized frequency of incoming wave.
We seek the solution inside the PML in the following form
(p, u, v)T = (P (x), U(x), V (x))T exp(iω(t− βy)).
Substituting this representation of the solution into (4) we obtain a system of ordinary
differential equations for P (x), U(x), V (x):
dP
dx
+
(
iω +
1
τ
)
U = 0,
dU
dx
+
(
iω +
1
τ
)
P − iωβV = 0, (5)
(
iω +
1
τ
)
V − iωβP = 0.
For the rest of the section we study the functions P (x), U(x) only. The solution of (5)
can be written as a combination of the two exponents:
P (x) = P1e
kx + P2e
−kx, U(x) = −P1 k
Ω
ekx + P2
k
Ω
e−kx, (6)
where
k =
√
Ω2 + β2ω2, Ω = iω +
1
τ
, (7)
and P1, P2 are constants which can be found from the boundary conditions for PML.
Our goal is to prove that the choice of the right boundary condition does not affect
the absorbtion properties of the PML provided the relaxation time τ is sufficiently small.
As the left boundary condition we take the value of the Riemann invariant which comes
to the PML from the computational domain. This means that this Riemann invariant is
continuous across the interface between the PML and computational domain. We remark
that this requirement is reasonable from the point of view of the theory of hyperbolic
equations. Thus we take
P (0) + U(0) = 1 + α.
4
Using the representation (6) of the solution we obtain the following relation between
P1 and P2:
P2 = −Ω− k
Ω + k
P1 + (1 + α)
Ω
Ω + k
.
Now expressions for P (x) and Q(x) can be transformed into
P (x) = P1
(Ω + k)ekx − (Ω− k)e−kx
Ω + k
+ (1 + α)
Ω
Ω + k
e−kx
U(x) = −P1k (Ω + k)e
kx + (Ω− k)e−kx
Ω(Ω + k)
+ (1 + α)
k
Ω + k
e−kx. (8)
Here only one constant P1 must be determined with the use of right boundary condition
for PML.
We are interested in the invariant P (x) − U(x) and its value at x = 0, because this
value gives an estimate of the amplitude of waves generated by the PML and propagating
into the main computational region. Such a wave propagates into main computational
region and must be damped by the PML. This invariant can be easily obtained from (8):
P (x)− U(x) = P1 (Ω + k)
2ekx − (Ω− k)2e−kx
Ω + k
+ (1 + α)
Ω− k
Ω + k
e−kx (9)
and its value at x = 0 is
P (0)− U(0) = P1 2Ωk
Ω + k
+ (1 + α)
Ω− k
Ω + k
. (10)
We shall study an asymptotic solution behavior assuming τ sufficiently small. Then
the following asymptotic formula can be used:
k =
1
τ
+ iω +
β2
3
ω2τ +O(τ 2) = Ω +
β2
3
ω2τ +O(τ 2). (11)
Now we derive the solution for the case of general right boundary condition for the
PML which we take in the form
aP (h) + bU(h) = 0,
where a, b are an arbitrary constants. Using (8) we obtain the value for P1:
P1 = − (1 + α)(aΩ + bk)e
−kh
(aΩ− bk)(Ω + k)ekh − (aΩ + bk)(aΩ− bk)e−kh . (12)
Hence the solution and the invariant P (x) − U(x) in particular can be obtained by sub-
stituting (12) into (8) and (9) accordingly.
Now the study of the asymptotic behavior should be based on the fact that
τ−γe
−h
τ → 0, if τ → 0,
5
where γ ≥ 0. Using this asymptotic behavior one can prove that
P1 v e
−2h
τ → 0, if τ → 0.
This leads us to the conclusion that the asymptotic behavior of the solution, invariant
P (x)− U(x) and its value at x = 0 are as follows:
P (x) v (1 + α) Ω
Ω + k
e−kx,
U(x) v (1 + α) k
Ω + k
e−kx,
P (x)− U(x) v (1 + α)Ω− k
Ω + k
e−kx,
P (0)− U(0) v (1 + α)Ω− k
Ω + k
.
In particular we have
P (0)− U(0) ' −(1 + α)β
2ω2
6
τ 2 +O(τ 3).
So the asymptotic value for the Riemann invariant which generates waves propagating
from the PML to the computational domain is of order τ 2 and does not depend on the
type of right boundary condition. Moreover the influence of the angle of incidence β and
frequency ω of incoming wave can be made negligible. This allows us to choose a very
simple numerical algorithm for the computations at the right PML boundary.
4 Numerical example
Here we show some numerical results as applied to the two-dimensional hyperbolic sys-
tems. The proposed idea is mostly easily implemented in the framework of one-step
Godunov-type methods. For background information see e.g. [11, 7]. The operator L in
(3) takes the following form:
Lij = −Fi+1/2,j − Fi−1/2,j
∆x
− Gi,j+1/2 −Gi,j−1/2
∆y
and the PML update formula is given by
Qn+1ij =
∆t
τ
Q∞
1 + ∆t
τ
+
Qnij +∆tLij(Qn)
1 + ∆t
τ
(13)
Alternatively, one can use advection schemes with Runge-Kutta time marching, e.g. [6].
In this case the relaxation step for the source term is executed after the all Runge-Kutta
6
stages are carried out. Effectively, we could regard this as a time-splitting procedure. The
whole time step dt is divided into two sub steps: i) carry out Runge-Kutta method ii)
apply the relaxation step as given by (13) but without the spatial operator:
Qn+1ij =
∆t
τ
Q∞ +QRKij
1 + ∆t
τ
(14)
where QRKij is the value of the vector of conservative variables, obtained after the Runge-
Kutta time stepping.
Below we use the one-step ADER3 scheme, see [12, 10] and references therein. This
scheme is uniformly third-order accurate in time and fifth order accurate in space.
We first tested our method on a two-dimensional explosion test problem for the non-
linear two-dimensional Euler equation from [11]. In the computations we take the width
of the PML layer to be equal to ten cells. The relaxation time is taken to be τ = 10−2
inside the layer when the equations are written in a conventional non-dimensional form.
This problem is an analog of the shock-tube problems in one space dimension. The ini-
tial condition consists of two regions of constant but different values of gas parameters
separated by a cylindrical surface. The solution involves a cylindrical shock wave leaving
the computational domain and is thus appropriate for assessing the robustness of the
method. Numerical results omitted here show that the proposed PML algorithm works
well without any stabilization or filtering used in [4, 5].
Secondly, we apply the method to a standard acoustic test problem [4, 5]. We solve
the two-dimensional linearized Euler equations of the form (1) with (again in the non-
dimensional form)
Q =

ρ
u
v
p
 , F (Q) =

u0 ρ0 0 0
0 u0 0 1/ρ0
0 0 u0 0
0 γp0 0 u0
Q, G(Q) =

v0 0 ρ0 0
0 v0 0 0
0 0 v0 1/ρ0
0 0 γp0 v0
Q
in a spatial domain of [−50, 50] × [−50, 50]. Here we take γ = 1.4, ρ0 = p0 = 1, v0 = 0,
u0 = 0.5
√
γ. The initial conditions include an acoustic pulse centered at (xa, ya) and a
vorticity and entropy pulses both centred at (xb, yb) and are given by
ρ = exp (−(ra/3)2 log 2) + 0.1 exp (−(rb/4)2 log 2),
p = γ exp (−(ra/3)2 log 2),
u =
√
γ 0.05(y − yb) exp (−(rb/4)2 log 2),
v = −√γ 0.05(x− xb) exp (−(rb/4)2 log 2)
where r2a,b = (x−xa,b)2+(y−ya,b)2, (xa, ya) = (−25, 0), (xb, yb) = (25, 0). Note that factors
proportional to γ appear due to the fact that our choice of non-dimensional variables is
different from that of [5, 4]. We use ∆x = ∆y = 1.
7
XY
-50 -25 0 25 50
-50
-25
0
25
50
X
Y
-50 -25 0 25 50
-50
-25
0
25
50
Figure 1: Density (left) and x component of velocity (right) for t = 25.
X
Y
-50 -25 0 25 50
-50
-25
0
25
50
X
Y
-50 -25 0 25 50
-50
-25
0
25
50
Figure 2: Density (left) and x component of velocity (right) for t = 50.
Figs. 1, 2 show the counter lines of density and x component of velocity for output
times: t = 25 and t = 50. Comparing our results with those reported in the literature [5, 4]
we conclude that their quality is comparable. No reflections from the boundary take place.
As time elapses, the pressure waves leave the computational domain. For example,
Fig. 3 illustrates pressure distribution along the x axis for t = 60. Here symbols correspond
to our numerical solution whereas the solid line represents the reference solution obtained
on a larger domain. As is seen, no spurious waves reflect from the PLM layer back to the
domain.
5 Conclusions
We have presented a new variant of the perfectly matched layer approach to construction
of non-reflecting boundary conditions. The method is exceedingly simple, robust, does
not involve altering of the governing equations and does not need any filters for stability.
Numerical results demonstrate that its performance is similar to the other PML schemes,
8
-100 -50 0 50 100
-0.02
0
0.02
0.04
0.06
0.08
ρ
x -100 -50 0 50 100
-0.01
-0.005
0
0.005
0.01
p
x
Figure 3: Density (left) and pressure (right) distribution along the x axis for t = 60.
presented in the literature.
Acknowledgements
The first author acknowledges the Italian Ministry of Education and Research (MIUR) for
the financial support provided during his stay as a visiting professor at the Department
of Civil and Environmental Engineering, University of Trento.
References
[1] J.-P. Berenger. A perfectly matched layer for the absorption of electro-magnetic
waves. J. Comput. Phys., 114:185–200, 1994.
[2] D.J. Bodony. Analysis of sponge zones for computational fluid mechanics. J. Comput.
Phys., 212(2):681–702, 2005.
[3] T. Colonius, S.K. Lele, and P. Moin. Boundary conditions for direct computation of
aerodynamic sound generation. AIAA Journal, 31:1574–1582, 1993.
[4] J.S. Hesthaven. On the analysis and construction of perfectly matched layers for the
linearized Euler equations. J. Comput. Phys., 142:129–147, 1998.
[5] F. Q. Hu. On absorbing boundary conditions for linearized Euler equations by a
perfectly matched layer. J. Comput. Phys., 129:201–219, 1996.
[6] G.S. Jiang and C.W. Shu. Efficient implementation of weighted ENO schemes. J.
Comput. Phys., 126:202–212, 1996.
9
[7] A.G. Kulikovskii, N. V. Pogorelov, and A. Yu. Semenov. Mathematical Aspects of
Numerical Solution of Hyperbolic Systems. Chapman and Hall, 2002. Monographs
and Surveys in Pure and Applied Mathematics, Vol. 118.
[8] V.S. Ryaben’kii, S.V. Tsynkov, and V.I. Turchaninov. Global discrete artificial
boundary conditions for time-dependent wave propagation. J. Comput. Phys.,
174:712–758, 2001.
[9] I.I. Sofronov. Artificial boundary conditions of absolute transparency for two- and
three-dimensional external time-dependent scattering problems. Eur. J. Appl. Math.,
9:561, 1998.
[10] V.A. Titarev and E.F. Toro. ADER schemes for three-dimensional nonlinear hyper-
bolic systems. J. Comput. Phys., 204(2):715–736, 2005.
[11] E.F. Toro. Riemann solvers and numerical methods for fluid dynamics. Springer-
Verlag, 1999. Second Edition.
[12] E.F. Toro, R.C. Millington, and L.A.M. Nejad. Towards very high order Godunov
schemes. In E. F. Toro, editor, Godunov Methods. Theory and Applications, pages
907–940. Kluwer/Plenum Academic Publishers, 2001.
10


Perfectly matched layers with high rate damping for hyperbolic systems

Unitn-eprints Research

PERFECTLY MATCHED LAYERS WITH HIGH RATEDAMPING FOR HYPERBOLIC SYSTEMSE. ROMENSKI 1, V.A. TITAREV and E. F. TOROLaboratory of Applied Mathematics, Faculty of Engineering,University of Trento, Italyromenski@ing.unitn.it; titarev@mail.ru; toro@ing.unitn.itAbstractWe propose a simple method for constructing non-reflecting boundary conditionsvia Perfectly Matched Layer approach. The basic idea of the method is to builda layer with high rate damping properties with are provided by adding the stiffrelaxation source terms to all equations of the system. No complicated modificationof the system to be solved is then required.1 IntroductionThe computational problem of wave propagation in infinite domains arises in acoustics,seismology and electromagnetic waves phenomena. The numerical study is usually refor-mulated in a bounded artificial computational domain in which infinity is modelled bysome properties of the numerical boundary or its neighborhood. The purpose of such areformulation of the problem is to avoid wave reflection from the numerical boundary.The problem of constructing numerical boundary conditions has been intensively stud-ied in recent decades. As a result, two main approaches have been developed. The firstapproach is the so called ”Absorbing Boundary Conditions” (ABC) approach and con-sists of formulating direct boundary conditions for the computational region, which wouldeliminate wave reflection from the boundary. We refer the reader to some recent papers[8, 9] and bibliography therein.Another approach is the ”Perfectly Matched Layer” (PML) approach in which thecomputational domain is surrounded by an additional boundary layer. The main ideais to let the waves propagate out of the computational region. Therefore, it is crucialthat the boundary layer does not generate waves propagating back to the computationalregion. Various versions of the PML approach have been proposed in the past, which arebased on introducing new artificial variables and differential equations for these variables.See [5, 4, 2] and references therein.In this paper we propose a simple, reliable and efficient PML method suitable forsolving wave propagation problems described by hyperbolic systems, either linear or non-linear. The basic idea of the method is to formulate a perfectly matched layer as a layer1On leave from Sobolev Institute of Mathematics, Russian Academy of Sciences, Novosibirsk, Russia1with high rate damping properties. These properties are provided by adding the stiffrelaxation source terms to all equations of the system. An asymptotic analysis for smallvalue of relaxation time τ allow us to make a conclusion that the amplitude of the wavereflected from the PML is of the second order on τ . Therefore the influence of the angleof incidence and wavelength of the wave incoming to the PML from the computationaldomain on the reflected wave can be made negligible. The main advantages of the methodproposed here are its simplicity and generality. No complicated modification of the systemto be solved is required. The structure of the source terms allows one to use practicallyany advection scheme without significant changes. No stability problem exists either.A similar idea to use source terms with soft damping rate of acoustic waves in the cer-tain and rather large buffer layer (comparable with computational domain) was proposedin [3]. In the quoted paper the mesh stretching and filtering together with non-reflectingboundary conditions have been used in order to provide outflow disturbance without re-flection for aerodynamic sound generation problems. The absorbing layers with relaxationdamping have been proposed also for electromagnetic equations [1] and for linearized Eulerequations [5]. In these papers absorbing technique requires additional reformulation thegoverning equations (splitting in the coordinate direction).The rest of the paper organized as follows. In Section 2 we formulate the idea of thenew perfectly matched layer for the general case of nonlinear systems of conservation lawsand the method of its numerical solving. In Section 3 we analyze the idea as applied tolinear acoustics. In Section 4 we discuss the application of the method to two-dimensionalequations and present some numerical examples. Conclusions are drawn in Section 5.2 Equations for perfectly matched layerWe study the wave propagation phenomena described by a hyperbolic system of conserva-tion laws in an infinite spacial domain. Generally speaking, the system can be nonlinear.In the three dimensional Cartesian coordinate system xi it can be written as∂∂tQ+∂∂xF (Q) +∂∂yG(Q) +∂∂zH(Q) = 0, (1)where Q is the vector of conserved variables, F (Q), G(Q), H(Q), are the flux vectors inthe coordinate directions.To solve the problem numerically we apply the PML strategy, which consists of defin-ing the finite computational domain D in which we intend to obtain a solution, and thenconstructing a surrounding absorbing boundary layer of a prescribed width. This bound-ary layer must not affect the basic computational domain, meaning that there are nowaves coming back to the basic domain from the boundary layer.Suppose that the values Q∞ of the conserved variables at infinity are known: Q→ Q∞if xi → ∞. We formulate the equations which describe the wave propagation in the2perfectly matched layer by adding a relaxation source term to (1):∂∂tQ+∂∂xF (Q) +∂∂yG(Q) +∂∂zH(Q) = −1τ(Q−Q∞), (2)where τ is the relaxation time which can be a function of variables Q. We take τ = ∞inside our computational domain whereas in the PML layer τ is taken to be small.The idea of adding such relaxation term is to provide rapid (exponential) dampingof all variables and waves inside the layer. It is clear intuitively that the decrease of thewave amplitude will be faster if the relaxation time tends to zero. Note that the structureand the type of equations does not change making it possible to use the same numericalmethod as for the original system (1).We now proceed to describe the numerical procedure to be used in the new PML.Suppose that we have a one-step numerical scheme for solving (1), which gives us a resultin the formQn+1 −Qn∆t= L(Qn).Here ∆t is the time step. The simplest numerical procedure is the one step implicitapproximation of the source term in the system (2). That is we use the following formula:Qn+1 −Qn∆t= L(Qn)− 1τ(Qn+1 −Q∞). (3)It gives us the result in the form(1 +∆tτ)Qn+1 =∆tτQ∞ +Qn +∆tL(Qn)¿From the latter formula we see that if the relaxation time τ is very small (in particularτ much less than ∆t) then the value of Qn+1 is very close to its value at infinity Q∞.3 Acoustic wave propagation in the PMLIn this section we study the method as applied to acoustic waves leaving the computationaldomain. We suppose that the wave of small amplitude can be obtained as a solution ofa linearized isentropic Euler equations which are simply the acoustic equations. Withoutloss of generality we can assume that the mass density and speed of sound are equal tounity. Then the acoustic equations with the added relaxation source terms read as follows:∂p∂t+∂u∂x+∂v∂y= −pτ,∂u∂t+∂p∂x= −uτ, (4)∂v∂t+∂p∂y= −vτ.3Here p, u, v are non-dimensional pressure, and components of velocities in x and y direc-tions, respectively. Inside the computational domain τ is supposed to be infinity (+∞),hence source terms vanish, whereas in the PML τ is finite and small enough to provide ahigh rate of decay of waves amplitude.We shall study the solution to system (4) in the PML (x, y) ∈ [0, h] × (−∞,+∞)supposing that the computational region is located from the left side of PML (x ≤ 0).Following the analysis in [4] we suppose that a harmonic plane wave of the form(p, u, v)T = (1, α, β)T exp(iω(t− αx− βy))propagates from the main computational domain to the PML. Here α, β (α2 + β2 = 1)represent the angle of wave incidence and ω is the normalized frequency of incoming wave.We seek the solution inside the PML in the following form(p, u, v)T = (P (x), U(x), V (x))T exp(iω(t− βy)).Substituting this representation of the solution into (4) we obtain a system of ordinarydifferential equations for P (x), U(x), V (x):dPdx+(iω +1τ)U = 0,dUdx+(iω +1τ)P − iωβV = 0, (5)(iω +1τ)V − iωβP = 0.For the rest of the section we study the functions P (x), U(x) only. The solution of (5)can be written as a combination of the two exponents:P (x) = P1ekx + P2e−kx, U(x) = −P1 kΩekx + P2kΩe−kx, (6)wherek =√Ω2 + β2ω2, Ω = iω +1τ, (7)and P1, P2 are constants which can be found from the boundary conditions for PML.Our goal is to prove that the choice of the right boundary condition does not affectthe absorbtion properties of the PML provided the relaxation time τ is sufficiently small.As the left boundary condition we take the value of the Riemann invariant which comesto the PML from the computational domain. This means that this Riemann invariant iscontinuous across the interface between the PML and computational domain. We remarkthat this requirement is reasonable from the point of view of the theory of hyperbolicequations. Thus we takeP (0) + U(0) = 1 + α.4Using the representation (6) of the solution we obtain the following relation betweenP1 and P2:P2 = −Ω− kΩ + kP1 + (1 + α)ΩΩ + k.Now expressions for P (x) and Q(x) can be transformed intoP (x) = P1(Ω + k)ekx − (Ω− k)e−kxΩ + k+ (1 + α)ΩΩ + ke−kxU(x) = −P1k (Ω + k)ekx + (Ω− k)e−kxΩ(Ω + k)+ (1 + α)kΩ + ke−kx. (8)Here only one constant P1 must be determined with the use of right boundary conditionfor PML.We are interested in the invariant P (x) − U(x) and its value at x = 0, because thisvalue gives an estimate of the amplitude of waves generated by the PML and propagatinginto the main computational region. Such a wave propagates into main computationalregion and must be damped by the PML. This invariant can be easily obtained from (8):P (x)− U(x) = P1 (Ω + k)2ekx − (Ω− k)2e−kxΩ + k+ (1 + α)Ω− kΩ + ke−kx (9)and its value at x = 0 isP (0)− U(0) = P1 2ΩkΩ + k+ (1 + α)Ω− kΩ + k. (10)We shall study an asymptotic solution behavior assuming τ sufficiently small. Thenthe following asymptotic formula can be used:k =1τ+ iω +β23ω2τ +O(τ 2) = Ω +β23ω2τ +O(τ 2). (11)Now we derive the solution for the case of general right boundary condition for thePML which we take in the formaP (h) + bU(h) = 0,where a, b are an arbitrary constants. Using (8) we obtain the value for P1:P1 = − (1 + α)(aΩ + bk)e−kh(aΩ− bk)(Ω + k)ekh − (aΩ + bk)(aΩ− bk)e−kh . (12)Hence the solution and the invariant P (x) − U(x) in particular can be obtained by sub-stituting (12) into (8) and (9) accordingly.Now the study of the asymptotic behavior should be based on the fact thatτ−γe−hτ → 0, if τ → 0,5where γ ≥ 0. Using this asymptotic behavior one can prove thatP1 v e−2hτ → 0, if τ → 0.This leads us to the conclusion that the asymptotic behavior of the solution, invariantP (x)− U(x) and its value at x = 0 are as follows:P (x) v (1 + α) ΩΩ + ke−kx,U(x) v (1 + α) kΩ + ke−kx,P (x)− U(x) v (1 + α)Ω− kΩ + ke−kx,P (0)− U(0) v (1 + α)Ω− kΩ + k.In particular we haveP (0)− U(0) ' −(1 + α)β2ω26τ 2 +O(τ 3).So the asymptotic value for the Riemann invariant which generates waves propagatingfrom the PML to the computational domain is of order τ 2 and does not depend on thetype of right boundary condition. Moreover the influence of the angle of incidence β andfrequency ω of incoming wave can be made negligible. This allows us to choose a verysimple numerical algorithm for the computations at the right PML boundary.4 Numerical exampleHere we show some numerical results as applied to the two-dimensional hyperbolic sys-tems. The proposed idea is mostly easily implemented in the framework of one-stepGodunov-type methods. For background information see e.g. [11, 7]. The operator L in(3) takes the following form:Lij = −Fi+1/2,j − Fi−1/2,j∆x− Gi,j+1/2 −Gi,j−1/2∆yand the PML update formula is given byQn+1ij =∆tτQ∞1 + ∆tτ+Qnij +∆tLij(Qn)1 + ∆tτ(13)Alternatively, one can use advection schemes with Runge-Kutta time marching, e.g. [6].In this case the relaxation step for the source term is executed after the all Runge-Kutta6stages are carried out. Effectively, we could regard this as a time-splitting procedure. Thewhole time step dt is divided into two sub steps: i) carry out Runge-Kutta method ii)apply the relaxation step as given by (13) but without the spatial operator:Qn+1ij =∆tτQ∞ +QRKij1 + ∆tτ(14)where QRKij is the value of the vector of conservative variables, obtained after the Runge-Kutta time stepping.Below we use the one-step ADER3 scheme, see [12, 10] and references therein. Thisscheme is uniformly third-order accurate in time and fifth order accurate in space.We first tested our method on a two-dimensional explosion test problem for the non-linear two-dimensional Euler equation from [11]. In the computations we take the widthof the PML layer to be equal to ten cells. The relaxation time is taken to be τ = 10−2inside the layer when the equations are written in a conventional non-dimensional form.This problem is an analog of the shock-tube problems in one space dimension. The ini-tial condition consists of two regions of constant but different values of gas parametersseparated by a cylindrical surface. The solution involves a cylindrical shock wave leavingthe computational domain and is thus appropriate for assessing the robustness of themethod. Numerical results omitted here show that the proposed PML algorithm workswell without any stabilization or filtering used in [4, 5].Secondly, we apply the method to a standard acoustic test problem [4, 5]. We solvethe two-dimensional linearized Euler equations of the form (1) with (again in the non-dimensional form)Q =ρuvp , F (Q) =u0 ρ0 0 00 u0 0 1/ρ00 0 u0 00 γp0 0 u0Q, G(Q) =v0 0 ρ0 00 v0 0 00 0 v0 1/ρ00 0 γp0 v0Qin a spatial domain of [−50, 50] × [−50, 50]. Here we take γ = 1.4, ρ0 = p0 = 1, v0 = 0,u0 = 0.5√γ. The initial conditions include an acoustic pulse centered at (xa, ya) and avorticity and entropy pulses both centred at (xb, yb) and are given byρ = exp (−(ra/3)2 log 2) + 0.1 exp (−(rb/4)2 log 2),p = γ exp (−(ra/3)2 log 2),u =√γ 0.05(y − yb) exp (−(rb/4)2 log 2),v = −√γ 0.05(x− xb) exp (−(rb/4)2 log 2)where r2a,b = (x−xa,b)2+(y−ya,b)2, (xa, ya) = (−25, 0), (xb, yb) = (25, 0). Note that factorsproportional to γ appear due to the fact that our choice of non-dimensional variables isdifferent from that of [5, 4]. We use ∆x = ∆y = 1.7XY-50 -25 0 25 50-50-2502550XY-50 -25 0 25 50-50-2502550Figure 1: Density (left) and x component of velocity (right) for t = 25.XY-50 -25 0 25 50-50-2502550XY-50 -25 0 25 50-50-2502550Figure 2: Density (left) and x component of velocity (right) for t = 50.Figs. 1, 2 show the counter lines of density and x component of velocity for outputtimes: t = 25 and t = 50. Comparing our results with those reported in the literature [5, 4]we conclude that their quality is comparable. No reflections from the boundary take place.As time elapses, the pressure waves leave the computational domain. For example,Fig. 3 illustrates pressure distribution along the x axis for t = 60. Here symbols correspondto our numerical solution whereas the solid line represents the reference solution obtainedon a larger domain. As is seen, no spurious waves reflect from the PLM layer back to thedomain.5 ConclusionsWe have presented a new variant of the perfectly matched layer approach to constructionof non-reflecting boundary conditions. The method is exceedingly simple, robust, doesnot involve altering of the governing equations and does not need any filters for stability.Numerical results demonstrate that its performance is similar to the other PML schemes,8-100 -50 0 50 100-0.0200.020.040.060.08ρx -100 -50 0 50 100-0.01-0.00500.0050.01pxFigure 3: Density (left) and pressure (right) distribution along the x axis for t = 60.presented in the literature.AcknowledgementsThe first author acknowledges the Italian Ministry of Education and Research (MIUR) forthe financial support provided during his stay as a visiting professor at the Departmentof Civil and Environmental Engineering, University of Trento.References[1] J.-P. Berenger. A perfectly matched layer for the absorption of electro-magneticwaves. J. Comput. Phys., 114:185–200, 1994.[2] D.J. Bodony. Analysis of sponge zones for computational fluid mechanics. J. Comput.Phys., 212(2):681–702, 2005.[3] T. Colonius, S.K. Lele, and P. Moin. Boundary conditions for direct computation ofaerodynamic sound generation. AIAA Journal, 31:1574–1582, 1993.[4] J.S. Hesthaven. On the analysis and construction of perfectly matched layers for thelinearized Euler equations. J. Comput. Phys., 142:129–147, 1998.[5] F. Q. Hu. On absorbing boundary conditions for linearized Euler equations by aperfectly matched layer. J. Comput. Phys., 129:201–219, 1996.[6] G.S. Jiang and C.W. Shu. Efficient implementation of weighted ENO schemes. J.Comput. Phys., 126:202–212, 1996.9[7] A.G. Kulikovskii, N. V. Pogorelov, and A. Yu. Semenov. Mathematical Aspects ofNumerical Solution of Hyperbolic Systems. Chapman and Hall, 2002. Monographsand Surveys in Pure and Applied Mathematics, Vol. 118.[8] V.S. Ryaben’kii, S.V. Tsynkov, and V.I. Turchaninov. Global discrete artificialboundary conditions for time-dependent wave propagation. J. Comput. Phys.,174:712–758, 2001.[9] I.I. Sofronov. Artificial boundary conditions of absolute transparency for two- andthree-dimensional external time-dependent scattering problems. Eur. J. Appl. Math.,9:561, 1998.[10] V.A. Titarev and E.F. Toro. ADER schemes for three-dimensional nonlinear hyper-bolic systems. J. Comput. Phys., 204(2):715–736, 2005.[11] E.F. Toro. Riemann solvers and numerical methods for fluid dynamics. Springer-Verlag, 1999. Second Edition.[12] E.F. Toro, R.C. Millington, and L.A.M. Nejad. Towards very high order Godunovschemes. In E. F. Toro, editor, Godunov Methods. Theory and Applications, pages907–940. Kluwer/Plenum Academic Publishers, 2001.10

http://eprints.biblio.unitn.it/974/1/UTM695.pdf

Perfectly matched layers with high rate damping for hyperbolic systems

Abstract

Similar works

Full text

Available Versions

University of Trento : UNITN-eprints

Unitn-eprints Research