Ultrasonics xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Ultrasonics journal homepage: www.elsevier.com/locate/ultras

A novel unsplit perfectly matched layer for the second-order acoustic wave equation Youneng Ma a,b, Jinhua Yu a,b,⇑, Yuanyuan Wang b a b

National Key Laboratory of ASIC and System of Fudan University, Shanghai 200030, China Department of Electronic Engineering, Fudan University, Shanghai 200433, China

a r t i c l e

i n f o

Article history: Received 14 January 2013 Received in revised form 18 March 2014 Accepted 31 March 2014 Available online xxxx Keywords: Auxiliary-differential equation Acoustic Perfectly matched layer Second-order equation

a b s t r a c t When solving acoustic field equations by using numerical approximation technique, absorbing boundary conditions (ABCs) are widely used to truncate the simulation to a finite space. The perfectly matched layer (PML) technique has exhibited excellent absorbing efficiency as an ABC for the acoustic wave equation formulated as a first-order system. However, as the PML was originally designed for the first-order equation system, it cannot be applied to the second-order equation system directly. In this article, we aim to extend the unsplit PML to the second-order equation system. We developed an efficient unsplit implementation of PML for the second-order acoustic wave equation based on an auxiliary-differential-equation (ADE) scheme. The proposed method can benefit to the use of PML in simulations based on second-order equations. Compared with the existing PMLs, it has simpler implementation and requires less extra storage. Numerical results from finite-difference time-domain models are provided to illustrate the validity of the approach. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Numerical simulations of wave propagation based on partial differential equations require the truncation of the computational domain due to the limited memory and computation time of computers. The truncation can lead to artificial reflections which can ruin the simulating results. Absorbing boundary conditions (ABCs) are routinely used at the truncated boundary to attenuate the spurious reflections. Many ABC techniques have been developed to complete this task during past years. Among them, the PML technique, first introduced by Berenger [1] in 1994, has been proven to be a very efficient scheme for the truncation of unbounded domain. Berenger’s original formulation is called the split-field PML, because he artificially split the wave solutions into the sum of two new artificial field components. Though first designed for the Maxwell equations in the context of electromagnetics, the PML has now been widely used in simulations of electromagnetic, elastic, as well as acoustic wave propagations [2,3] for its excellent absorbing efficiency. Following Berenger’ classical split-field scheme, a number of different implementations have been introduced [4–8], leading to the unsplit implementations of the PML. Though these implementations often ⇑ Corresponding author at: Department of Electronic Engineering, Fudan University, Shanghai 200433, China. Tel./fax: +86 21 65643526. E-mail address: [email protected] (J. Yu).

lead to equivalent reflection properties, they offer different mathematical representations which may reduce implementation costs and complexity of the PML [4,7]. The auxiliary-differential-equation (ADE) [7,8] approach is one of important techniques that have facilitated implementation of the PML. It provides a more flexible representation of the PML that can be extended to higher-order methods [7]. As the classical PML was primarily designed for first-order equation system, it cannot be applied to the second-order system directly. Over past years, in the field of acoustic simulation, the PML has been mostly applied to acoustic wave equations formulated as a first-order system, rarely to equations written as a second-order system in pressure. In certain cases, a second-order system can be rewritten as a group of first-order ones and the PML can then be directly applied. However, this can lead to complex implementation, since the second-order equation system is simpler and more convenient to be used in many simulations. What is more, some well-established equations, which are often used in simulations [9–11], are not appropriate to be rewritten. Therefore, it is meaningful to extend the PML to the second-order acoustic equation. Several attempts have been reported to extend the PML for second-order equations. In 2003, Komatitsch and Tromp [12] developed a split PML for the second-order seismic wave equation. This is first effort made to apply the PML to second-order equations. This method was applied to the acoustic equation in 2010 [13]. Though Komatitsch’ scheme can attenuate

http://dx.doi.org/10.1016/j.ultras.2014.03.016 0041-624X/Ó 2014 Elsevier B.V. All rights reserved.

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

