one of the great challenges for wave propagation is the efficient and stable computation of waves in unbounded domains .
the crucial point for these computations is that the numerical scheme avoids any reflections at the boundaries , even in case the diameter of the computational domain is just a fraction of a wavelength . since the eighties of the last century
, several numerical techniques have been developed to deal with this topic : infinite elements , dirichlet - to - neumann operators based on truncated fourier expansions , absorbing boundary conditions , etc .
the advantages and drawbacks of these different approaches have been widely discussed in literature , see e.g. . especially higher order absorbing boundary conditions ( abcs ) have gained increasing interest , since new methods do not involve high order derivatives [ 6,2124,26,27,25,28,8 ] .
an alternative approach to approximate free radiation is to surround the computational domain by an additional damping layer and guarantee within the formulation , that no reflections occur at its interface with the computational domain .
this so - called perfectly matched layer ( pml ) technique was first introduced by berenger   using a splitting of the physical variables and considering a system of first order partial differential equations ( pdes ) for electromagnetics . since then , there has been much research work on this technique which subsequently was applied to different pdes . in the framework of time - harmonic wave propagation , the pml can be interpreted as a complex - valued coordinate stretching .
therewith , a pml formulation for a linear pde in frequency domain can be considered as a straightforward approach . however , in time domain most pml formulations require a first order hyperbolic system , e.g. , .
the difficulty arising for the second order wave equation in time domain is , that an inverse fourier transform of its frequency representation will lead to convolution integrals , see e.g. . a method to avoid convolution integrals
e.g. , in   a pml method for the second - order elastodynamic equations has been formulated .
the basic idea of the formulation is to decompose the gradient operator in terms of components perpendicular and parallel to the interface , and then split the mechanical displacement in four variables .
however , as noticed in   the resulting equations need special treatment for the time stepping and additional memory is needed for the split - field variables . furthermore , such split - field pml methods suffer from numerical instability , see e.g. .
once a pml formulation has been obtained , the question of stability arises , which is a topic of strong ongoing research .
a stability analysis is not trivial and in general it has to be performed for each new formulation .
several works have analyzed the properties of the pml technique , such as   among others .
e.g. , in   a time - domain analysis of pml methods for wave equations in 2d by using the cagniard - de hoop method has been presented .
the main result is to validate the modified fundamental solution extended to the absorbing layers .
this method is easily applicable to the wave equation with any time - dependent point source .
however , the evaluation is not easy for general initial value problems of the wave equation , because those in general include not only propagating but also evanescent waves .
our stability analysis investigates the evolution of the energy over time and we are able to show decay of an upper bound on the energy for our formulation , thus achieving long term stability .
it should be noted that the pml formulation under consideration has been first published in , where a finite difference time domain ( fdtd ) method on space  time staggered grids with explicit time stepping has been used . in  
we have discussed this formulation in a finite element ( fe ) setting and have shown numerically that omitting terms in the time domain formulation ( calling it an almost pml ) provides long term stability when applying it to a computational setup with a thin damping layer . in this contribution , we further investigate in our fe formulation , use unstructured grids with tetrahedra , apply an implicit newmark time stepping method and provide a detailed discussion of our long time stability proof based on energy analysis .
furthermore , we perform rigorous numerical tests with different damping layer thicknesses and damping functions .
the results show that for a thick enough damping layer long time stability is achieved .
a key novel feature of our approach is that the analysis not only works for a special space discretization but applies to arbitrary space discretizations such as geometrically flexible unstructured finite element grids as demonstrated in the numerical case study and the application example of flow induced sound .
numerical case studies and an application to the computation of sound generated by the flow around a side view car mirror are presented in section 4 , and finally we summarize our main achievements in section 5 .
we consider the time dependent wave equation in an unbounded three dimensional rectangular domain dr3(1)1c2p2t2-p=0with the speed of sound c>0 and zero initial conditions . in order to model free radiation
, we surround the propagation region of interest prop by a pml of thickness li in each coordinate direction i and denote this region by pml . to obtain the correct pml formulation of the wave equation within pml
first we transform ( 1 ) to the frequency domain(2)jc2p-p=0with p the fourier transform of p and  the angular frequency . according
to   we introduce the complex change of variables inside pml(3)xi(xi)=xi+1j0xii(x)dx;xi{x , y , z},where the damping function i is positive inside pml and vanishes in prop .
hence , we obtain the following relationsxixi=1+ij=iandxi=1ixi.therewith , the modified helmholtz equation , which has to be solved in pml , reads as(4)xyzjc2p-yzx1xpx-xzy1ypy-xyz1zpz=0with(5)x=1+xj;y=1+yj;z=1+zj.before we proceed , we use the fact that l/xk=0 for kl to rewrite ( 4 ) by(6)xyzjc2p-xyzxpx-yxzypy-zxyzpz=0.now , we investigate in the first term of ( 6 ) and expand the terms i by ( 5 ) to obtain(7)xyzjc2p=jc21+xj1+yj1+zjp=1c2(j)2+j(x+y+z)+xy+xz+yz+xyzjp.for an inverse fourier transformation of ( 7 ) we recognize , that the last term will result in a time integral of p. therefore , we introduce the first auxiliary variable according to(8)v = pj.in a next step , we analyze the second term in ( 6 ) and start with(9)ii=yzxpx=1+yj1+zj1+xjpxnow , we perform the following rearrangements and use also ( 8)(10)ii=(y+j)(z+j)j(x+j)px=yzj+(y+z)+j(x+j)pxxx+jpx=yzj+(y+z-x)(x+j)px+px=1x+jyzvx+(y+z-x)px+px.the same procedure is performed on the third and fourth term in ( 6 ) .
having an inverse fourier transform in mind , we introduce a vectorial auxiliary variable u with the following relations(11)ux=1x+jyzvx+(y+z-x)pxuy=1y+jxzvy+(x+z-y)pyuz=1z+jxyvz+(x+y-z)pz.now , we are ready to apply the inverse fourier transform to ( 6 ) , ( 8) , ( 11 ) and achieve at the following coupled system of partial differential equations(12)1c22pt2+pt+p+v-p-u=0(13)ut+au+bp - cv=0(14)vt = pwith(15)=x+y+zc2;=xy+xz+yzc2;=xyzc2,(16)a=x000y000z;c=yz000xz000xy,(17)b=x-y-z000y-x-z000z-x-y .
the proper choice of the damping functions is of great importance , especially in order to obtain a very robust and efficient pml - technique . assuming a layer thickness of l , e.g. , in y - direction
, the damped wave will be totally reflected at the outer boundary of the pml - region and this reflected wave at the interface between propagation and pml region takes the value(18)pr=p0e-(2/c)cos0ly(y)dy = p0r.a reasonable choice of the reflection factor r has to be chosen by a trade - off between the necessity of sufficient reduction of reflected waves according to ( 18 ) and possible disturbances of the numerical solution by a too rapid damping in a narrow pml region . in our computations
in addition , in order to get rid of the dependence of the overall damping on the speed of sound c , we will choose y directly proportional to c. using ( 18 ) , we obtain the following relations:constant dampingyconst =- clnr2lcos.inverse distance dampingyinverse = cl - y .
as shown in , the inverse distance damping function leads to an optimal damping behavior for the acoustic wave equation in the frequency domain , and compared to other damping functions does not need specification of the reflection coefficient r. this is also confirmed by many of our computations for the wave equation in time .
however , we have also obtained instable results , especially for thin damping layers . in section 4
we provide rigorous numerical tests for different damping layer thicknesses and the two above defined damping functions .
for deriving boundedness of solutions ( p , v , u ) of ( 12)(14 ) in an appropriate norm ( related to the acoustic energy ) we proceed as follows:1.test the system with appropriate multipliers to derive energy estimates;2.combine these estimates to assess the time evolution of a scalar valued function (t ) , more precisely , to show that (t ) is nonincreasing over time ;  can be interpreted as a lyapunov function for the system ( 12)(14);3.prove that by a proper choice of the parameters defining  , the energy can be bounded by a fixed multiple of . test the system with appropriate multipliers to derive energy estimates ; combine these estimates to assess the time evolution of a scalar valued function (t ) , more precisely , to show that (t ) is nonincreasing over time ;  can be interpreted as a lyapunov function for the system ( 12)(14 ) ; prove that by a proper choice of the parameters defining  , the energy can be bounded by a fixed multiple of . in the course of this derivation , it will turn out that instabilities may emerge in regions where c is strictly positive .
thus we will end up with considering a reduced pml ( rpml ) where we just set c0 .
we consider the bounded computational domain  , which includes the propagation region prop and the attached pml layer pml . for the weak formulation ,
we introduce appropriate test functions ,q and  , and describe space and time dependence by considering p(t),v(t ) and u(t ) as elements of the function space v and h , respectively , for each time t0 .
we assume homogeneous dirichlet boundary conditions on p(t),v(t ) and the test functions ,q , but not on u(t),.p(t)==0on,pt(t)=v(t)=q=0on,i.e.,p(t),,v(t),qv = h01(),u(t),h = l2()3t0.in addition , we introduce the matrix d(19)d=y+z000x+z000x+y , so that we can express b by a - d .
therewith , the weak formulation of our coupled system of pdes ( 12)(14 ) becomes(20)2pt2(t)dx+pt(t)dx+p(t)dx+v(t)dx+p(t)dx+u(t)dx=0v,(21)ut(t)dx+(au(t))dx+(bp(t))dx-(cv(t))dx=0h,(22)vt(t)qdx-p(t)qdx=0qv , for all t0 .
the relation according to ( 22 ) will be used explicitly as v(t)/t = p(t ) .
in here ,,0 are scalars , a , b , cd are diagonal matrices b = a - d with a , c , d0 (    meaning positive semidefinite here ) and supp{}pml , {,,,a , b , c , d } ( see ( 15)(17 ) ) .
now we insert different test functions into ( 20 ) and ( 21 ) , and integrate with respect to time , using the fact that , e.g.0tp(s)pt(s)dxds=120tddt|p(s)|2dxds=12|p(t)|2dx-12|p(0)|2dx.setting(23)=p(s)/tin(20),=p(s)in(20)(with some possibly space dependent factor0c2),=p(t)in(21),=fu(t)in(21)(with some possibly space dependent diagonal matrix0f),and using the notation2=||2dx;2=||2dx,1,2=12dx;1,2=12dx , we therewith obtain(24)121cpt(t)2 - 121cpt(0)2+0tpt(s)2ds+12p(t)2 - 12p(0)2+0tv(s),pt(s)ds+12p(t)2 - 12p(0)2+0tu(s),pt(s)ds=0,(25)0t1c22pt2(s),p(s)ds+12p(t)|2 - 12p(0)|2+0tp(s)2ds+12v(t)2 - 12v(0)2+0tp(s)2ds+0tu(s),p(s)ds=0,(26)0tut(s),p(s)ds+0t(au(s)),p(s)ds+0ta1/2p(s))2ds-0td1/2p(s))2ds-12c1/2v(t))2 + 12c1/2v(0))2=0,(27)12f1/2u(t)2 - 12f1/2u(0)2+0t(fa)1/2u(s)2ds+0tfbp(t),u(s)ds-0tfcv(t),u(s)ds=0,for all t0 .
2nd step : conclude monotonicity of some lyapunov function . in a next step , we will compute the sum of the four above equations ( 24 ) , so that the term 0tu(s)/t,p(s)ds cancels out .
furthermore , we define the following function of time(29)(t)121cpt(t)2+p(t)2+f1/2u(t)2++p(t)2+v(t)2-c1/2v(t))2 + 2v(t),p(t)+2u(t),p(t)+2c2pt(t),p(t).remark  1here we wish to point out that the sixth term on the right hand side in ( 29 ) suggests that it is favorable to set c=0 from a stability point of view .
indeed , even if we would add the time integral of ( 20 ) with =v(s ) , which would give us a term 12v(t)2 on the left hand side for possible control of -c1/2v(t))2 , we would end up with a term 0tu(s),v(s)ds , which we can not control since we have no means for obtaining an estimate of 0tv(s)2ds , especially for long times when this expression gets larger than sups[0,t]v(s)2 . here
we wish to point out that the sixth term on the right hand side in ( 29 ) suggests that it is favorable to set c=0 from a stability point of view . indeed
, even if we would add the time integral of ( 20 ) with =v(s ) , which would give us a term 12v(t)2 on the left hand side for possible control of -c1/2v(t))2 , we would end up with a term 0tu(s),v(s)ds , which we can not control since we have no means for obtaining an estimate of 0tv(s)2ds , especially for long times when this expression gets larger than sups[0,t]v(s)2 . using ( 28 ) and the above definition according to ( 29 ) , we obtain the following relation for the sum of ( 24)(27)(30)(t)+0t(fa)1/2u(s)2ds+0t-c2pt(s)2ds+0tsign(-)||-|p(s)|2dxds+0tsign(i+b)||i+b|p(s)|2dxds=(0)-0t(i+a+fb)p(s),u(s)ds+0tfcv(t),u(s)ds(0)+120t(fa)1/2u(s)2ds+120t((fa)1/2)(i+a+fb)p(s)2ds+120t((fa)1/2)fcv(s)2dsor if c=0(31)(t)+0t(fa)1/2u(s)2ds+0t-c2pt(s)2ds+0tsign(-)||-|p(s)|2dxds+0tsign(i+b)||i+b|p(s)|2dxds=(0)-0t(i+a+fb)p(s),u(s)ds(0)+0t(fa)1/2u(s)2ds+0t12((fa)1/2)(i+a+fb)p(s)2ds , provided(32)u(s)n(fa),s(0,t)orr(i+a+fb)n(fa)(and for ( 30 ) additionally u(s)n(fc),s(0,t ) ) . in here , for some matrix m , we denote by r(m ) , n(m) , and m the range , the orthogonal complement of the nullspace , and the generalized inverse , respectively , and for some positive semidefinite matrix m we denote by m1/2 the square root , being defined by the relation m1/2m1/2=m . for deriving ( 30 ) and ( 31 )
, we have exploited the identity0t(i+a+fb)p(s),u(s)ds=0t((fa)1/2)(i+a+fb)p(s),((fa)1/2)u(s)ds , that holds under condition ( 32 ) , and whose right hand side can then be estimated via the cauchy 
schwarz inequality . studying ( 30 ) as well as ( 31 ) , we see that for(33)c=0,0c2,-0,i+b0,and(34)4fa(i+b)-(i+a+fb)20 in all of  , we get(t)(0).note that since here time  =  0 can be replaced by time = st , this shows monotone decrease of the function . 3rd step : estimate energy by . reconsidering the definition of (t ) 
 ( 29 ) , we see that in case c=0 it is an upper bound for a multiple of the energy(35)(t)121cpt(t)2+e1/2p(t)2+e1/2u(t)2with 0,e0,e0 , provided that additionally to ( 33 ) and ( 34)(36)i - f-10and1-cp2+2c2-(+)i - f-10,where cp is the constant in the poincar 