2

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

propagating waves effectively, the implementation costs of his method is very expensive since one third-order and two secondorder ADEs have to be solved. In 2009, Pinton et al. [11] developed an unsplit PML for the second-order acoustic equation and used it in ultrasonic imaging simulation. Though memory-efficient and seemingly simple, this mixed formulation of PML is not convenient to implement in practice for a deconvolving process involved and it is also difficult to extend to higher-order methods. In [14,15], the convolutional PML was formulated for the second-order elastic wave equation. As these two methods were designed for finite element methods (FEMs), it is not convenient to apply them in solving equations based on finite-difference time-domain (FDTD) methods or the pseudospectral time-domain (PSTD) [16] method. In this paper, a memory-efficient unsplit PML with a very easy implementation is developed for the second-order acoustic equation. The proposed method is formulated via the complex coordinate stretching scheme and introducing ADEs. It is a method suitable for numerical methods such as the FDTD, PSTD and FEM. The remainder of the paper is organized as follows. In Section 2 the background of the proposed method is introduced, and the formulation and implementation of the proposed method is presented. The memory requirement comparison between different schemes is also carried out. In Section 3, the proposed scheme is verified and compared with classical split PML for the secondorder equation in numerical experiments, followed by conclusions in Section 4.

2.2. The proposed method The traditional split-field formulation of PML can actually be viewed as an ADE-based scheme, since the splitting of pressure u is also a way to introduce several new auxiliary variables to simplify the formulation. Here we present a novel ADE-based unsplit ^ i as the implementation for (4). In the following part, we denote u FT of ui (i = 1, 2, 3). For simplicity of expression, we introduce a temporary variable p and let

p ¼ jx þ r x :

ð5Þ

Substituting p for jx on the right-side of (4), we have

^ ðp  rx Þ2 @ 2 u ^ @2u ^ @2u ^ 1 ðp  rx Þ2 r0x @ u 2 ^¼ ðjxÞ u þ þ 2þ 2: 2 3 2 2 @x @x @y @z c p p

Enforcing a partial expansion scheme on the right-side of the above equation and arrange its terms according to the order of p, Eq. (6) becomes

^ @2u ^ @2u ^ 1 @2u 2 ^¼ ðjxÞ u þ 2þ 2 2 2 @x @y @z c þ

r2x @ 2 u^ 2rx r0x p2 @x2

þ

p2

!

r0x @ u^ 2rx @ 2 u^

þ p @x p ! ^ @u r0 r2  x3 x p @x

^ @2u ^ @2u ^ 1 @2u 2 ^¼ ^1 ðjxÞ u þ 2 þ 2 u 2 2 c @x @y @z

2.1. Background

where u1 satisfies

The differential form of the linear acoustic wave equation in the frequency domain can be written as

^1 ¼ u

ð1Þ

where x denotes the angular frequency, c is the speed of the sound, ^ represents the Fourier transform (FT) of the acoustic pressure u, u and j is the imaginary unit. In the interest of notational simplicity, we only describe the solution in the x-axis using the stretched coordinate approach in [5,17]. In the PML region, Eq. (1) is formulated in complex stretched coordinates as

  ^ ^ @2u ^ 1 1 @ 1 @u @2u 2 ^¼ þ 2þ 2 ðj x Þ u 2 @y @z c sx @x sx @x

ð2Þ

^ @u @x

ð7Þ

p @x

þ2

rx @ 2 u^ p @x2

 rx

! ^ rx @ 2 u ^ 2r0x @ u r2 r0 @ u^ þ x3 x þ : 2 2 2 p @x p @x p @x

ð8Þ

ð9Þ

Multiplying both sides of (9) with p, we have

^ 1 ¼ r0x pu

^ ^ @u @2u ^2 þ 2r x 2  r x u @x @x

ð10Þ

where

! ^ rx @ 2 u ^ 2r0x @ u rx r0 @ u^  2x þ : p @x p @x2 p @x

^2 ¼ u

ð11Þ

Multiplying both sides of the above equation with p, Eq. (11) becomes

^ ^ @u @2u ^3 þ rx 2  rx u @x @x

ð12Þ

where

sx ¼ 1 þ rx =jx:

ð3Þ ^3 ¼ u

Here rx is the damping profile across the PML region. As demonstrated in [12,13], Eqs. (2) and (3) finally lead to 2

r0x @ u^

^ 2 ¼ 2r0x pu

where x-axis is stretched by

@x2

Introducing an auxiliary u1 and rearranging the above equation, we rewrite (7) as

2. Method

^ @2u ^ @2u ^ 1 @2u 2 ^¼ ðjxÞ u þ 2þ 2 2 2 c @x @y @z

ð6Þ

2

^ ^ @2u ^ @2u ^ 1 ðjxÞ r0x @ u ðjxÞ @2u 2 ^¼ ðjxÞ u þ þ 2þ 2 3 @x 2 @x2 c2 @y @z ðjx þ rx Þ ðjx þ rx Þ

ð4Þ

where rx0 is the partial derivative of rx with respect to x. It is difficult to derive the time-domain expression of (4) directly for its complex formulation. The general way is to introduce auxiliary variables. A split-field scheme is proposed in [12,13], where u is decomposed into 3 parts; each part is tackled separately, and the results are then combined. Though this split-field approach exhibits excellent absorbing efficiency, its numerical costs, especially the extra needed storage, are remarkably high. What is more, the splitting of pressure will often lead to cumbersome implementation when using FEM as discretization scheme [14].

r0x @ u^ p @x

ð13Þ

:

According to (13), we can easily get

^ 3 ¼ r0x pu

^ @u : @x

ð14Þ

Inserting (5) into (8), (10), (12), and (14), we obtain the following frequency domain equations in Cartesian coordinates

^ 3 ¼ r0x ðjx þ rx Þu

^ @u @x ^ ^ @u @2u ^3 þ rx 2  rx u @x @x

ð15bÞ

^ ^ @u @2u ^2 þ 2rx 2  rx u @x @x

ð15cÞ

^ 2 ¼ 2r0x ðjx þ rx Þu ^ 1 ¼ r0x ðjx þ rx Þu

ð15aÞ

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

3

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

( ) ^ @2u ^ @2u ^ @2u 2 2 ^¼c ^1 : ðjxÞ u þ þ u @x2 @y2 @z2

ð15dÞ

Transforming the above equations into the time-domain, we can finally obtain the time-domain expression as

ð@ t þ rx Þu3 ¼ r0x

@u @x

ð@ t þ rx Þu2 ¼ 2r0x

Secondly, enforcing a partial fraction expansion scheme on (19) and then adjusting the equation gives

^x ^z ^x @2u @2u @2u þ ðk þ lÞ þ ðk þ 2lÞ 2 2 @z @x@z !@x ! r0x @ u^x 2rx @ 2 u^x rx @ 2 u^z þ ðk þ  ðk þ 2lÞ þ l Þ p @x p @x2 p @x@z !! 2 0 rx @ u^x 2rx @ u^x  2 þ ðk þ 2lÞrx 2 p @x2 p @x

qðixÞ2 u^ x ¼ l

ð16aÞ

@u @2u þ rx 2  rx u3 @x @x

ð16bÞ

 r0x ðk þ 2lÞ

2

@u @ u þ 2rx 2  rx u2 @x @x ( ) @2u @2u @2u @ 2t u ¼ c2 þ þ  u 1 @x2 @y2 @z2 ð@ t þ rx Þu1 ¼ r0x

ð16cÞ

ð16dÞ