friedrichs inequalityv = h01():cp.namely , we can estimate(t)=121cpt(t)2+p(t)2+f1/2u(t)2++p(t)2+v(t)2 + 2v(t),p(t)+2u(t),p(t)+2c2pt(t),p(t)121cpt(t)2+p(t)2+f1/2u(t)2++p(t)2+v(t)2-v(t)2-p(t)2-(f - e)1/2u(t)2-(f - e)-1/2p(t)2 - 1-cpt(t)2-c1-p(t)2121cpt(t)2+e1/2p(t)2+e1/2u(t)2,where the last two inequalities hold if(37)f - e>0andi - e-(f - e)-10and+2c2(1-)-(+)p(t)2i - e-(f - e)-1p(t)i.e . , if(38)f - e>0andi - e-(f - e)-1-cp2+2c2(1-)-(+)i0 .
it remains to interpret conditions ( 32)(34 ) , ( 36 ) , ( 38 ) i.e. , to find ,f,,e,e such that they can be satisfied : according to ( 15)(17 ) we havec2=x+y+z , c2=yz+zx+xy , c2=xyz , aii=i , bii=i-j-k , dii=j+k;i,j,k{x,y,z}and setting c=0 , conditions ( 33 ) are satisfied for=x+y+z , sincec2(-)=(yz+zx+xy)(x+y+z)-xyz0.since f will be chosen as a diagonal matrix , condition ( 34 ) can be reformulated as04fa(i+a - d)-((i+a)+f(a - d))2=4(i+a)af-4adf-((i+a)+af - df)2=4(i+a)af-4adf-(i+a)2-a2f2-d2f2 - 2(i+a)af+2(i+a)df+2adf2=-(i+a)2+(2(i+a)(a+d)-4ad)f-(a - d)2f2,where the diagonal entries of the right hand side are ( with the abbreviation a = aii=i , d = dii=j+k , f = fii,=a+d)-(2a+d)2+(2(2a+d)(a+d)-4ad)f-(a - d)2f2=-(2a+d)2 + 2(2a2+ad+d2)f-(a - d)2f20f-ff+where(39)f=2a2+ad+d2(2a2+ad+d2)2-(2a+d)2(a - d)2(a - d)2.note that in here , the square root gives a real number since2a2+ad+d2|(2a+d)(a - d)|=|2a2-ad - d2|and we have 0f-f+ . therefore , f = fii0 can indeed be chosen so that ( 34 ) is satisfied .
condition ( 36 ) additionally requires f = fii1 which is enabled by the fact that f+1 due to(2a2+ad+d2)2-(2a+d)2(a - d)2(a - d)2-(2a2+ad+d2)=-a2 - 3ad , with strict inequality unless a vanishes , hence(40)2a2+ad+d2+(2a2+ad+d2)2-(2a+d)2(a - d)2(a - d)2,with strict inequality if a>0 .
( in case a = d=0 , condition ( 34 ) is trivially satisfied for any f and we can choose an arbitrary f > i . ) to see the right hand part of ( 36 ) we use the fact that we have chosen =c2 and ( 33 ) in(41)(+)--2c2=-=-0,so that the right hand part of ( 36 ) is implied by the left hand part . concerning ( 38 ) ,
we get from estimate ( 41 ) that it follows from(42)f - e0andi - e-(f - e)-1-cp22c21-i0 . in order to check the nullspace condition ( 32 ) and to identify regions
in which the matrices e,e or at least some of their entries can be chosen strictly positive , we distinguish between the following subdomains:(a)0={:x()=y()=z()=0}prop : on this subdomain the quantities ,,,a , b , c vanish .
since ( 13 ) together with u(0)=0 implies u(t)0=0 for all t>0 , we directly obtain that the respective part of (t ) , namely ( 29 ) with 2 replaced by 0l2(0)2 , just equals 0(t)121cpt(t)02+p(t)02 , so we can choose1,ei , ei , on0.(b)xyz={:x()>0andy()>0andz()>0}:since on this subdomain a is strictly positive definite , choosing a positive f we will automatically satisfy the nullspace condition ( 32 ) .
we can select f = diag(f11,f22,f33 ) with fii = f+=fi+ according to ( 39 ) with a=i>0 which implies fii>1 ( cf .
the observation following ( 40 ) ) . with this choice , =x+y+z>0 ande=12(f - i)>0,e=14i+12(f - i)-1(f - i)>0,=1+,onxyz , with 14c2cp22mini{1,2,3}(1
+ 12(fii-1))-1(fii-1)>0 , we also get ( 42 ) , which by ( 41 ) is sufficient for ( 38).(c)x={:x()>0andy()=z()=0 } ( analogously for y,z):here it can be seen that the nullspace condition ( 32 ) indeed restricts the first component of f : since n(fa)n(a)t = r{0}{0 } , which in general will not contain u(s ) , we have to make sure that ( i+a+fb)r{0}{0 } , which by fii1 ( from ( 36 ) ) , x ( from ( 33 ) ) is equivalent to f22=f33=1,=x>0 , so that we have to pute22=e33=e22=e33=0onx.still , application of the argument from ( b ) and direct consideration of ( 38 ) shows that we can choosee11=12(f11 - 1)>0,e11=141 + 12(f11 - 1)-1(f11 - 1)>0,=1+>0,onx , with 14c2cp22(1 + 12(f11 - 1))-1(f11 - 1)>0.(d)yz={:x()=0andy()>0andz()>0 } ( analogously for zx,xy):this case can be treated similarly to ( c ) , which yields a choice f11=1,fjj>1,=y+z>0,e11=e11=0onx.ejj=12(fjj-1)>0,ejj=141 + 12(fjj-1)-1(fjj-1)>0,=1+>0,onyz , for j=2,3 , with 14c2cp22mini{2,3}(1 + 12(fii-1))-1(fii-1)>0 .
since ( 13 ) together with u(0)=0 implies u(t)0=0 for all t>0 , we directly obtain that the respective part of (t ) , namely ( 29 ) with 2 replaced by 0l2(0)2 , just equals 0(t)121cpt(t)02+p(t)02 , so we can choose1,ei , ei , on0 . xyz={:x()>0andy()>0andz()>0 } : since on this subdomain a is strictly positive definite , choosing a positive f we will automatically satisfy the nullspace condition ( 32 ) .
we can select f = diag(f11,f22,f33 ) with fii = f+=fi+ according to ( 39 ) with a=i>0 which implies fii>1 ( cf .
the observation following ( 40 ) ) . with this choice , =x+y+z>0 ande=12(f - i)>0,e=14i+12(f - i)-1(f - i)>0,=1+,onxyz , with 14c2cp22mini{1,2,3}(1
+ 12(fii-1))-1(fii-1)>0 , we also get ( 42 ) , which by ( 41 ) is sufficient for ( 38 ) .
x={:x()>0andy()=z()=0 } ( analogously for y,z ) : here it can be seen that the nullspace condition ( 32 ) indeed restricts the first component of f : since n(fa)n(a)t = r{0}{0 } , which in general will not contain u(s ) , we have to make sure that ( i+a+fb)r{0}{0 } , which by fii1 ( from ( 36 ) ) , x ( from ( 33 ) ) is equivalent to f22=f33=1,=x>0 , so that we have to pute22=e33=e22=e33=0onx.still , application of the argument from ( b ) and direct consideration of ( 38 ) shows that we can choosee11=12(f11 - 1)>0,e11=141 + 12(f11 - 1)-1(f11 - 1)>0,=1+>0,onx , with 14c2cp22(1 + 12(f11 - 1))-1(f11 - 1)>0 .
yz={:x()=0andy()>0andz()>0 } ( analogously for zx,xy ) : this case can be treated similarly to ( c ) , which yields a choice f11=1,fjj>1,=y+z>0,e11=e11=0onx.ejj=12(fjj-1)>0,ejj=141 + 12(fjj-1)-1(fjj-1)>0,=1+>0,onyz , for j=2,3 , with 14c2cp22mini{2,3}(1 + 12(fii-1))-1(fii-1)>0 .
thus , we have proven the following stability result for a modified version of the pml system ( 12)(14):theorem  1consider the system 
 ( 12)(14 ) 
 with homogeneous dirichlet boundary conditions on p , homogeneous initial conditions on 
 u , and 
 c,,,,a , b , according to 
 ( 15)(17 ) , in 
 l(),l()33 , respectively , but with 
 c0.there exist a space dependent strictly positive function 
 >0 
 and positive semidefinite diagonal matrix functions 
 e,e 
 as well as a nonnegative monotonically non - increasing time dependent function 
  
 such that(t)121cpt(t)2+e1/2p(t)2+e1/2u(t)2holds for any solution 
 ( p , v , u)w(r+;v , h ) 
 of 
 ( 20)(22 ) , wherew(i , v , h)=(c2(i;v)c1(i;l2())c(i;v))(c2(i;l2())c1(i;v))c(i;h),v = h01(),h = l2()3,v=h-1( ) 
 ( the dual of 
 v).more precisely we have 
 e()=diag(e11(),e22(),e33( ) 
 and 
 e()=diag(e11(),e22(),e33( ) ) 
 with ( upon identification 
 1x,2y,3z)(43)eii>0,eii>0oni , where(44)i={:i()>0}{:x()=y()=z()=0 } consider the system 
 ( 12)(14 ) 
 with homogeneous dirichlet boundary conditions on p , homogeneous initial conditions on 
 u , and 
 c,,,,a , b , according to 
 ( 15)(17 ) , in 
 l(),l()33 , respectively , but with 
 c0 .
there exist a space dependent strictly positive function 
 >0 
 and positive semidefinite diagonal matrix functions 
 e,e 
 as well as a nonnegative monotonically non - increasing time dependent function 
  
 such that(t)121cpt(t)2+e1/2p(t)2+e1/2u(t)2holds for any solution 
 ( p , v , u)w(r+;v , h ) 
 of 
 ( 20)(22 ) , wherew(i , v , h)=(c2(i;v)c1(i;l2())c(i;v))(c2(i;l2())c1(i;v))c(i;h),v = h01(),h = l2()3,v=h-1( ) 
 ( the dual of 
 v ) .
more precisely we have 
 e()=diag(e11(),e22(),e33( ) 
 and 
 e()=diag(e11(),e22(),e33( ) ) 
 with ( upon identification 
 1x,2y,3z)(43)eii>0,eii>0oni , where(44)i={:i()>0}{:x()=y()=z()=0 } note that it suffices to impose homogeneous initial conditions on u in 0={:x()=y()=z()=0 } only.remark  2here we have assumed that a weak solution ( p , v , u)w(r+,v , h ) of ( 12)(14 ) exists .
indeed , existence and uniqueness at a first glance should be establishable for c0 along the usual lines ( cf .
, e.g. , ) : 1.galerkin approximation , i.e. , replacing vvh by a sequence of finite dimensional nested subspaces vnvnhn whose union is dense in vvh , and which satisfy hn for all vn ( cf .
( 23 ) ) ; local in time existence of the galerkin solutions ( pn , vn , un)w([0,t],v , h ) follows from linearity of the system;2.uniform energy estimates of the galerkin solutions in w([0,t];v , h ) , which can be done as above , replacing the test and ansatz spaces by their finite dimensional counterparts vn , hn , using the topology of the infinite dimensional spaces v , h on them ; the additional c(0,t;v ) bound on 2pnt2 follows directly from ( 20).3.these uniform estimates imply existence of weakly convergent subsequence whose weak limits by linearity solve the system ( 20)(22 ) in a weak sense provided the coefficients are in l( ) ; global in time existence follows from the uniform energy estimates.however , there is a gap in this proof : since , e.g. , e22=e33=e22=e33=0 in x={:x()>0andy()=z()=0 } , we only get boundedness of the l2(x ) norm of the first component of pn and of un , so only this component converges weakly , while the other two components might diverge .
still we get existence of a very weak solution , i.e. , ( p , v , u)(c2(0,t;(h2()))c1(0,t;l2()))c2(0,t;l2())c1(0,t;(h1()3) ) such that(45)2pt2(t)dx+pt(t)dx+p(t)dx+v(t)dx-p(t)dx+u(t)dx=0h2(),(46)ut(t)dx+(au(t))dx-p(t)(b)dx=0h1()3,(47)vt(t)=p(t)in this way.corollary  1let the assumptions of 
 theorem 1 
 be satisfied with additionally 
 bw1,()33 
 and 
 mini{1,2,3}infii()>0 , where 
 i={:i()>0}.then for any 
 t>0 

there exists a solution 
 ( p , v , u)(c2(0,t;(h2()))c1(0,t;l2()))c2(0,t;l2())c1(0,t;(h1()3) ) 
 of 
 ( 45)(47 ) , 
 t(0,t ) , for which a nonnegative monotonically non - increasing time dependent function 
  
 exists such that(t)1cpt(t)2for all 
 t(0,t).proofwe proceed as described in remark 2 with the additional assumption vnh2(),hnh1( ) on the test spaces .
along the lines of the proof of theorem ( 1 ) we obtain uniform boundedness of pntc(0,t;l2())=2vnt2c(0,t;l2( ) ) for the galerkin solutions ( pn , vn , un ) of ( 20)(22 ) . from this and the time differentiated version of ( 21 ) with =unt ( note that c=0 ) we obtain uniform boundedness of untc(0,t;(h1()3) ) .
altogether we have existence of a uniform constant k>0 such that for arbitrarily fixed finite time t:pntc(0,t;l2())=2vnt2c(0,t;l2())k , pnc(0,t;l2())=vntc(0,t;l2())k(1+t),vnc(0,t;l2())k(1+t)2,untc(0,t;(h1()3))kunc(0,t;(h1()3))k(1+t),2pnt2c(0,t;(h2()))k(1+t)2,where we have used ( 45 ) to obtain the latter estimate .
the function  here is set to 1min with  from theorem 1 , where min = min{1,infii } , i=14c2cp22(1 + 12(fii-1))-1(fii-1 ) , fii = min{2,2a2+ad+d2+(2a2+ad+d2)2-(2a+d)2(a - d)2(a - d)2}>1 
 a=i , d=j+k , j , ki,=a+d . 
 here we have assumed that a weak solution ( p , v , u)w(r+,v , h ) of ( 12)(14 ) exists . indeed , existence and uniqueness at a first glance should be establishable for c0 along the usual lines ( cf . , e.g. , ) : 1.galerkin approximation , i.e. , replacing vvh by a sequence of finite dimensional nested subspaces vnvnhn whose union is dense in vvh , and which satisfy hn for all vn ( cf .
( 23 ) ) ; local in time existence of the galerkin solutions ( pn , vn , un)w([0,t],v , h ) follows from linearity of the system;2.uniform energy estimates of the galerkin solutions in w([0,t];v , h ) , which can be done as above , replacing the test and ansatz spaces by their finite dimensional counterparts vn , hn , using the topology of the infinite dimensional spaces v , h on them ; the additional c(0,t;v ) bound on 2pnt2 follows directly from ( 20).3.these uniform estimates imply existence of weakly convergent subsequence whose weak limits by linearity solve the system ( 20)(22 ) in a weak sense provided the coefficients are in l( ) ; global in time existence follows from the uniform energy estimates .
galerkin approximation , i.e. , replacing vvh by a sequence of finite dimensional nested subspaces vnvnhn whose union is dense in vvh , and which satisfy hn for all vn ( cf .
( 23 ) ) ; local in time existence of the galerkin solutions ( pn , vn , un)w([0,t],v , h ) follows from linearity of the system ; uniform energy estimates of the galerkin solutions in w([0,t];v , h ) , which can be done as above , replacing the test and ansatz spaces by their finite dimensional counterparts vn , hn , using the topology of the infinite dimensional spaces v , h on them ; the additional c(0,t;v ) bound on 2pnt2 follows directly from ( 20 ) .
these uniform estimates imply existence of weakly convergent subsequence whose weak limits by linearity solve the system ( 20)(22 ) in a weak sense provided the coefficients are in l( ) ; global in time existence follows from the uniform energy estimates .
however , there is a gap in this proof : since , e.g. , e22=e33=e22=e33=0 in x={:x()>0andy()=z()=0 } , we only get boundedness of the l2(x ) norm of the first component of pn and of un , so only this component converges weakly , while the other two components might diverge .
still we get existence of a very weak solution , i.e. , ( p , v , u)(c2(0,t;(h2()))c1(0,t;l2()))c2(0,t;l2())c1(0,t;(h1()3) ) such that(45)2pt2(t)dx+pt(t)dx+p(t)dx+v(t)dx-p(t)dx+u(t)dx=0h2(),(46)ut(t)dx+(au(t))dx-p(t)(b)dx=0h1()3,(47)vt(t)=p(t)in this way .
let the assumptions of 
 theorem 1 
 be satisfied with additionally 
 bw1,()33 
 and 
 mini{1,2,3}infii()>0 , where 
 i={:i()>0}. then for any 
 t>0 
 there exists a solution 
 ( p , v , u)(c2(0,t;(h2()))c1(0,t;l2()))c2(0,t;l2())c1(0,t;(h1()3) ) 
 of 
 ( 45)(47 ) , 
 t(0,t ) , for which a nonnegative monotonically non - increasing time dependent function 
  
 exists such that(t)1cpt(t)2for all 
 t(0,t ) .
we proceed as described in remark 2 with the additional assumption vnh2(),hnh1( ) on the test spaces . along the lines of the proof of theorem ( 1 ) we obtain uniform boundedness of pntc(0,t;l2())=2vnt2c(0,t;l2( ) ) for the galerkin solutions ( pn , vn , un ) of ( 20)(22 ) . from this and the time differentiated version of ( 21 ) with =unt ( note that c=0 ) we obtain uniform boundedness of untc(0,t;(h1()3) ) .
altogether we have existence of a uniform constant k>0 such that for arbitrarily fixed finite time t:pntc(0,t;l2())=2vnt2c(0,t;l2())k , pnc(0,t;l2())=vntc(0,t;l2())k(1+t),vnc(0,t;l2())k(1+t)2,untc(0,t;(h1()3))kunc(0,t;(h1()3))k(1+t),2pnt2c(0,t;(h2()))k(1+t)2,where we have used ( 45 ) to obtain the latter estimate .
the function  here is set to 1min with  from theorem 1 , where min = min{1,infii } , i=14c2cp22(1 + 12(fii-1))-1(fii-1 ) , fii = min{2,2a2+ad+d2+(2a2+ad+d2)2-(2a+d)2(a - d)2(a - d)2}>1 
 a=i , d=j+k , j , ki,=a+d . 
 according to our stability analysis , we have no estimate on the l2-wrt time norm of ||v|| ( cf .
therefore , we will also numerically investigate the case , where c is set to zero ( i.e. we omit cv in ( 13 ) ) and will denote this reduced formulation by rpml .
we are aware of the fact that this rpml formulation will not achieve perfect matching . however , as numerical results will demonstrate , the additional error compared to the full pml formulation is small , and it strongly increases the stability in case of thin damping layers .
furthermore , we want to state that the rpml is a true pml in case of 1d as well as 2d computations . in these cases , we do not need the additional scalar auxiliary variable v and so c is not present .
e.g. , in 2d we just need the auxiliary vector variable u=(ux , uy)t leading to a total number of just three unknowns in the pml region .
this can be also seen by analyzing ( 12)(14 ) , e.g. , assuming waves in the xy - plane .
then  and c get zero ( z=0 ) resulting in just three scalar equations for p , ux , uy . so an error just occurs in 3d , when waves propagate towards corners , where all three damping coefficients are active .
the setup of our numerical example for investigating the pml formulation is displayed in fig . 1 .
on the surface of the sphere we set a dirichlet boundary condition for all nodes according to the following functionf(t)=ddte-2(f0t-1)2.we choose the propagation region to be /2 (  being the wavelength given by =c / f0 with c the speed of sound ) and use pml regions with /8,/4 and /2 .
the propagation as well as pml region are discretized with tetrahedra finite elements of first order having an average edge length of about /16 .
thus , the pml is discretized in thickness direction with two finite elements for the /8 layer , with four finite elements for the /4 layer , and with eight finite elements for the /2 layer , respectively .
furthermore , we use an implicit newmark time stepping scheme with =0.25 and =0.5 
  and we choose a time step size of t=1/(80f0 ) . for the first investigations concerning the accuracy
, we use the inverse distance damping functions as given in section 2 . for this simple setup
however , to just measure the error introduced by the pml , we compute a reference solution using the same computational mesh in the propagation region and use instead of the pml region an additional large propagation region .
we consider the evolution of the l2-error at each time step tn(48)errorpml(tn)=1npropi=1nproppi(tn)-pnref(tn)2.in ( 48 ) 
 pi denotes the solution at each fe node i using the pml formulation , piref the reference solution at fe node i and nprop the number of nodes within the propagation region . in a first investigation
, we computed the acoustic pressure over time applying the pml formulation and compare the results with our reference solution at the three observation points p1,p2 and p3 .
4 . from the three graphs , we can clearly see the improvement in the solutions with increasing pml thickness . the main difference , as expected , is given for the pressure at observation point p3 , which is the corner point between propagation region and pml region ( see fig .
, the results demonstrate that the reduced pml formulation denoted by rpml just has slightly worse properties than the full pml formulation and the differences vanishes almost completely for a damping layer thickness of /2 ( see fig .
4 ) . based on the obtained results we further computed the l2-norm of the error as given in ( 48 ) . in fig .
5 we show the obtained error introduced by the pml and rpml with a thickness of /8 and /4 . again
, we observe the decrease in the error over time when increasing the pml thickness .
we had to stop the computation of the error at t  =  40  s , since at this time the reference solution already showed first reflections back to the propagation region .
therefore , in order to investigate the stability , we also computed the overall acoustic energy at each time step tn within the propagation region by the following formula(49)eacoustic(tn)=prop02v(tn)v(tn)+p(tn)220c2dx.in ( 49 ) 
 0 denotes the mean density of the fluid and v the acoustic particle velocity . for this long term stability study ,
we perform computations for all three cases of layer thicknesses ( /8,/4 and /2 ) and use the two different damping functions  constant and inverse distance  as introduced in section 2 .
the acoustic particle velocity v is related to the acoustic pressure p via the linear momentum conservation0vt=-p.therewith , in a time discrete setting using a trapezoidal scheme , we obtainv(tn)=v(tn-1)-t0p(tn),which is inserted into ( 49 ) . fig .
6 displays the computed acoustic energy within the propagation region for the reference solution and the pml as well as rpml formulations .
the first main observation is that for the /8 damping layer the results for both damping functions get instable .
canceling the critical term in the pml formulation ( denoted by rpml ) as revealed by the stability proof , leads to stable results ( see fig .
the second main observation is that both damping functions results in stable long time computations for layer thicknesses of /4 and /2 .
furthermore , applying the rpml formulation achieves in all cases a stronger damping behavior of the acoustic energy for increasing time .
furthermore , when comparing the energy curves with the reference solution , we may clearly state that the inverse distance damping function ( both for pml and rpml ) provides the most accurate results . as a practical application of our pml formulation
, we consider the computation of flow induced sound by applying lighthill s analogy , which results in an inhomogeneous wave equation(50)1c22pt2-p=(t).here , t denotes lighthill s tensor , which for incompressible flows and isentropic state can be approximated by(51)tfvvwith v the flow velocity and f the mean density of the fluid . for all details about our formulation
, our focus is on the application of the pml formulation and so we just briefly discuss the flow computation and then the acoustic results .
7 displays the geometry of the side view mirror and the computed flow field for a characteristic time step .
we have computed the flow field by applied a les ( large eddy simulation ) on a grid with 6.5 million hexahedra applying ansys fluent very similarly as reported in . in order to preserve the acoustic energy ,
we have first computed the acoustic sources within our fe formulation on the fine flow grid . in a second step , we have performed a conservative interpolation to the coarser acoustic grid .
8 displays the acoustic field around the side view mirror on two slices for a characteristic time step .
one can clearly see the propagating waves which are then absorbed within the pml region .
we have used a quite coarse acoustic grid with about 340.000 elements having about four elements within the pml thickness and applying an inverse distance damping function . as in our numerical test case ( see section 4 ) , we used a fully implicit newmark time discretization scheme and observed no instabilities for our pml formulation .
therewith , we can state that also for complex computational setups with a broadband acoustic field the pml formulation works effective .
in addition , the few additional auxiliary unknowns within the pml region result in a fast computation . finally , fig .
9 displays the computed sound pressure level ( spl)1 over the frequency range at the indicated monitoring point compared to measurements .
the setup of our numerical example for investigating the pml formulation is displayed in fig . 1 .
on the surface of the sphere we set a dirichlet boundary condition for all nodes according to the following functionf(t)=ddte-2(f0t-1)2.we choose the propagation region to be /2 (  being the wavelength given by =c / f0 with c the speed of sound ) and use pml regions with /8,/4 and /2 .
the propagation as well as pml region are discretized with tetrahedra finite elements of first order having an average edge length of about /16 .
thus , the pml is discretized in thickness direction with two finite elements for the /8 layer , with four finite elements for the /4 layer , and with eight finite elements for the /2 layer , respectively .
furthermore , we use an implicit newmark time stepping scheme with =0.25 and =0.5 
  and we choose a time step size of t=1/(80f0 ) . for the first investigations concerning the accuracy
, we use the inverse distance damping functions as given in section 2 . for this simple setup
however , to just measure the error introduced by the pml , we compute a reference solution using the same computational mesh in the propagation region and use instead of the pml region an additional large propagation region .
we consider the evolution of the l2-error at each time step tn(48)errorpml(tn)=1npropi=1nproppi(tn)-pnref(tn)2.in ( 48 ) 
 pi denotes the solution at each fe node i using the pml formulation , piref the reference solution at fe node i and nprop the number of nodes within the propagation region . in a first investigation
, we computed the acoustic pressure over time applying the pml formulation and compare the results with our reference solution at the three observation points p1,p2 and p3 .
4 . from the three graphs , we can clearly see the improvement in the solutions with increasing pml thickness . the main difference , as expected , is given for the pressure at observation point p3 , which is the corner point between propagation region and pml region ( see fig .
, the results demonstrate that the reduced pml formulation denoted by rpml just has slightly worse properties than the full pml formulation and the differences vanishes almost completely for a damping layer thickness of /2 ( see fig .
4 ) . based on the obtained results we further computed the l2-norm of the error as given in ( 48 ) . in fig .
5 we show the obtained error introduced by the pml and rpml with a thickness of /8 and /4 . again
, we observe the decrease in the error over time when increasing the pml thickness .
we had to stop the computation of the error at t  =  40  s , since at this time the reference solution already showed first reflections back to the propagation region .
therefore , in order to investigate the stability , we also computed the overall acoustic energy at each time step tn within the propagation region by the following formula(49)eacoustic(tn)=prop02v(tn)v(tn)+p(tn)220c2dx.in ( 49 ) 
 0 denotes the mean density of the fluid and v the acoustic particle velocity . for this long term stability study ,
we perform computations for all three cases of layer thicknesses ( /8,/4 and /2 ) and use the two different damping functions  constant and inverse distance  as introduced in section 2 .
the acoustic particle velocity v is related to the acoustic pressure p via the linear momentum conservation0vt=-p.therewith , in a time discrete setting using a trapezoidal scheme , we obtainv(tn)=v(tn-1)-t0p(tn),which is inserted into ( 49 ) . fig .
6 displays the computed acoustic energy within the propagation region for the reference solution and the pml as well as rpml formulations .
the first main observation is that for the /8 damping layer the results for both damping functions get instable .
canceling the critical term in the pml formulation ( denoted by rpml ) as revealed by the stability proof , leads to stable results ( see fig .
the second main observation is that both damping functions results in stable long time computations for layer thicknesses of /4 and /2 .
furthermore , applying the rpml formulation achieves in all cases a stronger damping behavior of the acoustic energy for increasing time .
furthermore , when comparing the energy curves with the reference solution , we may clearly state that the inverse distance damping function ( both for pml and rpml ) provides the most accurate results .
as a practical application of our pml formulation , we consider the computation of flow induced sound by applying lighthill s analogy , which results in an inhomogeneous wave equation(50)1c22pt2-p=(t).here , t denotes lighthill s tensor , which for incompressible flows and isentropic state can be approximated by(51)tfvvwith v the flow velocity and f the mean density of the fluid .
for all details about our formulation we refer to . here , our focus is on the application of the pml formulation and so we just briefly discuss the flow computation and then the acoustic results .
7 displays the geometry of the side view mirror and the computed flow field for a characteristic time step .
we have computed the flow field by applied a les ( large eddy simulation ) on a grid with 6.5 million hexahedra applying ansys fluent very similarly as reported in . in order to preserve the acoustic energy ,
we have first computed the acoustic sources within our fe formulation on the fine flow grid . in a second step ,
8 displays the acoustic field around the side view mirror on two slices for a characteristic time step .
one can clearly see the propagating waves which are then absorbed within the pml region .
we have used a quite coarse acoustic grid with about 340.000 elements having about four elements within the pml thickness and applying an inverse distance damping function . as in our numerical test case ( see section 4 ) , we used a fully implicit newmark time discretization scheme and observed no instabilities for our pml formulation .
therewith , we can state that also for complex computational setups with a broadband acoustic field the pml formulation works effective .
in addition , the few additional auxiliary unknowns within the pml region result in a fast computation . finally , fig .
9 displays the computed sound pressure level ( spl)1 over the frequency range at the indicated monitoring point compared to measurements .
we have introduced an efficient pml formulation for the second order wave equation . by using auxiliary variables ,
we achieve at a time domain formulation without convolution integrals . for our un - split pml formulation
we observed a term within our formulation which we could not bound due to its adverse sign , and which may be a source for long time instabilities .
the rigorous numerical tests revealed that for layer thicknesses of /4 and /2 the pml formulation for both investigated damping functions ( constant and inverse distance ) achieves long time stability .
omitting the critical term , we arrive at our rpml formulation , which strongly improves the stability as demonstrated by the numerical examples .
therewith , we can summarize the three main results as follows : ( 1 ) we can prove long time stability of our formulation ; just one critical term containing the mixed products of the damping functions , may disturb the stability ; omitting this term yields long time stability for arbitrary layer thickness , ( 2 ) numerical tests show that the full pml formulation is long time stable for thick enough damping layers ( in our examples at least /4 ) , and ( 3 ) numerical tests demonstrate that for thin damping layers one should use the rpml formulation to achieve long time stability .
 we consider the second order wave equation in an unbounded domain and propose an advanced perfectly matched layer ( pml ) technique for its efficient and reliable simulation . in doing so , we concentrate on the time domain case and use the finite - element ( fe ) method for the space discretization . 
 our un - split - pml formulation requires four auxiliary variables within the pml region in three space dimensions . for a reduced version ( rpml ) 
 , we present a long time stability proof based on an energy analysis . 
 the numerical case studies and an application example demonstrate the good performance and long time stability of our formulation for treating open domain problems . 
Introduction
PML formulation
Stability analysis
Results
Numerical studies
Practical application: side view car mirror
Conclusion