where @ t represents the time derivative. In this formulation, the pressure u is not split which can lead to easier implementation. Considering that the original split PML [12,13] needs the solution of three second-order and one first-order auxiliary equations, it is worthwhile to mention that this formulation contains only three first-order auxiliary Eqs. (16a)–(16c), which can reduce the extra need of memory and facilitate its implementation to a large extend. We also mention that the differential equations in (16) must be tackled in sequence in each time step during the implementation. Since u3 depends only on u, it should be updated firstly. Then u2, which depends on u and u1, can be calculated, followed by u1 and u. The PML in the y-, z-axis could be formulated in a same manner as for the x-axis. As shown above, the proposed PML formulation can be roughly divided into three steps: 1. Derive the stretched function of the original function using the complex-coordinate-stretching technique (Eqs. (1)–(4)); 2. Enforce a partial fraction expansion on the stretched function and arrange its terms (Eqs. (5) and (6)); 3. Introduce ADEs step by step (Eqs. (8)–(14)). In this process, the introduced ADEs are of first order. Generally, only three ADEs are needed in each axis. Though we only describe how to derive the PML formulation of acoustic equation above. It is also possible to extend the method to the other second-order equations using the above procedures. Take the second-order elastic equation in isotropic media as an example. The original elastic equation is given by

q

^x @ 2 ux @2u @ 2 ux @ 2 uz ¼ ðk þ 2lÞ 2 þ l 2 þ ðk þ lÞ 2 @x @z @x@z @t

ð17aÞ

q

@ 2 uz @ 2 uz @ 2 uz @ 2 ux ¼ l 2 þ ðk þ 2lÞ 2 þ ðk þ lÞ 2 @x @z @x@z @t

ð17bÞ

  ^x ^x ^z 1 @ 1 @u @2u 1 @2u þ l 2 þ ðk þ lÞ qðixÞ u^x ¼ ðk þ 2lÞ @z sx @x sx @x sx @x@z 2

ð18Þ

Referring to [12], Eq. (18) can lead to: 2 ^x r0x @ u^x ðjxÞ @2u þ 3 @x 2 @x2 ðjx þ rx Þ ðjx þ rx Þ ^x ^z @2u jx @2u : þ l 2 þ ðk þ lÞ @z jx þ rx @x@z

qðixÞ2 u^x ¼ ðk þ 2lÞ 

ðjxÞ

2

!

p3 @x

ð20Þ

;

where p = jx + rx. Thirdly, introducing auxiliary variables step by step (similar to the processes from (7) to (15)). The following equations can be finally obtained:

^ 3x ¼ r0x pu

^x @u @x

ð21aÞ

^ 2x ¼ rx pu

^x ^x @2u @u ^ 3x  2r0x  rx u @x2 @x

ð21bÞ

^z @2u ^ 1x ¼ ðk þ lÞrx pu þ ðk þ 2lÞ @x@z

qðixÞ2 u^ x ¼ l

^ @x

0 @ ux þ2 x

r

^ @2u rx 2x  rx u^2x @x

^x ^z ^x @2u @2u @2u ^ 1x : þ ðk þ lÞ þ ðk þ 2lÞ 2  u @z2 @x@z @x

ð21cÞ

ð21dÞ

2.3. Implementation The proposed formulations can be discretized using schemes such as PSTD, FEM and FDTD conveniently, which we can see from its expression. Here we only present an FDTD update equations of (16) (grid positions are suppressed for brevity)

@ n u @x

u3nþ1=2 ¼ C 1 un1=2 þ C2 3 nþ1=2

u2

nþ1=2 u1

n1=2

¼ C 1 u2

¼

n1=2 C 1 u1

þ 2C 2

ð22aÞ

@ n @2 n nþ1=2 u þ C3 u  u3 @x @x

!

@ @2 nþ1=2 þ C 2 un þ C 3 2 un  u2 @x @x

ð22bÞ ! ð22cÞ

" nþ1

n

n1

¼ 2u  u

2 2

þ Dt c

# @2 n @2 n @2 n nþ1=2 ; u þ 2 u þ 2 u  u1 @x @y @z

ð22dÞ

where

C1 ¼

1  rx Dt=2 r0x Dt r x Dt ; C2 ¼ ; C3 ¼ : 1 þ rx Dt=2 1 þ rx Dt=2 1 þ rx Dt=2

The spatial derivatives in the x-axis can be discretized using a fourth-order-accuracy finite difference scheme: n n n n @un uiþ2;j;k þ 8uiþ1;j;k  8ui1;j;k þ ui2;j;k ¼ @x 12Dx n

n

n

n

n

ð23aÞ n

uiþ2;j;k þ 16uiþ1;j;k  30ui;j;k þ 16ui1;j;k  ui2;j;k @2u  : @x2 12Dx2 ð19Þ

!

Eq. (21) can be then be transformed directly into time domain and the PML equations are obtained.

u

where ux, uz are displacements in the horizontal and vertical direction respectively, q is the density of the media and k, l are elastic coefficients. We only show how to derive the corresponding PML equation of (17a) in the x-axis using the proposed procedures. PML equations of (17b) and in other axes can be obtained similarly. Firstly, stretching x-axis in (17a) yields

r2x @ u^x

ð23bÞ

Other spatial derivatives can be approximated similarly. Note that the constants C1, C2 and C3 can be calculated and stored first, and then be used at each step. It should be pointed

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

4

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

out that, in each time step, with several temporary memory variables introduced, partial derivatives in (22) need to be calculated only once at each step and what is more, the partial derivatives can be approximated using the arbitrary finite-order difference quite conveniently. 2.4. Memory requirement comparison In order to perform the comparison between the proposed method and other existing ones in terms of memory and simplicity, we assume a vacuum explicit FDTD computational domain of N  N cells surrounded by an M-cell thick PML on each side of the simulation domain. It is easy to find that N > 2M. In the PML regions, the proposed method requires 12MN variables (auxiliary variables only) in total; the approach in [11] requires 16MN variables; the split-field PML in [12,13] requires 40MN  16M2 variables. Thus the proposed method requires 28MN  16M2 less than the split-field method and 4MN less than the method in [11]. Furthermore, as we can see from their numerical expression, the proposed scheme and the method in [12,13] can be discretized by the arbitrary finite-order difference and the PSTD algorithm conveniently. However, the method in [11] is not convenient to discretize, for it contains a tricky hybrid partial derivative term related to both the time and space. This tricky term cannot be discretized using the PSTD [16] which calculates the partial derivative globally with high efficiency using the Fast Fourier transformation. It is also difficult to tackle it when using other numerical methods. Compared with the method in [12,13], the proposed method has the advantage of simpler implementation. In each axis, it only needs to calculate three first-order and one second-order partial differential equations, whereas the method in [12] needs to calculate one first-order and three second-order partial derivative equations. In a word, the proposed method is not only memoryefficient and suitable for different numerical methods, but also easy to implement.

Fig. 1. The homogeneous computational model. S is the point source. A and B are receivers, placed two cells away from the PML region. The total domain has an area of 2.31  2.31 cm and is discretized into 150  150 grids with an additional of Mcell thick PML.

The PML are tested with the thickness of 30, 50 and 70 cells with R0 = 5  103, 1  105, 1  105 respectively. Fig. 2 presents snapshots of the acoustic pressure propagation in the computational domain with various time steps where M is 70. This figure clearly indicates that propagating waves can be well absorbed by the proposed PML. Fig. 2(a) illustrates waves emit from the source. Fig. 2(b) and (c) presents that few spurious reflections generated when the wave comes across the interface between PML domain and the interior simulation domain. Fig. 2(d) indicates that the waves are almost all absorbed by the PML in the latter time. In order to evaluate the efficiency of the PML quantitatively, the total energy E is computed at each time step



N X N X

u2 ði; jÞ

ð26Þ

i¼1 j¼1

3. Numerical experiments 3.1. Validation of absorbing efficiency in homogeneous media To validate the efficiency of the PML system, a 2-D numerical simulation is carried out, where Eq. (16) are discretized as shown in Section 2.3. The simulation field is demonstrated in Fig. 1. It has a dimension of 2.31  2.31 cm, and the field is divided into 150  150 grid elements with uniform cell spacing. An addition of M-cell thick PML is placed on its all four sides with the same cell spacing. S is the point source and A, B are receivers. With the very left bottom grid set as (1, 1) and the cell spacing as the unit length, their grid locations are indicated in Fig. 1. In the simulation, a time step of 50 ns is used and the speed of the sound c is set as 1540 m/ s. In the PML, we use the empirical damping profile

rx ¼

3c lnð1=R0 Þ  x2 2d d

ð24Þ

where d is the width of the PML layer and R0 is the reflection coefficient, which was carefully selected in our experiment; x is the horizontal coordinate with origin chosen at the interface between the PML layer and interior region. ry has the same profile. The point source is a sinusoidal pulse modulated by a Gaussian shaped envelope, commonly used by ultrasound imaging systems [9], given by

SðtÞ ¼  cosð2pftÞ exp 

2pðt  t 0 Þ2

!

s2

where f = 106 Hz, t0 = 0.8s, t = 3  106.

ð25Þ

where u(i, j) is the experimental value of the acoustic pressure at the grid (i, j), N is the number of grid points in each coordinate direction. The relative error R at each time step is also calculated according to the following expression

 n   u  unref   R ¼ 20log10  maxðuref Þ

ð27Þ

where max() represents the maximum value over all time steps. The reference solution uref is calculated by extending the computational domain large enough (1200  1200) so that there is no reflection from its outer boundaries during the time-stepping span of interest. We compare the results of our method with those of the classic split PML in [12,13]. Figs. 3 and 4 illustrate the error relative to the maximum reference field as a function of the time at receivers A and B respectively for various PML thicknesses as computed via the proposed method and the original one. A comparison of the maximum error of our method and the original one for different PML thickness at the receiver A is shown in Table 1. It can be found that both methods lead to a low relative error which confirms their excellent absorbing efficiency. Obviously, their results are quite comparable to each other. It should also be noted that the relative error can be reduced by simply adding extra layers. We can see that when the layer number rises from 30 to 50, the largest relative error improves nearly 10 dB. However, it is also worthwhile to mention that the adding of PML layers does not contribute significantly to the improvement of absorbing efficiency when the layers of PML are large enough. This is because at this moment, the relative error is mainly determined by the discretization error of the equation.

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

5

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

Fig. 2. Snapshots of the pressure propagation in the homogeneous computational domain with 70-cell PML at four different time steps: (a) 150, (b) 300, (c) 400 and (d) 700.

Fig. 3. The relative error at the receiver A in the homogeneous computational domain via the proposed method (a) and the split PML (b) as a function of the time.

Fig. 4. The relative error at the receiver B in the homogeneous computational domain via the proposed method (a) and the split PML (b) as a function of the time.

As the point source generates a certain amount of energy into the calculation domain, the energy will be absorbed by the PML and the total energy will decay as the time passes. Fig. 5 displays the evolution of energy as a function of the time in the calculation domain for the proposed method and the split PML with a 50-cell thick PML. The figure indicates that the total energy has decayed dramatically at the latter time after acoustic waves enter the

Table 1 The maximum error in dB for different PML thickness at the receiver A in the computational domain. M

30

50

70

Split PML Proposed method

47.53 dB 47.58 dB

57.07 dB 57.08 dB

60.21 dB 60.53 dB

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

6

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

PML. Again, it can be clearly found that the results of two methods coincide with each other, which demonstrates their comparable absorbing efficiency. 3.2. Validation of absorbing efficiency in layered media In this experiment, we validate the absorbing efficiency in layered media. The domain is depicted in Fig. 6. It has an area of 2.31  2.31 cm, and is divided into 150  150 grid elements with uniform cell spacing. An addition of 50-cell thick PML is placed on its all four sides with the same cell spacing. S is the point source and A, B are receivers. With the very left bottom grid set as (1, 1) and the cell spacing as unit length, their grid locations are indicated in Fig. 6. The size of the figure is shown with the cell spacing as the unit length in this figure. A time step of 50 ns is used and the speed of the sound c is defined as

8 > < 700 m=s; for 1 6 y 6 100; cðx; yÞ ¼ 900 m=s; for 101 6 y 6 150; > : 1300 m=s; for y > 150: Fig. 5. The energy decay in the homogeneous computational domain using the split PML and the proposed one as a function of time.

ð28Þ

In the simulation, R0 is chosen as 1  105 for both the proposed method and the split PML. Fig. 7 illustrates the snapshots in the computational domain at three different time steps. In this figure, we can hardly identify any reflection, which confirms the good absorbing efficiency of the proposed PML. To further evaluate the performance of the proposed PML formulation, we also calculate the relative error of pressure at receivers A and B against a reference solution obtained using an enlarged domain with 950  950 grids. Figs. 8 and 9 are relative error time histories for A, B respectively. We can again see from these two figures that the results of the proposed PML overlays the results of the Split PML, which confirms their comparable absorbing efficiency. The largest relative error of A and B is around 56.4 dB and 60.2 dB respectively, which are quite small. 3.3. Implementation costs and discussion

Fig. 6. The layered computational domain. The domain has three layers of different speed of sound c. The total domain has an area of 2.31  2.31 cm and is discretized into 150  150 grids with an additional of M-cell thick PML. The size of the figure is shown with the cell spacing as the unit length.

In the above experiments, we can notice that the proposed unsplit PML performs equally with the split PML in terms of absorbing efficiency. Since they are equivalent to each other in the frequency domain, the comparative results are theoretically expected. However, their different mathematical expression in time domain does lead to significant difference in the implementation costs. The proposed unsplit formulation can actually give rise to great reduction in extra memory and easier implementation compared with the original split formulation. Take a domain of 90  90 grids with an additional of 30 cells of PML as an example. The computational costs are illustrated in Table 2 for both

Fig. 7. Snapshots of the pressure propagation in the layered computational domain at three different time steps: (a) 400, (b) 900 and (c) 1500.

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

Y. Ma et al. / Ultrasonics xxx (2014) xxx–xxx

7

only represent a small fraction of the total simulation and that it does absorbs the waves almost perfectly. 4. Conclusions

Fig. 8. The relative error at the receiver A in the layered computational domain using the split PML and the proposed one.

In this article, we presented a novel unsplit-field PML for the acoustic wave equation written as second-order system based on the ADE technique. The proposed PML can be used as an ABC in numerical simulations based on the second-order acoustic wave equations such as the Westervelt equation [9,10] and the full-wave equation [11]. It only needs three first-order auxiliary differential equations in each axis, which facilitates its implementation greatly. Numerical tests show that the proposed scheme has almost the same absorbing efficiency as the classical split PML [12], yet it requires far less storage with easier implementation. Finally, we want to mention that though the proposed secondorder implementation of the PML was only applied to the acoustic wave-equations in this article, it can be extended to other secondorder equations as well. Acknowledgements This work is supported by The National Natural Science Foundation of China (81101049 and 61271071), Shanghai Pujiang Talent Program (12PJ1401200), Doctoral Fund of Ministry of Education (20110071 120019), and Fudan University ASIC and System State Key Laboratory Project (11MS007). References

Fig. 9. The relative error at the receiver B in the layered computational domain. using the split PML and the proposed one.

Table 2 Implementation costs of PMLs for a domain of 90  90 grids with an additional 30-cell PML.

Implementation time Total extra variables

Proposed method

Split PML

29.97 s 54,000

36.44 s 158,000

methods. Here the implementation time is the calculating time elapsed in the PML region only, and the total extra variables are total extra variables needed in the PML domain. As we can see from Table 2, the proposed method is faster than the original one by about 17%. What is more, it only needs 34.18% storage of the latter. We also remind that the larger the non-PML domain compared with the PML region is, the more storage our method can save. In addition, we want to mention that both methods are stable over 15,000 time steps. However, the proposed method has much simpler implementation and requires far less extra memory. It can also be found that the PML for second-order acoustic equation needs more cells than that for the first-order system to attenuate propagating waves satisfactorily in the FDTD simulation. However, this is acceptable in practical situations considering that the PML

[1] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994) 185–200. [2] X.J. Yuan, D. Borup, J. Wiskin, M. Berggren, Simulation of acoustic wave propagation in dispersive media with relaxation losses by using FDTD with relaxation losses with PML absorbing boundary condition, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (1999) 14–23. [3] Q.H. Liu, J. Tao, The perfectly matched layer for acoustic waves in absorptive media, J. Acoust. Soc. Am. 102 (1997) (1997) 2072–2082. [4] A. Giannopoulos, Unsplit implementation of higher order PMLs, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 60 (2012) 1479–1485. [5] W.C. Chew, W.H. Weedon, A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates, Microw. Opt. Technol. Lett. 7 (1994) 599–604. [6] T. Yu, B.H. Zhou, B. Chen, An unsplit formulation of the Berenger’s PML absorbing boundary condition for FDTD meshes, IEEE Microw. Wireless Compon. Lett. 13 (2003) 348–350. [7] S.D. Gedney, B. Zhao, An auxiliary differential equation formulation for the complex-frequency shifted PML, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58 (2010) 838–847. [8] O. Ramadan, Auxiliary differential equation formulation: an efficient implementation of the perfectly matched layer, IEEE Microw. Wireless Compon. Lett. 13 (2003) 69–71. [9] A. Karamalis, W. Wein, N. Navab, Fast ultrasound image simulation using the Westervelt equation, MICCAI 2010, Part I, LNCS 6361, 2010, pp. 243–250. [10] J. Huijssen, A. Bouakaz, M.D. Verweij, N. de Jong, Simulations of the nonlinear acoustic pressure field without using the parabolic approximation, IEEE Sympos. Ultrason. 2 (2003) 1851–1854. [11] G.F. Pinton, J. Dahl, S. Rosenzweizg, G.E. Trahey, A heterogeneous nonlinear attenuating full-wave model of ultrasound, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 60 (2009) 1479–1485. [12] D. Komatitsch, J. Tromp, A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation, Geophys. J. Int. 154 (2003) 146–153. [13] X. Li, PML condition for the numerical simulation of acoustic wave, in: 2010 International Conference on Computing, Control and Industrial Engineering vol. 2, 2010, pp. 129–132. [14] R. Matzen, An efficient finite element time-domain formulation for the elastic second-order wave equation: a non-split complex frequency shifted convolutional PML, Int. J. Numer. Meth. Eng. 88 (2011) 951–973. [15] Y.F. Li, O.B. Matar, Convolutional perfectly matched layer for elastic secondorder wave equation, J. Acoust. Soc. Am. 127 (2010) 1318–1327. [16] Q.H. Liu, The PSTD algorithm: a time-domain method requiring only two cells per wavelength, Microw. Opt. Technol. Lett. 15 (1997) 158–165. [17] F. Collino, T. Chrysoula, Application of the perfectly matched layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66 (2001) 294–307.

Please cite this article in press as: Y. Ma et al., A novel unsplit perfectly matched layer for the second-order acoustic wave equation, Ultrasonics (2014), http://dx.doi.org/10.1016/j.ultras.2014.03.016

A novel unsplit perfectly matched layer for the second-order acoustic wave equation.

When solving acoustic field equations by using numerical approximation technique, absorbing boundary conditions (ABCs) are widely used to truncate the...
1MB Sizes 0 Downloads 3 Views