Some contributions in time-harmonic dissipative acoustic problems
Full text
Andr´es Prieto Aneiros PhD Dissertation SOME CONTRIBUTIONS IN TIME-HARMONIC DISSIPATIVE ACOUSTIC PROBLEMS Departamento de Matem´atica Aplicada Facultade de Matem´aticas
¿C´omo me vas a explicar, di, la dicha de esta tarde, si no sabemos porqu´e fue, ni c´omo, ni de qu´e ha sido, si es pura dicha de nada? En nuestros ojos visiones, visiones y no miradas, no percib´ıan tama˜nos, datos, colores, distancias. Pedro Salinas
Contents Preface v I Porous materials 1 1 Porous models 3 1.1 Introduction .................................... 4 1.2 Rigid porous models ............................... 6 1.2.1 Darcy’s like model ............................ 7 1.2.2 Allard-Champoux model ......................... 13 1.3 Poroelastic models ................................ 15 1.3.1 Classical Biot’s model .......................... 15 1.3.2 Non-dissipative poroelastic model (closed pores) ............ 18 1.3.3 Non-dissipative poroelastic model (open pores) ............. 21 1.3.4 Dissipative poroelastic model (open pore) . .............. 22 2 Finite element solution of acoustic propagation in rigid porous media 25 2.1 Introduction .................................... 26 2.2 Models for fluid-porous vibrations ........................ 27 2.3 Associated nonlinear eigenvalue problems .................... 30 2.4 Statement of the weak formulation ....................... 34 2.5 Finite element discretization ........................... 35 2.6 Matrix description ................................ 37 2.7 Numerical results ................................. 38 2.8 Conclusions .................................... 41 3 Finite element solution of new displacement/pressure poroelastic models in acoustics 43 3.1 Introduction .................................... 44 3.2 Statement of the problem ............................ 44 3.3 Weak formulation ................................. 47 3.4 Finite element discretization ........................... 48 3.5 Matricial description ............................... 51 i
ii Contents 3.6 Numerical solution of cell problems ....................... 53 3.7 Numerical results ................................. 55 3.8 Conclusions .................................... 57 II Perfectly Matched Layers 61 4 A non reflecting porous material: the Perfectly Matched Layers 63 4.1 Introduction .................................... 64 4.2 Wave equation ................................... 66 4.2.1 Time domain equations .......................... 66 4.2.2 Time-harmonic equations ........................ 67 4.3 Cartesian Perfectly Matched Layers ....................... 67 4.3.1 Time-domain equations .......................... 68 4.3.2 A physical interpretation ......................... 68 4.3.3 Time-harmonic equations ........................ 70 4.4 Plane wave analysis of the PMLs ........................ 73 5 An optimal PML in Cartesian coordinates 79 5.1 Introduction .................................... 80 5.2 The time-harmonic acoustic scattering problem ................ 80 5.3 Finite element discretization. ........................... 82 5.4 Determination of the absorbing function .................... 85 5.5 Comparison with classical absorbing functions ................. 89 5.6 Numerical tests .................................. 92 5.7 Conclusions .................................... 94 5.8 Computation of the element matrices ...................... 95 6 An exact bounded PML in radial coordinates 103 6.1 Introduction ....................................104 6.2 Scattering problem ................................105 6.3 Statement of the PML equation .........................106 6.4 PML fundamental solution ............................108 6.5 PML integral representation formula ......................115 6.6 Addition theorem .................................118 6.7 Existence and uniqueness of solutions for the PML equation .........120 6.8 Coupled fluid/PML problem ...........................122 6.9 Discretization and numerical results .......................125 6.10 Appendix .....................................128 6.10.1 Technical results .............................128 6.10.2 Some classical results about the Hankel functions ...........135
Contents iii III Computational applications on dissipative acoustics 139 7 Validation of acoustic dissipative models 141 7.1 Introduction ....................................142 7.2 Statement of the problem. Mathematical modeling ..............144 7.2.1 The Allard-Champoux model ......................144 7.2.2 The wall impedance model ........................145 7.2.3 Computing the wall impedance .....................146 7.3 Planar unbounded wall ..............................148 7.3.1 Plane waves with oblique incidence ...................148 7.3.2 Spherical waves ..............................152 7.4 Curved wall ....................................156 7.4.1 The Perfectly Matched Layer ......................157 7.4.2 Finite-element discretization .......................159 7.4.3 Verification of the numerical methods ..................161 7.4.4 Numerical validation of the wall impedance model for non-planar geometries ..................................164 7.5 Conclusions ....................................166 8 Numerical simulation of locally reacting panels 169 8.1 Introduction ....................................170 8.2 Modelling the panel ................................171 8.2.1 Wall-like impedance ...........................172 8.2.2 Porous veil and micro-perforated plates .................173 8.2.3 Thin porous layer .............................174 8.2.4 Multilayer panel with a rigid back ....................175 8.3 Variational formulation ..............................176 8.3.1 Wall-like impedance ...........................177 8.3.2 Porous veil .................................178 8.3.3 Thin porous layer .............................178 8.3.4 Multilayer panel with a rigid back ....................179 8.4 Finite element discretization ...........................179 8.5 Numerical validation ...............................180 8.6 Numerical results for an absorbing box for reducing noise in rooms ......185 8.7 Conclusions ....................................187 Further research 191 Acknowledgments 192 Resumo en galego 195 Bibliography 203
iv Contents
Preface Nowadays the numerical methods have a fundamental role as a tool to reduce the design time and developed costs of new products in fields such as Aerospace, Mechanic, Naval Engineering, etc. From this point of view, sometimes the quick evolution of computers is not enough in all the cases to solve the real-life problems of engineering efficiently and in a practical time. Hence the computational capability of actual computers has to be completed with efficient and renewed numerical methods. One of the problems which has become more relevant from a social point of view is the reduction of acoustic pollution produced by cars, planes, air-conditioned systems, etc., as it is reflected in national and European laws, which are more restrictive in the last years. In this context, it arises the necessity to solve more complex acoustic propagation problems which cannot be tackled with by numerical techniques based on classical methods. Obviously, prototype essays are fundamental to asses the feasibility of the proposed technologies, nevertheless the high cost of prototype production makes necessary that this kind of experiments has to be done in an advanced phase of design, with a product close to the final one. These two factors are the reasons why computational acoustics becomes a scientific field of great importance nowadays and numerical simulation is a relevant tool to do analysis of products and to study innovative systems with comfortable acoustic properties with competitive cost and saving developing time. The sophistication of the acoustic materials related with the real-life problems in the last decades have caused that the mathematical models used to solve acoustic propagation problems have been enriched from a mathematical point of view, and consequently they require the use of new advanced computational and numerical techniques. Among these new models, we focus in this thesis on those derived from the porous materials and, from a computational point of view, on the Perfectly Matched Layer (PML) technique, which allows solving numerically acoustic propagation problems in unbounded domains. In any case, the complexity of the acoustic models and the geometrical configuration of the problems require their resolution to be done by numerical methods as, for instance, the finite element methods. The study presented through this thesis is in the context of the frequency domain, i.e., under the assumption of time-harmonic dependency of the time variable of the acoustic fields. In fact, our attention is focused in acoustic propagation problems in the low range of frequencies, where the discretization by finite element methods is suitable and non excessively expensive from a numerical point of view. v
Chapter 1 Porous models Contents 1.1 Introduction .............................. 4 1.2 Rigid porous models .......................... 6 1.2.1 Darcy’s like model .......................... 7 1.2.2 Allard-Champoux model ....................... 13 1.3 Poroelastic models ........................... 15 1.3.1 Classical Biot’s model ........................ 15 1.3.2 Non-dissipative poroelastic model (closed pores) ......... 18 1.3.3 Non-dissipative poroelastic model (open pores) .......... 21 1.3.4 Dissipative poroelastic model (open pore) ............. 22 3
4Chapter 1. Porous models 1.1 Introduction The main goal of this chapter is to introduce various mathematical models which can be used to characterize the acoustic behavior of porous materials. By porous material we mean a material consisting of a solid matrix which is completely saturated by a fluid. The acoustical behavior of porous media depends not only on the fluid but also on the stiffness of its solid skeleton. Acoustic behavior of porous media is of utmost importance because they exhibit good properties as sound absorbers. Such kind of materials, as glasswools or ceramic foams (see Figure 1.1) are used in isolation systems in buildings, vehicles or airplanes. Figure 1.1: A particular porous material: ceramic foams. When considering macroscopic models for porous media, they can be classified depending on whether the solid part is rigid or elastic: A) Rigid porous models. If the solid matrix is rigid, which is the simplest case, the porous material can be considered as an equivalent fluid, with complex mass density and bulk modulus. These parameters can be obtained through semi-empirical or experimental laws. Delany and Bazley [56] presented a first model in 1970, which has been widely used to describe sound propagation in fibrous materials. This model, was subsequently improved by Morse and Ingard [88], Attenborough [13] and Allard and Champoux [5], among others. All models of this kind have a common characteristic: they are stated under the assumption of time-harmonic dependency. Hence, the coefficients that appear in the models depend on the frequency. Generally, the equations of these models can be derived from the classical compressible fluid equations with some slight modifications on the coefficients, for instance, adding a damping term or becoming the mass density or the bulk modulus complex valued.
1.1. Introduction 5 B) Poroelastic models. For the more realistic case when the elastic deformation of the skeleton is taken into account, the theoretical basis for the mechanical behavior was mainly established by Biot [38]. His theory describes the propagation of elastic waves in fluid-saturated porous media. Adaption of this theory to the acoustic context was done, for example, by Allard et al. [4] and Shiau [98] (see also the book by Allard [3] and the references therein). In spite of the fact that the classical Biot’s model was introduced in the time domain, it has been also adapted to the frequency domain by Broubard and Lafarge [76] among others, which allows writing a fluid-equivalent formulation for poroelastic materials. One of the main drawbacks when Biot’s model is analyzed is that the coefficients are not properly defined and, in general, their determination is not clear although several experimental procedures have been proposed, as can be seen in Biot and Willis [40]. This gives motivation to undertake a derivation of models, rigorously from a mathematical point of view, by using homogenization theory. It can be done for the both cases: rigid or elastic solid matrix. For rigid matrix, Darcy’s model is obtained. Ene and Sanchez- Palencia [59] were the first who gave a derivation of this model from the Stokes system, using a formal multi-scale method, while Tartar [97] made that derivation rigorously in the case of 2D periodic porous media. This methodology allows obtaining not only the homogenized model but also the mathematical expression of the coefficients appearing in it. Derivation of macroscopic models for poroelastic materials depends strongly on the connectivity of the fluid part. When the domain occupied by the fluid is connected, the material is named open pore material; otherwise, it is named closed pore material. Fundamental references are papers by Gilbert and Mikeli´c [65], and by Clopeau et al. [49], where the classical dissipative Biot’s model was derived by homogenization, using two-scale convergence methods. They also contain a number of references to papers on the dissipative Biot’s law. Moreover, the same procedure has been applied, for the first time, by Ferr´ın and Mikeli´c [62] to derive macroscopic models for non-dissipative poroelastic material with open or closed pore. The outline of this chapter is as follows. In Section 1.2 we introduce the two rigid models that will be used through this work: the Darcy’s like model (Subsection 1.2.1) and the Allard-Champoux model (Subsection 1.2.2). In the case of Darcy’s like model, we also give a sketch of its derivation by homogenization techniques. In Section 1.3, we first introduce the classical Biot’s model (Subsection 1.3.1) and then we recall three poroelastic models obtained by homogenization techniques through the Subsections 1.3.2-1.3.4: nondissipative and open pore, non-dissipative and close pore, and dissipative and open pore. In every model the auxiliary cell problems which define the coefficients of the models are also detailed.
6Chapter 1. Porous models 1.2 Rigid porous models In the last decades, simplified models, where absorptive materials are characterized by normal wave impedance, were widely used to study wave propagation in rigid lined ducted systems. More recently, when the solid skeleton is assumed to be rigid, the porous material has been considered as an equivalent fluid, with complex mass density and bulk modulus. These parameters can be obtained through empirical or experimental laws. A first model, introduced by Delany and Bazley (see [56]), was presented for the first time in 1970; it has been widely used to describe sound propagation in fibrous materials. Subsequently, this model was improved by Morse and Ingard [88], Johnson et al. [73], Attenborough [12], Allard et al. [6], Champoux and Stinson [45] or Allard and Champoux [5], among others. But models simulating a slow flow fluid through porous media can be also derived rigorously from a mathematical point of view by means of homogenization techniques. By doing so, if we consider a rigid porous medium, we obtain the Darcy’s model. To the best of our knowledge, Ene and Sanchez-Palencia [59] were the first who gave derivation of it from the Stokes system, using a formal multiscale method. This derivation was made rigorous in the case of 2D periodic rigid porous media by Tartar (see appendix in [97]) and subsequently generalized by Mikeli´c, among others (see [86] and references therein). As we have said in the introduction of this chapter, by means of this methodology we obtain not only the homogenized model but also the mathematical expression of the coefficients appearing in it. For instance, in the case of rigid porous media, the most important coefficient in Darcy’s law is permeability, which can be computed by solving a boundary-value problem in a unit cell of the periodic porous medium. For poroelastic media, generalized Biot models were also derived from the first principles by using homogenization techniques (see Gilbert and Mikeli´c [65], Clopeau et al. [49] or Ferr´ın and Mikeli´c [62]). Classically, from a macroscopic point of view, rigid porous media are characterized modifying the conservative mass and momentum laws that model compressible fluids (see [12] or [88] among others), in order to take into account the friction phenomena and the energy exchange between the walls of pores of the solid material and the enclosed fluid. From these modifications, a dissipative term arises in the momentum conservation law, which models the friction phenomena, while the conservative mass law is modified to take into account that the movement of the fluid is restricted to a part of the porous material. Another fundamental feature consists in the isothermal character of the movement. Whereas, in non-dissipative acoustics, the motion is assumed to be isentropic, when a rigid porous material is modelled, we consider that the temperature is constant, i.e, we assume that the movement in the porous material produce a change in the entropy but not in the temperature. Anyway, the pressure and velocity fields computed with these models must be understood in a macroscopic sense, i.e., both the displacement and the pressure fields are only averaged estimates in a control volume element of the real pressure and velocity inside the porous material. Among the different models existing in the bibliography for rigid porous media, we
1.2. Rigid porous models 7 focus on a Darcy’s like model (obtained by homogenization techniques) and on another one obtained through empirical laws by Allard and Champoux [5], for harmonic motions. The latter generalizes the former, in the sense that both are equivalent for low frequencies when the porosity is near to one. The first one is derived from Darcy’s law, by adding the inertial effect (see for instance [13]). In fact, if the porous material have a periodic structure, the Darcy’s model can be obtained by applying a two-scale homogenization technique to the Stokes equations (see [59] and [97]). In this case, the only parameter characterizing the porous media (named flow resistivity), can be related to the permeability tensor, which can be computed explicitly by solving a Stokes problem defined in the fluid part of the unit cell. A finite element solution of the acoustic propagation in rigid porous media simulated with these two models has been carried out by Berm´udez et al. [27]. 1.2.1 Darcy’s like model If the porous material is homogeneous from a macroscopic point of view (i.e, if the pores that compose the micro-structure of the material are distributed uniformly) then, according to Darcy’s law, it is necessary to introduce two new parameters: the flux resistivity tensor, σ, which gives information about the resistance to the fluid movement exerted by the rigid skeleton of the porous material, and the porosity coefficient φ, which is the ratio between the volume of the fluid part and the total volume of the porous material. If we suppose that the fluid filling the pores of the material is compressible, homogeneous and isothermal, then the equations governing the movement of a porous material with rigid solid part are ρ∂V ∂t +ρ(grad V)V+grad P+σV =0,(1.1) ∂ρ ∂t +V·grad ρ+ρ φdiv V=0,(1.2) P=ρRθ, (1.3) where σis the flow resistivity tensor ([σ]=kg/(m3s)), φis the porosity coefficient and θ, P,ρand Vare the temperature, the pressure, the mass density and the velocity averaged in the part occupied by the fluid in a control volume element, respectively. If the acceleration term is negligible in comparison with the dissipative term (which depends on the flow resistivity), the linear momentum conservation law (1.1) becomes the Darcy’s law (see [21]) V=−σ−1grad P. (1.4) Usually, some authors (see [91]) write the above model in terms of the permeability tensor K([K]=m 2). The permeability tensor can be written in terms of the flow resistivity and the porosity coefficient as follows: K=ηφσ−1,
8Chapter 1. Porous models where ηis the viscosity of the fluid filling the pores of the material. Since it is possible to check that Konly depends on the geometry of the pores, if we assume that the structure of the solid part of the porous media is uniformly periodic, the permeability tensor can be computed applying homogenization techniques. Now, if we write the pressure in terms of mass density, ρ, and entropy, s, since we are assuming that the flow of the fluid is isothermal, the value of the sound speed is given by c=∂P ∂ρ (ρ, s)=γ∂P ∂ρ (ρ, θ0)=γRθ0,(1.5) where θ0is the constant temperature of the porous material and γis the ratio of specific heats at constant pressure and constant volume. Similar to the classical arguments used in the linearization of the equations which govern a compressible fluid (see for instance Section 11.2 in [25]), we are going to write the linearized equations for this rigid porous model. More precisely, we linearize the mass and momentum conservation laws in the neighborhood of a state at rest with constant mass density ρ0an pressure P0. Since the movement is isothermal, P≈P0+∂P ∂ρ (ρ0,θ 0)(ρ−ρ0)+∂P ∂θ (ρ0,θ 0)(θ−θ0)=P0+c2 γ(ρ−ρ0).(1.6) If we neglect the term V·grad ρand integrate (1.2) in time, we have ρ=ρ0−ρ0 φdiv U, where Uis the displacement field. So, using (1.6) we can write the pressure in terms of the displacement, P=P0−ρ0c2 φγ div U, whereas, the displacement field satisfies the equation ρ0 ∂2U ∂t2−ρ0c2 φγ grad(div U)+σ∂U ∂t =0.(1.7) On one hand, if the porous material is isotropic from a macroscopic point of view, i.e, if the flow resistivity tensor σis proportional to the identity, σ=σI, then (1.7) can be rewritten as ρ0 ∂2U ∂t2−ρ0c2 φγ grad(div U)+σ∂U ∂t =0. On the other hand, if the porous material is fibrous, i.e., if the acoustic propagation in the parallel and orthogonal directions to the fibers have different characteristics, then the porous material would be modelled as an orthotropic medium. In this case, if we assume
1.2. Rigid porous models 9 that the principal directions of σmatch with the coordinates axes, then it is diagonal, i.e., σ=σ1e1⊗e1+σ2e2⊗e2+σ3e3⊗e3. Finally, in the frequency domain the governing equations of the movement in the Darcy’s like model are (ω2ρ0I+iωσ)u+ρ0c2 φγ grad(div u)=0,(1.8) p=P0−ρ0c2 φγ div u,(1.9) where uand pare the phasors of the displacement and the pressure fields in the frequency domain, respectively. Let us remark that we are using e−ωt for the time convention of the harmonic dependency. Computing the permeability tensor To finish this subsection, we show how to compute the permeability tensor for the rigid porous media by using homogenization techniques. In fact, we summarize the main results of the classic two-scale technique to determine the permeability tensor in the Darcy’s law, where the fluid is considered viscous and incompressible (see [86]). Assuming that the porous media have rigid solid part, we have seen that the equations governing the movement involve the permeability tensor, which only depends on the geometry of the pores of the periodic solid structure. If the material can be constructed with a periodic pattern from a unique cell with size , the permeability tensor can be computed by solving a Stokes problem in only one of these cells. We denote by Ω ⊂Rn(n= 2 or 3) the domain occupied by the porous material whereas Ωis the domain occupied by its fluid part. Since the porous material has a periodic structure, if we denote by Y=(0,1)nthe unit cell of reference and by YFand YS its fluid and solid part, respectively, both the fluid part Ωand the solid part of the porous medium Ω \Ωcan be rewritten as Ω= j∈T Y Fj,Ω\Ω= j∈T Y Sj, where Y Fjand Y Sjare respectively the fluid and the solid part of the cell Y jof size and homeomorphic to Y(see Figure 1.2). The index set Tis defined by T={k∈ Zn:(YS+k)=Y Sk⊂Ω}. For the sake of simplicity in the exposition, we avoid to treat mathematically the exterior boundaries and consider that Ω = (0,L)nwith periodic boundary conditions on ∂Ω. Since the geometrical structure of the porous medium is fixed, we are in a position to precise what is the macroscopic problem that governs the slow flow of a viscous incompressible flow through the porous medium. With this aim, we use the following steady Stokes
10 Chapter 1. Porous models Solid part Fluid part Fluid part Solid part Figure 1.2: Rigid porous structure and unit cell. system in Ω(the fluid part of Ω): −ηΔv(x)+grad p(x)=0in Ω,(1.10) div v(x)=0 inΩ ,(1.11) v=0on ∂Ω\∂Ω,(1.12) pand vare L−periodic, (1.13) where ηis the dynamical viscosity, vthe velocity and pthe pressure. If we define the functional space Wby W=u∈H1(Ω)n:u=0on ∂Ω\∂Ω,uis L−periodic}, then the weak formulation of the problem (1.10)-(1.13) is: Find v∈Wwith div v=0in Ωand p∈L2(Ω)such that ηΩ grad v·grad ψdVx−Ω pdiv ψdVx=0,∀ψ∈W.(1.14) On the one hand, the elemental theory of elliptic problems guarantees the existence and uniqueness of a velocity field solution of the above weak problem. On the other hand, the pressure field pis unique up to a constant (see [86]), which is usually fixed such that the following condition holds: Ω pdVx=0. Our aim is to take limits when →0 in the system of equations (1.10)-(1.13). As a first step, we study the a priori estimates for the velocity vand the pressure p. With
1.2. Rigid porous models 11 this purpose, we extend vby zero in the rigid part of the porous medium and define an extension of the pressure as ˜p:= ⎧ ⎪ ⎨ ⎪ ⎩ pin Ω, 1 Y FjY Fj pdVxin Y Sj,for each j∈T, where Y Fjand Y Sjare the fluid and solid part of the cell Y j, respectively. If we take into account these extensions for the pressure and the velocity, the following a priori estimates are satisfied (see [86]): vL2(Ω)n≤C2,(1.15) grad vL2(Ω)9≤C, (1.16) ˜p−1 |Ω|Ω ˜pdVxL2(Ω) ≤1 |YF|p−1 ΩΩ pdVxL2(Ω)≤C, (1.17) where Cis a positive constant not necessarily the same at each occurrence. Because of these a priori estimates, we can assume that the pressure and the velocity fields have the asymptotic expansions v(x,y)=2v0(x,y)+3v1(x,y)+4v2(x,y)+..., (1.18) p(x,y)=p0(x,y)+p1(x,y)+2p2(x,y)+..., (1.19) where y=x/. These expansions are written in two scales, i.e., taking into account the macroscopic and microscopic levels of the geometry. Moreover, since the geometry of the medium is periodic, it is natural to assume a periodic dependency in the scale of variable y. From the two scales of the expansions, the derivatives are transformed and, hence, the differential operators can be rewritten as grad =gradx+−1grady,(1.20) div = divx+−1divy,(1.21) Δ=Δ x+2−1divxgrady+−2Δy,(1.22) where the subindexes denote the spatial variable involved in the differentiation. If we substitute the expressions (1.18)-(1.19) in (1.10)-(1.13), the lower order terms satisfy the following equations: •term O(−1), gradyp0(x,y)=0in Ω ×YF,(1.23) •term O(1), −ηΔyv0(x,y)+gradyp1(x,y)+gradxp0(x,y)=0in Ω ×YF,
18 Chapter 1. Porous models open pores). Besides, according to the different treatment of the dissipative effects in the models, we can also derive two different kind of them named dissipative and non-dissipative models. Firstly we focus our attention in the description of the non-dissipative model with closed pores, following the ideas presented by Ferr´ın and Mikelic in [49] and [62]. 1.3.2 Non-dissipative poroelastic model (closed pores) If we again assume that the porous medium is homogeneous and isotropic from a macroscopic point of view, then the equations governing the pressure field, P, and the displacement field, U, in the porous medium are ρ∂2U ∂t2−div (A[E(U)]) + 1 ˆc(B−φI)grad (div (B−φI)U)=0,(1.35) P=−1 ˆcdiv (B−φI)U.(1.36) We recall that φis the porosity, ρFand ρSare, respectively, the density of the fluid and of the solid part, and the coefficient ˆc, the symmetric tensor B, and the linear operator A depend on the geometrical shape of the pores of the porous material, but not on the spatial variables. Derivation of the model by two-scale homogenization technique For the sake of completeness in the exposition, in this subsection we describe how to obtain the systems of equations (1.35)-(1.36) by using a two-scale homogenization technique. Moreover, we also state the cell problems which define the coefficients that arise in this model. In the case of periodic porous media with elastic solid part and closed pores, we apply the same homogenization technique reviewed in Section 1.2.1 for porous media with rigid solid part. With this technique we obtain a model with a unique unknown displacement field to describe the motion in the porous material (see [48], [49] or [50]). In what follows, we use the same technique to derive the equations (1.35)-(1.36). Firstly, we use the same notation used for the fluid and solid domains and the cells as the ones introduced in Subsection 1.2.1. We only recall that Ω Fand Ω Sare, respectively, the fluid and solid part of the porous domain. Under the small deformations assumption, the macroscopic problem is given by the following fluid-structure problem, where the displacements are described in terms of the linearized equations in the fluid and solid domains:
1.3. Poroelastic models 19 ρS ∂2w ∂t2−div(S(w)) = 0in Ω S×(0,T),(1.37) S(w):=C[E(w)] = 2μE(w)+λ(tr(E(w)))Iin Ω S×(0,T),(1.38) ρF ∂2u ∂t2+grad p=0in Ω F×(0,T),(1.39) p=−ρFc2div uin Ω F×(0,T),(1.40) u·ν=w·νon Γ×(0,T),(1.41) S(w)ν=−pνon Γ×(0,T),(1.42) w=∂w ∂t =0in Ω S×{0},(1.43) u=∂u ∂t =0in Ω F×{0},(1.44) p,w,and uare L−periodic. (1.45) Here, uand pare, respectively, the fluid displacement and pressure fields, wis the solid displacement field, S(w) is the linear approximation of the Piola-Kirchhoff stress tensor and E(w)=1 2(grad w+(grad w)t) is the infinitesimal strain tensor associated to the solid displacement field. Finally ρFand ρSare the mass densities in the fluid and in the solid part, respectively, λand μare the Lam´e coefficients of the solid part, and cis the sound speed in the fluid part of the porous material. Now we describe the construction of the homogenized problem using the two-scale method. We recall that xis the slow spatial variable whereas yis the fast one. Since we have two scales, the spatial derivatives are transformed by Eqs. (1.20)-(1.22), introduced in Subsection 1.2.1, in terms of operators gradx,grady, divx, and divy, where the subscripts denote the spatial variable involved in the differential operator. Following the ideas presented by Ferr´ın & Mikeli´c [62], as a first step, we extend the fluid and solid displacement field to a unique field ˜ udefined in Ω by the expression ˜ u(x,t)=u(x,t)ifx∈Ω F, w(x,t)ifx∈Ω S,(1.46) and, analogously, a pressure field ˜pdefined as ˜p(x,t)=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ −ρFc2div u(x,t)+ρFc2 |Ω|Ω F div u(x,t)dVxif x∈Ω F, ρFc2 |Ω|Ω F div u(x,t)dVxif x∈Ω S. (1.47) If we denote by u0(x,t) and ˜p0(x,y,t) the limit fields when tends to zero for the displacement ˜ uand the pressure ˜p, respectively, and we define P(x,t) such that ˜p0(x,y,t)=χYF(y)P(x,t)+B(t),
20 Chapter 1. Porous models where B(t) only depends on time t∈(0,T), then u0(x,t) and P(x,t) satisfy the limit problem ρ∂2u0 ∂t2−divxAE(u0)−divx(PB)+|YF|gradxP=0,(1.48) −|YF| ρFc2 ∂P ∂t =|YF|divx∂u0 ∂t −B:∂ ∂t Ex(u0)−∂P ∂t Y divyw0dVy,(1.49) where ρ=ρF|YF|+ρS|YS|. Vector field w0, and tensor Band linear operator Aare defined from problems stated in only one of the cells of the porous medium. More precisely, let w0 be the solution of the following problem in the solid part of the cell Y: −divyCEy(w0)=0in YS,(1.50) −CEy(w0)ν=νon ∂YS\∂Y, (1.51) Y w0dVy=0,(1.52) where νis the unit normal vector to ∂YS\∂Y exterior to YS. Moreover, if for each 1 ≤ i, j ≤n, we define the vector wij as the solution of the following problem: divyCei⊗ej+ej⊗ei 2+Eywij =0in YS,(1.53) Cei⊗ej+ej⊗ei 2+Ey(wij) ν=0on ∂YS\∂Y, (1.54) YS wij dVy=0,(1.55) then the symmetric tensor Bis given by B=YS CEyw0dVy,(1.56) whereas every component Aklij of the linear operator A((A[E])kl =AklijEij), is given by Aklij =YS Cei⊗ej+ej⊗ei 2+Eywij dVykl ,1≤i, j, k, l ≤n. (1.57) We can easily check that Aklij =Alkij =Alkji. Hence, if U=u0denotes the macroscopic displacement and we integrate (1.49) with respect to time, we obtain ρ∂2U ∂t2−div (A[E(U)]) −(B−φI)grad P=0,(1.58) φdiv U=B:E(U)+ˆcP, (1.59)
1.3. Poroelastic models 21 where the differential operators are assumed to be written with respect to the macroscopic variable xand ˆc=−φ ρFc2+YS divyw0dVy.(1.60) If we rewrite the terms arising in the equations, taking into account that B:E(U)= div BU, and we eliminate the pressure in (1.48), then we have ρ∂2U ∂t2−div (A[E(U)]) + 1 ˆc(B−φI)grad (div (B−φI)U)=0,(1.61) P=−1 ˆcdiv (B−φI)U.(1.62) 1.3.3 Non-dissipative poroelastic model (open pores) Now, we consider a porous medium with elastic solid part and open pores. We apply the same homogenization technique that we have used in the previous section to obtain other new model. If we again assume that the porous medium is homogeneous and isotropic, from a macroscopic point of view, then the equations satisfied by the pressure field Pand the displacement Uin the porous medium are (ρI−ρFA)∂2U ∂t2−div (A[E(U)]) −(A+B−φI)grad P=0,(1.63) ˆc∂2P ∂t2+1 ρF div (Agrad P)=−div (A+B−φI)∂2U ∂t2,(1.64) where we recall that φis the porosity coefficient, ρFand ρSare, respectively, the mass density of the fluid and solid part, ρ=ρF|YF|+ρS|YS|. As we have made in the previous subsection, we are going to summarize how to obtain this model by using the two-scale homogenization technique. If we state the same fluid-structure problem as in (1.37)-(1.45) and define the same extended pressure, ˜p, and displacement field ˜u, given by (1.47) and (1.46), respectively, then u0(x,t) and P(x,t) satisfy the limit problem ρ∂2u0 ∂t2−AgradxP+ρF ∂2u0 ∂t2= divxAE(u0)+ divx(PB)−|YF|gradxP, (1.65) −|YF| ρFc2 ∂P ∂t =divx(|YF|I−A)∂u0 ∂t −1 ρF At 0 gradxP(x,τ)dτ −B:∂ ∂t Ex(u0)−∂P ∂t Y divyw0dVy,(1.66) where we recall that ρ=ρF|YF|+ρS|YS|.
22 Chapter 1. Porous models Tensor Band linear operator Aare defined by (1.56) and (1.57), respectively, also using the same cell problems (1.50)-(1.52) and (1.53)-(1.55) stated in the previous subsection. We recall that the coefficient ˆc, the tensors Aand B, and the operator Adepend on the geometrical shape of the pores of the material and on the properties of the fluid and solid parts, but not on the spatial variables. In order to define tensor A, we need to introduce a new cell problem. For each 1 ≤j≤n, we define the scalar field ξjas the solution of the following problem in the fluid part of the unit cell: −Δyξj=0 inYF,(1.67) ∂ξj ∂ν=ej·νon ∂YF\∂Y, (1.68) ξjis 1 −periodic in Y, (1.69) YF ξjdVy=0,(1.70) then each component Aij of tensor Ais given by Aij =YFδij −∂ξj ∂yidVy,1≤i, j ≤n. If U=u0denotes the macroscopic displacement and we integrate equation (1.66) with respect to time, we obtain ρ∂2U ∂t2−Agrad P+ρF ∂2U ∂t2−div (A[E(U)]) −(B−φI)grad P=0,(1.71) div φ∂2U ∂t2−1 ρF Agrad P−A∂2U ∂t2=B:E∂2U ∂t2+ˆc∂2P ∂t2.(1.72) where the differential operators are assumed to be written with respect to the macroscopic variable xand ˆcis again given by (1.60). If we rewrite the terms which arise in the equations taking into account that B:E(U) = div BU, and eliminate the pressure in the first of these equations, we finally obtain the system of equations (1.63)-(1.64). In spite of the fact that this model only uses a unique displacement field U, which is an advantage with respect to the classical Biot’s model, the inclusion of the pressure gradient in equation (1.64) does not allow decoupling the equations as we have done in the model for closed pores. For the sake of completeness in the exposition we also summarize the model when a dissipative behavior is assumed in a porous medium with open pores. 1.3.4 Dissipative poroelastic model (open pore) The dissipative generalized Biot model has been derived, by using homogenization techniques, in Clopeau et al. [49] for the case where fluid viscosity ηis of order O(ε2), εbeing the size of the elementary cell. The equations are given by
1.3. Poroelastic models 23 ρI∂2U ∂t2−d dt t 0 A(t−τ)grad P(x,τ)+ρF ∂2U ∂τ2(x,τ)dτ −div (A[E(U)]) −(B−φI)grad P=0, and div φ∂U ∂t −t 0 A(t−τ)1 ρF grad P(x,τ)+∂2U ∂τ2(x,τ)dτ= div B∂U ∂t +ˆc∂P ∂t , where we recall that ρFis the fluid density, ρSis the density of the solid skeleton, ρ= φρF+(1−φ)ρSand Aklij := YS Cei⊗ej+ej⊗ei 2+Eywij dVykl , B:= YS CEy(w0)dVy, ˆc:= YS divyw0dVy, Aij(t):=YF wij y,ρF ηtdVy, which implies, in particular, Aij(0) = !YFwij(y,0) dVy=φδij. Let us remark that tensor B, linear operator A, and coefficient ˆccoincide with the expressions given for the non-dissipative models, since the vector fields w0and wij are solution of the same boundary-value problems (1.50)-(1.52) and (1.53)-(1.55), respectively.
24 Chapter 1. Porous models
Chapter 2 Finite element solution of acoustic propagation in rigid porous media Contents 2.1 Introduction .............................. 26 2.2 Models for fluid-porous vibrations ................. 27 2.3 Associated nonlinear eigenvalue problems ............. 30 2.4 Statement of the weak formulation ................. 34 2.5 Finite element discretization ..................... 35 2.6 Matrix description ........................... 37 2.7 Numerical results ........................... 38 2.8 Conclusions ............................... 41 25
26 FEM for rigid porous media models in acoustics 2.1 Introduction In this chapter, only the case of rigid frame porous material will be considered. Two models will be taken into account: the above mentioned sort of Darcy’s model and the Allard-Champoux model (see [5]), both presented in Section 1.2. As we have shown in the previous chapter, the main difference between them lies in the frequency dependence of the coefficients, namely, the mass density and bulk modulus. In fact, through this chapter we will focus our attention on the numerical computation of the resonance frequencies and the frequency response of some acoustic problems which involve these two models, which will be always stated in bounded domains. With this purpose we have implemented a finite element method. Because its easy implementation and its effectiveness in handling complex geometries, the finite element method has become popular to solve such problems. Some examples of the finite element method applied to sound propagation in poroelastic media are in the papers by Easwaran et al [58], Panneton and Atalla [92], G¨oransson [68] or Atalla et al [11]. All of them, take Biot’s general theory as the starting point. Other kind of problems concerning porous materials, related to vibration modes, were solved by Berm´udez et al [36]. More precisely, a finite element method introduced by Raviart and Thomas [96] will be used to solve numerically the two models which are formulated in displacements. In fact, if we consider a tetrahedron partition of the computational domain, the degrees of freedom of the Raviart-Thomas elements are the normal displacement in each face of the tetrahedra. So the divergence of the displacement field is conserved in the continuous and discrete problem and, as it has been proven in Berm´udez et al [26] (see also [35] and references therein), these finite elements do not produce spurious modes. Numerical experiments using both models will be presented for different three-dimensional examples. More precisely, we solve the source problem associated with an external harmonic excitation which allows us to know the response of the porous material. We also solve the nonlinear spectral problem associated with it. The outline of this chapter is as follows. In Section 2.2 we present the two models associated with the problem consisting of a finite two-layer system with rigid porous materials. They will be stated in the frequency domain leading to the response problem and to a nonlinear eigenvalue problem. In Section 2.3 the free vibration problem associated with this nonlinear eigenvalue problem is analyzed in order to obtain a deeper insight of the overdamped vibration frequencies. In Section 2.4 weak formulations for both problems are presented and an analysis of overdamped vibration frequencies is made. In Section 2.5 the finite element method is introduced, whereas in Section 2.6 the corresponding matrix description is shown. Finally, in Section 2.7, numerical results for some 3D examples are given for both the response and the spectral problems.
2.2. Models for fluid-porous vibrations 27 2.2 Models for fluid-porous vibrations Let us consider a coupled system consisting of an acoustic fluid (i.e. compressible barotropic inviscid) and a porous medium contained in a three-dimensional cavity. Let ΩF and ΩAbe the domains occupied by the fluid and the porous medium, respectively (see Figure 2.1). The boundary of ΩF∪ΩA, denoted by Γ, is the union of two parts, ΓDand ΓE.Γ Ddenotes the rigid walls of the cavity. Let νthe outward unit normal vector to Γ. We assume the interface between the fluid and the porous media, denoted by ΓI, is the union of surfaces, Γ0,Γ1,...,ΓJ. Let nbe the unit normal vector to this interface pointing outwards ΩA. Figure 2.1 shows a vertical cut of the domain for a better understanding of the notation. ΩA ΩF n ν ΓD ΩF Γ0Γ1 ΓE Figure 2.1: 3D domain and vertical cut. For studying the response of the coupled system (fluid-porous medium), subject to harmonic forces acting on ΓE, we consider two different models for the vibrations in the porous medium: Darcy’s like model and Allard-Champoux model. Both models assume the skeleton of the porous media is rigid. Firstly, the governing equations for free small amplitude motions of an acoustic fluid filling ΩFare given in terms of displacement and pressure fields by ρF ∂2UF ∂t2+grad PF=0in ΩF,(2.1) PF=−ρFc2div UFin ΩF,(2.2) where PFis the pressure, UFthe displacement field, ρFthe density and cthe acoustic speed in the fluid. Secondly, let us recall the rigid porous models introduced in Section 1.2 of the previous chapter. The Darcy’s like model only has slight differences with respect to the above fluid model. One of them consists of an additional damping term, named Darcy’s term (see [3]). Moreover, the interstitial fluid flow is supposed to be isothermal, a standard assumption
34 FEM for rigid porous media models in acoustics λ=−iω βFβF(ρF+σ/iω) cosh(aAβA)sin(aFβF)+βAρFsin(aAβA)cos(aFβF)=0 j =0,k=0 j =1,k=0 j =1,k=1 j =1,k=2 j =2,k=0 j=3,k=2 j=3,k=3 Figure 2.4: Curves (2.52) and (2.54) for real values of βFand λ=−iω; case aA> πρFc σ"1 3−φγ . 2.4 Statement of the weak formulation For the sake of simplicity, we restrict our attention to the case where the porous medium is isotropic, i.e., μ(ω)=μ(ω)I. Let us define the set Vof kinematically admissible virtual displacements, V={(vF,vA)∈H:vF·n=vA·non ΓI}, where H={(vF,vA)∈H(div,ΩF)×H(div,ΩA):vF·ν= 0 on Γ ∩∂ΩF, vA·ν= 0 on Γ ∩∂ΩA}, and H(div,Ω) = v∈(L2(Ω))3: div v∈L2(Ω)%, where L2(Ω) denotes the space of square integrable functions. To get a weak formulation of the eigenvalue problem (2.32)-(2.39), equation (2.32) is multiplied by the conjugate of a virtual fluid displacement ¯ vFsatisfying the Dirichlet condition (2.38) and then integrated in ΩF. By using a Green’s formula and equation (2.34), we obtain ΩF ρFc2div uFdiv ¯ vF−ΓI pF¯ vF·n=ω2ΩF ρFuF·¯ vF.
2.5. Finite element discretization 35 In an analogous way, equations (2.33), (2.35) and (2.39) yield ΩA μ(ω) div uAdiv ¯ vA+ΓI pA¯ vA·n=ω2ΩA ρ(ω)uA·¯ vA. Now, by adding both equations and using the kinetic constraint (2.37), we can write the following pure displacement eigenvalue problem: Find a complex angular frequency ωand a pair of displacements (uF,uA)∈V,with uF and uAnot both identically zero, satisfying ΩF ρFc2div uFdiv ¯ vF+ΩA μ(ω) div uAdiv ¯ vA= ω2ΩF ρFuF·¯ vF+ΩA ρ(ω)uA·¯ vA,(2.55) for all (vF,vA)∈V. As it is typical in displacement formulations (see [26]), ω= 0 is an eigenfrequency of this problem in both the Darcy’s like model and the Allard-Champoux model, with an infinite-dimensional eigenspace given by Z={(uF,uA)∈V: div uF= 0 in ΩF,div uA= 0 in ΩA}. This eigenspace consists of pure rotational fluid motions inducing neither variations of pressure in the fluid nor in the porous medium. They are mathematical solutions of the eigenvalue problem with no physical entity because they do not correspond to vibration modes of the coupled system. They arise because no irrotational constraint is imposed to the fluid and porous displacements (see [26]). 2.5 Finite element discretization Fluid and porous displacements belong to the same class of functional spaces, H(div,ΩF) and H(div,ΩA), respectively; hence the same type of finite elements should be used for each of them to discretize the variational problem (2.55). Let Thbe a regular tetrahedral partition of ΩF∪ΩAsuch that every tetrahedra is completely contained either in ΩFor in ΩA. We also assume that the faces of tetrahedra lying on ΓD∪ΓEare completely contained either in ΓDor in ΓE. To approximate the fluid and porous displacements, the lowest order Raviart-Thomas elements (see [96]) are used to avoid spurious modes typical of displacement formulations (see [74]). They consist of vector valued functions which, when restricted to each tetrahedron, are incomplete linear polynomials of the form uh(x1,x 2,x 3)=(a+dx1,b+dx2,c+dx3), a,b,c,d∈C. These vector fields have constant normal components on each of the four faces of a tetrahedron (Figure 2.5) which define a unique polynomial function of this type. Moreover, the
36 FEM for rigid porous media models in acoustics global discrete displacement field uhis allowed to have discontinuous tangential components on the faces of the tetrahedra of the partition Th. Instead, its constant normal components must be continuous through these faces (these constant values being the degrees of freedom defining uh). Because of this, div uhis globally well defined in the domain, ΩF∪ΩA. Figure 2.5: Raviart-Thomas finite element. Then, for fluid displacements we use the Raviart-Thomas space (see [96]) Rh(ΩF):={u∈H(div,ΩF):u|T∈R 0(T),∀T∈T h,T⊂ΩF}, and an analogous space for porous medium displacements: Rh(ΩA):={u∈H(div,ΩA):u|T∈R 0(T),∀T∈T h,T⊂ΩA}, where R0(T):=u∈P 1(T)3:u(x1,x 2,x 3)=(a+dx1,b+dx2,c+dx3), a,b,c,d∈C%. Then, the discrete analogue of Vis Vh:= {(uF,uA)∈Rh(ΩF)×Rh(ΩA):uF·n=uA·nfor each face on ΓI,uF·ν= 0 on Γ ∩∂ΩF,uA·ν= 0 on Γ ∩∂ΩA}. With this finite element space we define an approximate problem to (2.55): Find a complex number ωhand a pair of displacements (uh F,uh A)∈Vhnot both identically zero, such that ΩF ρFc2div uh Fdiv ¯ vh F+ΩA μ(ωh) div uh Adiv ¯ vh A= ω2 hΩF ρFuh F·¯ vh F+ΩA ρ(ωh)uh A·¯ vh A,(2.56) for all (vh F,vh A)∈Vh.
2.6. Matrix description 37 2.6 Matrix description In the previous section, a discrete formulation of our eigenvalue problem has been stated. Now a matrix description is given and it is shown that it is a well posed symmetric nonlinear generalized eigenvalue problem involving sparse matrices. Let uh Fand vh Fdenote the column vectors of components of uh Fand vh F, respectively, in the standard finite element basis associated with Rh(ΩF). Similarly, let uh Aand vh Adenote the column vectors of components of uh Aand vh A, respectively, in the standard finite element basis associated with Rh(ΩA). Then the problem (2.56) can be written in matrix form as RF0 0RA(ωh)uh F uh A=ω2 hMF0 0MA(ωh)uh F uh A,(2.57) where vh F ∗RFuh F=ΩF ρFc2div uh Fdiv ¯ vh F, vh A ∗RA(ωh)uh A=ΩA μ(ωh) div uh Adiv ¯ vh A, vh F ∗MFuh F=ΩF ρFuh F·¯ vh F, vh A ∗MA(ωh)uh A=ΩA ρ(ωh)uh A·¯ vh A. RFand MFare the standard stiffness and mass matrices of the fluid, respectively, while RA(ωh) and MA(ωh) are the corresponding ones for the porous medium. Notice that every matrix is highly sparse because only a maximum of seven entries per row can be different from zero (this corresponds to the number of faces of two adjacent tetrahedra). Matrices RFand RA(ωh) in the eigenvalue problem (2.57) are singular; however, by performing a translation in the eigenvalues, it can be written in an equivalent more convenient way: RF+MF0 0RA(ωh)+MA(ωh)uh F uh A=(ω2 h+1)MF0 0MA(ωh)uh F uh A. Now, matrix RF+MFis clearly positive definite, hence non-singular and symmetric. However, the matrix RA(ωh)+MA(ωh) is singular if there exists ωhsuch that ρ(ωh) is null. When we use the Darcy’s like model, the dynamic density is null only if ωh=iσ ρF(which is an eigenvalue associated to null divergence displacements in the porous medium) while the bulk modulus is positive. In the case of Allard-Champoux model, there exists a frequency ωhsuch that μ(ωh) is null but this is not true for ρ(ωh). Thus, except for this special case, the matrix on the left hand side is non-singular and, consequently, it can be used to build a well posed generalized eigenvalue problem to help us to solve the non-linear eigenvalue problem (2.57). With this aim, we can define a function
38 FEM for rigid porous media models in acoustics S:C→Csuch that S(ωh)=λh, where λhis the least modulus eigenvalue of the following problem, RF+MF0 0RA(ωh)+MA(ωh)uh F uh A=λhMF0 0MA(ωh)uh F uh A. The function Sis well defined because the generalized eigenvalue problem is well posed. Furthermore, both matrices of this problem are symmetric and highly sparse and, hence, convenient for computational purposes. Finally, calculation of the eigenvalues of the problem (2.57) is equivalent to find the roots ωhof the nonlinear equation, S(ωh)−(ω2 h+1)=0. A similar problem arising from finite element analysis of dissipative acoustic models can be found in [36]. 2.7 Numerical results In this section we present some numerical results obtained with a computer code implementing the numerical method given in this chapter. This code allows us to compute the response diagram of enclosures as those in Figures 2.1 and 2.2, consisting of several layers of fluids and porous media, and also to solve the nonlinear eigenvalue problem (2.55) by means of a secant method combined with an inverse power method. This type of method has been already used in other works devoted to solve nonlinear eigenvalues problems (see [36]). In order to validate our method, we have considered the following data: fluid is air with density ρF=1.225 kg/m3,c= 343 m/s whereas properties of the porous material are summarized in σ=10 4kg/(m3s), φ=0.71, γ=1.4, Npr =0.702 and P0= 101320 Pa. Concerning dimensions of enclosures shown in Figures 2.1 and 2.2, they are as follows: length and width are 4 m whereas height is 2 m for the first layer of free fluid, 0.05 m for the second layer of porous material and 0.1 m for the third layer of free fluid in the case of the enclosure in Figure 2.1. These two enclosures have been decomposed in tetrahedra (see Figure 2.6). Depending on which one is considered and on the degree of mesh refinement (parameter nrefers to the number of divisions introduced for each layer of the enclosure in Figure 2.6), the meshes are denoted as it is shown in Table 2.1. Firstly we consider the enclosure in Figure 2.2. In Table 2.2 we show the first complex eigenfrequencies (in Hz) for three different meshes: mesh 2, mesh 3 and mesh 4. One can see that they have a small imaginary part and a real part “close” to the response peaks. We also include the extrapolated complex eigenfrequencies computed by the least square method, and the exact ones corresponding to both Darcy’s like model and Allard-Champoux model, obtained by solving their respective nonlinear system of equations, for any pair of integers
2.7. Numerical results 39 Figure 2.6: Mesh 2 corresponding to the enclosure shown in Figure 2.2 Mesh 0 Mesh 1 Mesh 2 Mesh 3 Mesh 4 Sample Figure 2.1 Figure 2.1 Figure 2.2 Figure 2.2 Figure 2.2 n2 4 6 8 10 d.o.f 1072 8128 21600 50688 98400 Table 2.1: Name and degrees of freedom for the different meshes. jand k, namely, β2 F+ω2 c2=π2 b2j2+π2 d2k2, β2 A+ω2ρ(ω) K(ω)=π2 b2j2+π2 d2k2, βFρ(ω) cosh(aAβA) sinh(aFβF)=−βAρFsinh(aAβA) cosh(aFβF). An excellent agreement can be observed between exact and computed values, even for the coarser mesh. This shows the effectiveness of the method. From this Table we have calculated the order of convergence of the method and found that it is approximately O(h2), hbeing a parameter associated with the mesh size, which is optimal for the lowest order Raviart-Thomas finite elements we have used. On the other hand, according to the analysis in Section 2.3, there exist overdamped modes. In spite of the fact that these overdamped modes are not the magnitudes of interest, from the numerical point of view it is important to know if they are well approximated by the finite element method. Otherwise they could be a source of spectral pollution. Table 2.3 includes the computed and exact purely imaginary eigenfrequencies of higher modulus for the same three meshes described above.
40 FEM for rigid porous media models in acoustics Mode Mesh 2 Mesh 3 Mesh 4 Extrapolated Exact ωF 100 265.791+0.382 i 265.891+0.287 i 265.961+0.242 i 266.708+0.162 i 266.026+0.164 i ωF 110 375.722+0.786 i 375.894+0.593 i 375.977+0.502 i 376.144+0.338 i 376.132+0.342 i ωF 001 523.913+1.943 i 524.685+1.284 i 525.038+0.965 i 525.644+0.332 i 525.656+0.389 i ωF 200 531.251+1.704 i 531.418+1.284 i 531.506+1.087 i 531.733+0.772 i 531.680+0.737 i ωF 020 528.894+1.684 i 530.052+1.273 i 530.622+1.081 i 531.807+0.730 i 531.680+0.737 i Table 2.2: Rigid porous medium and air. Darcy’s like model. Mode Mesh 2 Mesh 3 Mesh 4 Extrapolated Exact ωF 00 8163.265 i 8163.265 i 8163.265 i 8163.265 i 8163.265 i ωF 10 8156.630 i 8157.003 i 8157.277 i 8163.163 i 8158.435 i ωF 11 8149.990 i 8150.735 i 8151.285 i 8163.382 i 8153.605 i Table 2.3: Rigid porous medium and air. Darcy’s like model. Overdamped modes. When the Allard-Champoux model is used we obtain similar results to those previously obtained for the Darcy’s like model. They are shown in Tables 2.4 and 2.5. Mode Mesh 2 Mesh 3 Mesh 4 Extrapolated Exact ωF 100 264.498+0.768 i 264.545+0.582 i 264.570+0.495 i 264.641+0.336 i 264.622+0.340 i ωF 110 373.830+1.586 i 373.915+1.212 i 373.968+1.033 i 374.195+0.692 i 374.082+0.712 i ωF 001 519.622+4.457 i 519.968+3.274 i 520.166+2.687 i 520.829+1.440 i 520.584+1.606 i ωF 200 528.662+3.436 i 528.559+2.652 i 528.550+2.266 i 528.535+1.461 i 528.598+1.561 i ωF 020 526.277+3.395 i 527.194+2.630 i 527.670+2.254 i 528.816+1.475 i 528.598+1.561 i Table 2.4: Rigid porous medium and air. Allard-Champoux model. Nevertheless, when calculating overdamped modes with the Allard-Champoux model some difficulties appear due to the highly oscillating eigenfunction associated with it. This means that a very refined mesh must be used in order to get a suitable approximation. In fact, only with mesh 4 a good accuracy has been achieved. Figure 2.7 shows the oscillation in the eigenvector in a plane near the interface between the fluid and the porous material. Finally, we consider the enclosure shown in Figure 2.1. The response curves are drawn in Figure 2.8 when the model is solved with mesh 0 and mesh 1. In these curves log10 ||p||L2 is plotted for frequencies ranging from 50 to 1000 Hz. Several response peaks can be observed in these curves depending on the refinement of the mesh. We notice that the finer the mesh, the smaller the number of peaks in the response diagram. Moreover, Table 2.6 shows the computed (complex) resonance frequencies for the damped coupled system shown in Figure 2.1 and the (real) eigenfrequencies of a similar undamped enclosure where the porous material has been replaced with air. As one can see, all are quite similar.
2.8. Conclusions 41 Mode Mesh 2 Mesh 3 Mesh 4 Exact ωF 100 4983.698 i 4439.175 i 4439.603 i 4438.380 i ωF 110 4979.907 i 4436.174 i 4436.828 i ωF 001 4976.004 i 4435.718 i 4435.633 i 4435.279 i Table 2.5: Rigid porous medium and air. Allard - Champoux model. Overdamped modes. Figure 2.7: Eigenvector for an overdamped mode with Allard-Champoux model 2.8 Conclusions In this chapter a three-dimensional finite element method has been implemented to solve the system of equations modelling the macroscopic behavior of a porous material with rigid solid frame. It allows us to compute both the response to a harmonic excitation and the free vibrations of a three-dimensional multilayer system consisting of different layers composed of free fluids and rigid porous media. The finite element used is the lowest order face element introduced by Raviart and Thomas, with the advantage of eliminating the spurious modes. For rigid porous media we have considered two models: a Darcy’s like model and the Allard-Champoux model. These two models are equivalent when frequency is much lower than flow resistivity. When solving the problem of free vibrations, the computer program predicts very well the exact complex eigenfrequencies for the Darcy’s like model in the case of a test example. This is true even for the overdamped modes. On the other hand, when using the Allard-Champoux model, the eigenfrequencies with non-null real part are well approximated whereas calculation of overdamped modes is much more complicated due to the highly os-
42 FEM for rigid porous media models in acoustics Figure 2.8: Response curve with Mesh 0 (left) and Mesh 1 (right) Undamped peaks Complex eigenfrequencies Damped peaks Damped peaks Mesh 0 Mesh 1 269.392 266.189+1.495 i 250.5 250.5 588.840 594.737+8.535 i 572.5 581.0 751.297 751.811+1.506 i 720.5 748.0 808.175 797.242+1.750 i 785.5 - 924.519 936.395+5.309 i 921.0 - Table 2.6: Resonance vibration frequencies and complex eigenfrequencies cillating eigenfunction associated with them, as observed in Figure 2.7. This forces us to use very fine meshes.
Chapter 3 Finite element solution of new displacement/pressure poroelastic models in acoustics Contents 3.1 Introduction .............................. 44 3.2 Statement of the problem ...................... 44 3.3 Weak formulation ........................... 47 3.4 Finite element discretization ..................... 48 3.5 Matricial description ......................... 51 3.6 Numerical solution of cell problems ................ 53 3.7 Numerical results ........................... 55 3.8 Conclusions ............................... 57 43
50 FEM for poroelastic models in acoustics This bubble function is a polynomial of degree four, null on the surface of tetrahedron T and taking value one at barycenter of T. The approximating space associated with the MINI element consists of continuous vector valued functions whose components, restricted to each tetrahedron, are sum of a bubble function and a polynomial of degree one, i.e., uh i(x1,x 2,x 3)|T=ax1+bx2+cx3+d+eα(x1,x 2,x 3), a,b,c,d,e∈C. The degrees of freedom for functions in this space are the values of the vector field at vertices and barycenters of tetrahedra (see Figure 3.3). Value at vertex Value at barycente r Figure 3.3: MINI finite element. Then, for porous displacements, we use the MINI space Mh(ΩA):=u∈H1(ΩA)3:u|T∈(P1(T)⊕Pb(T))3,∀T∈T h,T⊂ΩA%, where Pb(T)={aα :a∈C}. To approximate porous medium pressure, continuous piecewise linear finite elements are used. They consist of scalar valued functions which, when restricted to each tetrahedron, are polynomials of the form ph(x1,x 2,x 3)|T=ax1+bx2+cx3+d, a, b, c, d ∈C. Thus, the porous medium pressure is approximated in the finite-dimensional space, Lh(ΩA):=p∈H1(ΩA):p|T∈P 1(T),∀T∈T h,T⊂ΩA%. We recall that the degrees of freedom defining phare its values at vertices of tetrahedra. Finally, in order to approximate the interface pressure we use piecewise constant functions on the triangles of the mesh lying on the interface ΓI. In other words, for interface pressure we use the space Ch(ΓI):=p∈L2(ΓI):p|∂T ∈P 0(∂T),∀T∈T h,∂T∩ΓI=∅%.
3.5. Matricial description 51 The degrees of freedom of this finite element space are the (constant) values on triangles in ΓI. Consequently, the discrete analogue to Vis Vh=Rh(ΩF)×Mh(ΩA)×Lh(ΩA)×Ch(ΓI) while the corresponding to V0is V0h={(vF,q F,vA,q A)∈Rh(ΩF)×Mh(ΩA)×Lh(ΩA)×Ch(ΓI): vF·ν= 0 on (ΓD∪ΓE)∩∂ΩF,vA=0on ΓD∩∂ΩA}. With these finite element spaces we can define the approximate problem to (3.23)-(3.24) by For an angular frequency ωfixed, find (uh F,p h F,uh A,p h A)∈Vhsatisfying uh F·ν=0on faces in ΓW∩∂ΩF, uh A=0at vertices in ΓW∩∂ΩA, uh F·ν=fon faces in ΓE, and furthermore, ΩF ρFc2div uh Fdiv ¯ vh F−ω2ΩF ρFuh F·¯ vh F−ω2ΩA (ρI−ρFA)uh A·¯ vh A+ ΩA AH[D(uh A)] : D(¯ vh A)+ΩA div (A+BH−φI)t¯ vh Aph A+ΩA ˆcph A¯qh A+ ΩA 1 ρFω2Agrad ph A·grad ¯qh A+ΩA div (A+BH−φI)uh A¯qh A= ΓI ph F(¯ vh F·n−¯ vh A·n) (3.25) and ΓI ¯qh F(uh A·n−uh F·n)=0,(3.26) for all (vh F,qh F,vh F,qh A)∈V0h. 3.5 Matricial description In the previous Section, a discrete formulation of our source problem has been established. Now a matrix description is given and, assuming that it is well posed, we show that it is equivalent to a reduced linear system whose unknowns are the degrees of freedom of the interface pressure.
52 FEM for poroelastic models in acoustics Let Uh Fand Vh Fdenote the column vectors of nodal components of fluid displacement fields uh Fand vh F, in the standard finite element basis associated with Rh(ΩF), excluding those corresponding to faces on (ΓW∪ΓE)∩∂ΩF. Similarly, let (Uh A,Ph A) and (Vh A,Q h A) denote the column vectors of nodal components of the pair fields (uh A,p h A) and (vh A,qh A), in the standard finite element basis associated with Mh(ΩA)×Lh(ΩA),excluding those corresponding to vertices in ΓW∩∂ΩAfor Uh Aand Vh A. Lastly, let us call Ph Fand Qh Fthe vectors of nodal components of the interface pressure fields ph Fand qh F, in the space of finite elements Ch(ΓI). Then the discretization problem can be written in matrix form as ⎛ ⎜ ⎜ ⎝ RF−ω2MF00−DF 0Ru,A−ω2Mu,AC∗DA 0Cω −2Rp,A+Mp,A0 −D∗ FD∗ A00 ⎞ ⎟ ⎟ ⎠⎛ ⎜ ⎜ ⎝ Uh F Uh A Ph A Ph F ⎞ ⎟ ⎟ ⎠=⎛ ⎜ ⎜ ⎝ Bh 0 0 0 ⎞ ⎟ ⎟ ⎠(3.27) where (Vh F)∗RFUh F=ΩF ρFc2div uh Fdiv ¯ vh F, (Vh A)∗Ru,AUh A=ΩA A[E(uh A)] : E(¯ vh A), (Qh A)∗Rp,APh A=ΩA 1 ρF Agrad ph A·grad ¯qh A, (Vh F)∗MFUh F=ΩF ρFuh F·¯ vh F, (Vh A)∗Mu,AUh A=ΩA (ρI−ρFA)uh A·¯ vh A, (Qh A)∗Mp,APh A=ΩA ˆcph A¯qh A, (Qh A)∗CUh A=ΩA div (A+B−φI)uh A¯qh A, (Vh F)∗DFPh F=ΓI ph F¯ vh F·n, (Vh F)∗DAPh F=ΓI ph F¯ vh A·n, and the right hand side Bhcomes from boundary data f.RFand MFare the standard stiffness and mass matrices of the fluid, while Ru,A,Mu,Aand Rp,A,Mp,Aare the corresponding ones for the porous medium in the case of displacement and pressure, respectively. Notice that every matrix depending on fluid fields is highly sparse because only a maximum of seven entries per row can be different from zero (this corresponds to the number of faces of two adjacent tetrahedra).
3.6. Numerical solution of cell problems 53 Now, we are going to choose ωsuch that it is not an eigenvalue of the discrete problem. Specifically, we assume that the entire matrix of linear system (3.27) and matrices RF− ω2MFand KA=Ru,A−ω2Mu,AC∗ Cω −2Rp,A+Mp,A(3.28) are non-singular. On the other hand, we notice that matrix ω−2Rp,A+Mp,Ais clearly positive definite and hence non-singular. In order to improve the resolution of linear system (3.27), we are going to take into account the non-singularity of diagonal matrices. Firstly, we can rewrite this system as (RF−ω2MF)Uh F−DFPh F=Bh,(3.29) (Ru,A−ω2Mu,A)Uh A+C∗Ph A+DAPh F=0,(3.30) CUh A+(ω−2Rp,A+Mp,A)Ph A=0,(3.31) −D∗ FUh F+D∗ AUh A=0.(3.32) Since matrix RF−ω2MFis non-singular, because of the choice of frequency ω, and ω−2Rp,A+ Mp,Ais positive definite, we can obtain from (3.29) and (3.31) Uh F=(RF−ω2MF)−1(Bh+DFPh F),(3.33) Ph A=−(ω−2Rp,A+Mp,A)−1CUh A.(3.34) If we take into account that KAis non-singular and ω−2Rp,A+Mp,Ais positive definite, we can conclude that matrix Ru,A−ω2Mu,A−C∗(ω−2Rp,A+Mp,A)−1Cis also non-singular. Then, it results from (3.30) and (3.34) that Uh A=−(Ru,A−ω2Mu,A−C∗(ω−2Rp,A+Mp,A)−1C)−1DAPh F.(3.35) Collecting equations (3.33) and (3.35), we can write a simple linear system from (3.32) whose unique unknown is the nodal interface pressure vector Ph F, namely, D∗ A(Ru,A−ω2Mu,A−C∗(ω−2Rp,A+Mp,A)−1C)−1DA +D∗ F(RF−ω2MF)−1DF%Ph F=D∗ F(RF−ω2MF)−1Bh, where the involved matrix is non-singular because this system is equivalent to the full system (3.27) which has unique solution. 3.6 Numerical solution of cell problems In this section we solve the cell problems allowing us to obtain values of the macroscopic coefficients for the poroelastic model considered in this chapter. More precisely, we are going to solve problems (1.50)-(1.52), (1.53)-(1.55) and (1.67)-(1.70) presented in Chapter 1 for the fluid and solid cells shown in Figure 3.4. This is a first step to determine the
54 FEM for poroelastic models in acoustics components of operator A, matrices Aand B, and scalar ˆcappearing in the poroelastic model (3.3)-(3.4). We use a finite element method to solve these boundary-value problems. More precisely, we employ continuous piecewise linear finite elements on tetrahedral meshes to approximate wij,w0and ξi. These meshes can be seen in Figure 3.4. O X Y Z O X Y Z Figure 3.4: Meshes of the fluid and solid part of the unit cell (YFand YS, respectively) We assume that the solid part of the poroelastic material is glasswool of type R, with the following properties: Young modulus 8.7×1010 N/m2 Poisson coefficient 0.15 Density 2500 kg/m3 Assuming that the poroelastic material is completely saturated by air, we have obtained the following macroscopic coefficients for generalized Biot model (3.3)-(3.4): A=⎛ ⎝0.7576 0.3222E-4 0.8560E-4 0.8526E-4 0.7581 0.8646E-4 0.5362E-4 0.2928E-4 0.7573 ⎞ ⎠, ˆc= -0.4955E-11,
3.7. Numerical results 55 B=⎛ ⎝-0.2050 -0.1385E-3 -0.5167E-3 -0.1385E-3 -0.2061 0.4572E-3 -0.5167E-3 0.4572E-3 -0.2047 ⎞ ⎠, A=10 11 ×⎛ ⎝M11 M12 M13 M12 M22 M23 M13 M23 M33 ⎞ ⎠, where M11 =⎛ ⎝0.1589 -0.3337E-4 -0.8505E-4 -0.3337E-4 0.1195E-1 -0.1311E-3 -0.8505E-4 -0.1311E-3 0.1213E-1 ⎞ ⎠, M12 =⎛ ⎝-0.1057E-4 0.1813E-1 -0.5855E-4 0.1813E-1 -0.2185E-4 -0.1051E-5 -0.5855E-4 -0.1051E-5 0.2270E-4 ⎞ ⎠, M13 =⎛ ⎝0.1115E-5 -0.1254E-4 0.1834E-1 -0.1254E-4 0.8837E-6 0.1045E-5 0.1834E-1 0.1045E-5 0.4812E-5 ⎞ ⎠, M22 =⎛ ⎝0.1190E-1 -0.1144E-3 0.2185E-3 -0.1144E-3 0.1586 0.7522E-5 0.2185E-3 0.7522E-5 0.1204E-1 ⎞ ⎠, M23 =⎛ ⎝0.1190E-5 0.4927E-4 -0.1044E-3 0.4927E-4 -0.1069E-5 0.1851E-1 -0.1044E-3 0.1851E-1 -0.1966E-5 ⎞ ⎠, M33 =⎛ ⎝0.1212E-1 0.5964E-2 -0.1349E-3 0.5964E-2 0.1203E-1 -0.2908E-3 -0.1349E-3 -0.2908E-3 0.1584 ⎞ ⎠. Moreover, porosity of the porous sample is φ=0.648. Finally, by progressive refinement of meshes one can see that matrices of macroscopic parameters Aand Btend to isotropic matrices. Similarly, tensor Ahas the same structure as the elasticity tensor for isotropic elastic materials. Thus, due to the symmetry properties of the cell, the macroscopic behavior of the studied porous medium is isotropic. 3.7 Numerical results In this Section we present some numerical results obtained with a computer code developed by us which implements the numerical methods proposed in this chapter. This code allows us to compute the response diagram of the enclosure shown in Figure 3.1, consisting of a fluid and a poroelastic medium.
56 FEM for poroelastic models in acoustics In order to validate our method, we are going to build a simple example which can be reduced to a one-dimensional problem and then solved exactly. We recall that equations satisfied by pressure and displacement fields in the poroelastic media are ω2(ρI−ρFA)uA+ div (A[D(uA)])+(A+B−φI)grad pA=0in ΩA, −ω2ˆcpA+1 ρF div(Agrad pA)=ω2div ((A+B−φI)uA)inΩ A. If we assume that every linear operator is a multiple of the identity operator, we can find a solution (uA,p A) of the form uA(x1,x 2,x 3)=uA(x3)e3,pA(x1,x 2,x 3)=pA(x3), and rewrite the above three-dimensional problem as a one-dimensional problem, namely, ω2(ρ−ρFa)uA+su A+(a+b−φ)p A= 0 in (0,a A), −ω2ˆcpA+a ρF p A=ω2(a+b−φ)u Ain (0,a A), where the prime denotes derivative with respect to zand we have supposed that A[E(uA)] e3= su Ae3,A=aI,B=bI. After some algebraic manipulations, we obtain a ω2ρF p A−ˆc−(a+b−φ)2 s−a(ρ−ρFa) ω2sρFp A−ω2(ρ−ρFa)ˆc spA= 0 in (0,a A), uA=s ω2(a+b−φ)(ρ−ρFa)−a ω2ρF p A+ˆc−(a+b−φ)2 sp Ain (0,a A). Let us assume a similar assumption for fluid displacement and interface pressure, i.e., uF(x1,x 2,x 3)=uF(x3)e3and pF(x1,x 2,x 3)=pF(x3). We also suppose that ΩF=(0,b)× (0,d)×(˜aF,0), being ˜aF=−aF, and ΩA=(0,b)×(0,d)×(0,a A). Then the coupled fluid-poroelastic problem can be written as −ω2ρFuF+p F=0in(˜aF,0),(3.36) pF=−ρFc2u Fin (˜aF,0),(3.37) a ω2ρF p A−ˆc−(a+b−φ)2 s−a(ρ−ρFa) ω2sρFp A−ω2(ρ−ρFa)ˆc spA= 0 in (0,a A),(3.38) uA=s ω2(a+b−φ)(ρ−ρFa)−a ω2ρF p A+ˆc−(a+b−φ)2 sp Ain (0,a A),(3.39) −pF(0) = su A(0)+(a+b−φ)pA(0),(3.40) uF(0) = uA(0),(3.41) uF(˜aF)=0,(3.42) uA(aA)=0,(3.43) ap A(aA)=0,(3.44) ap A(0) = 0,(3.45) uF(˜aF)=f. (3.46)
3.8. Conclusions 57 The general solution of the ordinary differential system (3.36)-(3.38) is uF(x3)=C1e−ikFx3+C2eikFx3,x 3∈(−aF,0),(3.47) pA(x3)=C3e−ikA,1x3+C4eikA,1x3+C5e−ikA,2x3+C6eikA,2x3,x 3∈(0,a A),(3.48) where kF=ω cand {kA,1,−kA,1,k A,2,−kA,2}are the four roots of the polynomial equation −a ω2ρF k4 A+ˆc−(a+b−φ)2 s−a(ρ−ρFa) ω2sρFk2 A+ω2(ρ−ρFa)ˆc s=0. If we take into account boundary and interface conditions (3.40)-(3.46) and expressions (3.47) and (3.48), amplitudes Cj,1≤j≤6 can be calculated by solving a linear system of equations. We have considered that fluid is air with ρF=1.225 kg/m3and c= 343 m/s, whereas properties of the porous material are summarized in s=9.18633 ×1010N/m2,φ=0.95, a=0.67857, b=−0.05, ˆc=−6.59172 ×10−6ms2/kg and ρ=1.26163 ×102kg/m3. With respect to the dimensions of the enclosure, length and width are 1 m while height is 1 m for the first layer of free fluid and 1 m for the second layer of porous material whereas the normal displacement on ΓEis f= 60. We have computed the solution to this problem with three different uniform meshes, named mesh 1, mesh 2 and mesh 3 of 2548, 8140 and 18788 degrees of freedom, respectively. In Figure 3.5, we show the L2-norm of the relative errors for fluid displacement, uh F−uF0,ΩF/uF0,ΩF, porous displacement, uh A−uA0,ΩA/uA0,ΩA, porous pressure, ph A−pA0,ΩA/pA0,ΩA, and interface pressure, ph F−pF0,ΩF/pF0,ΩF, against meshsize, h. As it can be seen, convergence of order 2 is achieved for poroelastic fields and interface pressure. In addition, convergence of order 1 is achieved for fluid displacement. As a real life test, we are going to compute the solution of the coupled problem, using the data obtained in the previous Section by solving cell problems. Figure 3.6 shows the response curves for enclosure in Figure 3.1, when solved with mesh 1 having 2548 degrees of freedom. In this curve the logarithm of “energy”, log10 1 2ΩF ρFc2|div uF|2+ΩA A[E(uA)] : E(¯ uA)+ΩA div ((A+B−φI)uA)¯pA, is plotted for angular frequencies ωranging from 50 to 2000 rad/s. Several response peaks can be observed in this curve. In fact, response peaks of the curve shows the resonance frequencies for the coupled system shown in Figure 3.1. 3.8 Conclusions We have considered a mathematical model for acoustic propagation in periodic nondissipative porous media with elastic solid frame and open pore. Parameters of this model have been computed by solving some partial differential equations in the unit cell obtained
58 FEM for poroelastic models in acoustics 0.125 0.175 0.25 10−3 10−2 10−1 h (mesh−size) Relative error Relative error in the fluid displacement Relative error Order h 0.125 0.175 0.25 10−5 10−4 10−3 h (mesh−size) Relative error Relative error in the porous displacement Relative error Order h2 0.125 0.175 0.25 10−5 10−4 10−3 h (mesh−size) Relative error Relative error in the porous pressure Relative error Order h2 0.125 0.175 0.25 10−5 10−4 10−3 h (mesh−size) Relative error Relative error in the interface pressure Relative error Order h2 Figure 3.5: Curves of convergence for fluid and porous fields. by homogenization methods. Then a three-dimensional finite element method has been proposed and implemented for numerical solution of the coupling between a fluid and the above porous medium. In order to validate the proposed methodology and to assess convergence properties, the computer code has been used for a test example having analytical solution. Then, as an application, we have computed the response curve for an enclosure containing air and a layer of porous material.
3.8. Conclusions 59 Figure 3.6: Response curve.
66 Chapter 4. A non reflecting porous material: the Perfectly Matched Layers or when a singular absorbing function is used in the construction of the perfectly matched layer. Finally, in Section 4.4 we consider a simple problem: the propagation of plane waves with oblique incidence in a two-dimensional unbounded domain. We show that a PML method based on a non-integrable absorbing function allows recovering the exact solution in the domain of interest. 4.2 Wave equation For the sake of completeness, the derivation of the two-dimensional linear wave equation is outlined, since the same steps will be done for the PML equation. 4.2.1 Time domain equations First, we state the system of equations of the wave motion in terms of the pressure and the velocity fields. Find the velocity field Vand pressure field Pin the whole space R2which satisfy the system of equations: ∂ρ ∂t +ρ0div V=0,(4.1) ρ0 ∂V ∂t +grad P=0,(4.2) P=˜ P(ρ),(4.3) where ρ0is the fluid mass density at rest, which is supposed to be a constant, and ρis the mass density. From the fluid constitutive equation (4.3), we deduce that ∂P ∂t =c2∂ρ ∂t, where cis the sound speed. By using (4.1), we obtain 1 c2 ∂P ∂t +ρ0div V=0.(4.4) If we derive with respect to the time variable, then we have 1 c2 ∂2P ∂t2+ρ0div ∂V ∂t =0. Now, by applying the divergence operator in (4.2), ρ0div ∂V ∂t +ΔP=0,
4.3. Cartesian Perfectly Matched Layers 67 and replacing this expression in (4.4), it results 1 c2 ∂2P ∂t2−ΔP=0.(4.5) 4.2.2 Time-harmonic equations We take into account a simple toy problem which illustrates the construction of perfectly matched layers in Cartesian coordinates. If there exists an acoustical source in the fluid domain, a non-null right-hand term Farises in the above equation. If we assume that Fis time-harmonic, i.e., F(x,t)=Re(e −iωtf(x)), then P(x,t)=Re(e −iωtp(x)) is a solution of equation (4.5) provided −ω2 c2p−Δp=f. (4.6) Moreover, from a physical point of view, if we postulate that no waves are reflected from infinity, pmust satisfy the Sommerfeld’s radiation condition uniformly in all the directions (see [72]). Hence, the source problem stated in the whole space R2and using the pressure field as unique unknown is: Given an acoustical source f, find the pressure field psatisfying −ω2 c2p−Δp=fin R2,(4.7) lim r→∞ r1 2∂p ∂r −ikp=0,(4.8) where r=|x|and k=ω cis the wave number. If the harmonic excitation inside the fluid domain is a monopole supported in d= (d1,d 2), then the right-hand term is f=−iωρ0Qδd, where Qis the volume velocity and δd is the Dirac’s delta (see [63]). In this case the solution of problem (4.7)-(4.8) is p(x)=ωρ0Q 4H(1) 0(k|x−d|),(4.9) where H(1) 0is the Hankel function of first kind and order zero. 4.3 Cartesian Perfectly Matched Layers Now we are going to introduce the Cartesian perfectly matched layers to deal with the same problem. As a first step, we assume that the domain filled by the PML is unbounded. In spite that this fact does not avoid to state a problem in a bounded domain, it allows us to illustrate the derivation of the PML problem.
68 Chapter 4. A non reflecting porous material: the Perfectly Matched Layers 4.3.1 Time-domain equations We assume that the equations of linear acoustics (4.1)-(4.3) are satisfied in ΩF=(−a, a)×(−b, b), whereas the perfectly matched layers are situated in ΩA=R2\ΩF. The two-dimensional Cartesian PML equations involve absorbing coefficients σjdefined for each 1 ≤j≤2, such that they are monotonically increasing, non negative, smooth in ΩAand null inside the fluid domain ΩF. Moreover σjonly depends on the spatial variable xj. Next, we state which are the equations governing the pressure field inside ΩA. With this aim, we use the ‘splitting’ technique developed originally by Berenger ([22], [72]). First, we suppose that the pressure field in ΩAis the addition of two terms which involve two new unknowns without any physical meaning, Pj,1≤j≤2, such that Eq. (4.1) can be rewritten in terms of the pressure and the velocity fields as ∂Pj ∂t +ρ0c2∂Vj ∂xj =0,1≤j≤2,(4.10) P=P1+P2, where Vjdenotes the j-th component of the velocity field. Analogously, if we write every component of Eq. (4.2), we obtain ρ0 ∂Vj ∂t +∂P ∂xj =0,1≤j≤2.(4.11) Inside the domain ΩA, where the PML is situated, we rewrite the equations adding a dissipative or damping term involving the fictitious pressure component Pjand the velocity field Vj, modifying (4.10)-(4.11). Consequently, the system to be solved in ΩAis ∂Pj ∂t +σjPj+ρ0c2∂Vj ∂xj =0,1≤j≤2,(4.12) ρ0∂Vj ∂t +σjVj+∂P ∂xj =0,1≤j≤2,(4.13) P= 2 j=1 Pj.(4.14) 4.3.2 A physical interpretation If we denote by σthe diagonal tensor of coefficients σj,j=1,2, associated with the j-th component, equation (4.13) can be rewritten as ρ0∂V ∂t +σV =−grad P, (4.15)
4.3. Cartesian Perfectly Matched Layers 69 which is equal to the equation introduced in Chapter 1 for the Darcy’s like model in the porous media (where σwas the flux resistivity tensor). The unique difference between that model and the system of equations of the PML is the absorbing functions introduced for the pressure field in (4.12). Integrating (4.12) in time, we obtain Pj=−ρ0c2t 0 eσj(s−t)∂Vj ∂xj ds. Adding P1and P2, we obtain an explicit expression for the pressure P, P=−ρ0c2t 0 2 j=1 eσj(s−t)∂Vj ∂xjds. (4.16) If we derive the previous equality two times respect to the time variable, we obtain ∂2P ∂t2=−ρ0c2'div ∂V ∂t +σV +t 0 2 j=1 σ2 jeσj(s−t)∂Vj ∂xjds −2 2 j=1 σj ∂Vj ∂xj(.(4.17) Finally, if we take into account (4.15)-(4.17), we obtain a time-formulation of the PML in terms of the pressure field, ∂2P ∂t2−c2ΔP−c2t 0 2 j=1 σ2 jeσj(s−t)∂ ∂xjs 0 eσj(τ−s)∂P ∂xj dτds +2c2 2 j=1 σj ∂ ∂xjt 0 eσj(s−t)∂P ∂xj ds=0. As a particular case, if the absorbing functions are constant and equal, i.e., if σ1(x1)= σ2(x2)=σ, in the domains where they are not null, then we can rewrite the previous PML equation in terms of the pressure field in a simpler form. In fact, the system of equations (4.12)-(4.14) leads to ρ0∂V ∂t +σV=−grad P, (4.18) ∂P ∂t +σP =−ρ0c2div V.(4.19) Integrating (4.18) with respect to the time variable, we have V=−1 ρ0t 0 eσ(s−t)grad P ds, (4.20)
70 Chapter 4. A non reflecting porous material: the Perfectly Matched Layers and, substituting this expression in (4.19) and deriving with respect to time, ∂2P ∂t2+σ∂P ∂t −c2ΔP+ div c2σt 0 eσ(s−t)grad Pds =0.(4.21) As a conclusion, we observe that the PML equation written in terms of the pressure has a damping term depending linearly on σand, as in the case of materials with memory, the pressure field at time t=Tdepends on the pressure fields in the interval 0 ≤t<T. Moreover, let us remark that, obviously, if σ= 0 then we recover the wave equation. 4.3.3 Time-harmonic equations Now we rewrite the system of equations (4.12)-(4.14) in the frequency domain. If we assume that P(x,t) = Re(p(x)e −iωt),P j(x,t) = Re(pj(x)e −iωt),V(x,t) = Re(v(x)e −iωt), then the PML equations (4.12)-(4.14) turn into pj=−ρ0c2 σj−iω ∂vj ∂xj ,(4.22) vj=1 ρ0(σj−iω) ∂p ∂xj ,(4.23) p=p1+p2.(4.24) By substituting (4.23) in (4.22), pj=−c21 σj−iω ∂ ∂xj1 σj−iω ∂p ∂xj ,1≤j≤2, and then, using this expression in (4.24), we obtain −ω2 c2p− 2 j=1 iω σj−iω ∂ ∂xjiω σj−iω ∂p ∂xj =0.(4.25) Following the ideas introduced in [47], we can define the complex change of variable ˆxj(xj)=xj+i ωxj 0 σj(s)ds, 1≤j≤2,(4.26) so that formally we recover the weights that arise in (4.25) using its derivatives, i.e., ∂xj ∂ˆxj =−iω σj(xj)−iω.
4.3. Cartesian Perfectly Matched Layers 71 Hence the PML equation (4.25) can be seen as the Helmholtz equation in a new complex coordinate system ˆ x=(ˆx1,ˆx2), since we can rewrite Eq. (4.25) formally as −ω2 c2p− 2 j=1 ∂2p ∂ˆx2 j =0. So, the construction of the PML equation can be understood as a complex stretching of coordinates. However, let us remark that we will not use this formal property for any derivation of the theoretical results on PML through this thesis. Finally, we introduce some notation for the weights of the PML equation that will be used through the next chapters. We define γ(xj),1≤j≤2as γj(xj):=σj(xj)−iω −iω =1+ i ωσj(xj). Now, if we assume an external harmonic source, F(x,t) = Re(f(x)e −iωt), such that supp f⊂ΩF, the equation to solve in the PML domain as well as in the fluid domain is −ω2 c2p− 2 j=1 1 γj ∂ ∂xj1 γj ∂p ∂xj =fin R2.(4.27) Infinite thickness layer with a bounded absorbing function In order to close the system of equations of the time-harmonic source problem with PML layers, (4.27), we need a radiation boundary condition, uniform in all the directions, which substitutes the classical Sommerfeld’s condition. We take a homogeneous Dirichlet condition at infinity. Hence, the source problem stated in the whole space R2written in terms of the pressure field is the following: Given an acoustic source f, with compact support contained in ΩF,find the pressure field pthat satisfies −k2p− 2 j=1 1 γj ∂ ∂xj1 γj ∂p ∂xj =fin R2,(4.28) lim |x|→∞ p(x)=0,(4.29) where k=ω/c is the wavenumber. For example, if we suppose that the source is a monopole supported at the point d= (d1,d 2)∈ΩF, then f=−iωρ0Qδd. Straightforward computations show that a solution of (4.28)-(4.29) is p(x)=ωρ0Q 4H(1) 0(kˆrd(x)),(4.30)
72 Chapter 4. A non reflecting porous material: the Perfectly Matched Layers where ˆrd(x)=(ˆx1(x1)−d1)2+(ˆx2(x2)−d2)2. From the complex stretching of coordinates, since ˆx1(x1)=x1if x1∈(−a, a) and ˆx2(x2)=x2if x2∈(−b, b), it is immediate to check that the fundamental solutions given by (4.9) and (4.30) coincide in the fluid domain ΩF. In general, the exact solution of the original scattering problem is recovered from the PML problem if the physical domain of interest is surrounded by an unbounded PML layer (see [52]). Finite-thickness layer with a singular absorbing function Now we consider the more realistic case of a bounded PML domain. We assume that the linear acoustic equations (4.1)-(4.3) are satisfied in ΩF=(−a, a)×(−b, b), whereas the PML is situated in the bounded domain ΩA=[−a,a ]×[−b,b ]\ΩF, where a<a and b<b . We denote by ΓD=∂ΩA\∂ΩFthe exterior artificial boundary where the PML domain is truncated (see Figure 4.1). −b∗ b −aa ∗ b∗ −b −a∗a ΓD ΩA ΩF x2 x1 Figure 4.1: Cartesian PML stated in a two-dimensional bounded domain. In this case the analogous to problem (4.28)-(4.29) is the following: Given an acoustic source fwith compact support in ΩF,find the pressure filed pthat satisfies −k2p− 2 j=1 1 γj ∂ ∂xj1 γj ∂p ∂xj =fin ΩA∪ΩF,(4.31) p=0 on ΓD.(4.32)
4.4. Plane wave analysis of the PMLs 73 It is well-known that since the PML domain has been truncated, spurious reflections arise in the solution of the classical PML problem with standard bounded absorbing functions. So it does not coincide any more with the original solution of the Helmholtz problem stated in the unbounded domain. This drawback can be avoided by using singular absorbing functions. More precisely, we state the same two-dimensional Cartesian PML equations involving the absorbing coefficients σj. But now, besides assuming that the absorbing functions are non negative, smooth, monotonically increasing and null in ΩF, we suppose that they are not integrable, so that lim |x1|→|a|σ1(x1)=+∞, lim |x2|→|b|σ2(x2)=+∞, Through the rest of this part of the thesis, we will focus our attention on PMLs based on this kind of singular absorbing functions. In this case, using the same source as in the previous section, namely f=−iωρ0Qδd, then it is easy to check that a solution of problem (4.31)-(4.32) is p(x)=ωρ0Q 4H(1) 0(kˆrd(x)),(4.33) where ˆrd(x)=(ˆx1(x1)−d1)2+(ˆx2(x2)−d2)2. As in the previous case where the PML domain was unbounded, since ˆx1(x1)=x1if x1∈(−a, a) and ˆx2(x2)=x2if x2∈(−b, b), it is immediate to check that the fundamental solutions given by (4.9) and (4.33) coincide in the fluid domain ΩF. This kind of results are proved theoretically in Chapter 6 (see also [33]) for the PML equation written in radial coordinates for time-harmonic scattering problems. 4.4 Plane wave analysis of the PMLs We consider a simple problem which will provide valuable information for the design of an efficient PML method: the propagation of two-dimensional acoustic plane waves with oblique incidence. With this problem, we also illustrate the ideas introduced in the previous subsection. Consider the following time-harmonic problem posed in the right half-space: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Δp+k2p=0,x 1>0,x 2∈R, p(0,x 2)= e ik2x2,x 2∈R, lim x1→+∞∂p ∂x1−ik1p=0, (4.34) where k=ω/c is the wave number, k1=kcos θand k2=ksin θ, with θbeing the incidence angle. It is well-known that the solution of this problem is the plane wave p(x1,x 2)= e i(k1x1+k2x2).
74 Chapter 4. A non reflecting porous material: the Perfectly Matched Layers We introduce a PML in the vertical strip a<x 1<a ∗, to truncate the unbounded domain in the x1-direction (see Figure 4.2). The strip 0 <x 1<ais the so called ‘physical domain’, i.e., the domain where we are interested in computing the solution of problem (4.34). θ x1=0 x1=ax 1=a∗ PML Figure 4.2: PML in the x1-direction for plane waves with oblique incidence. We consider a variable absorption coefficient σin the PML. This coefficient is allowed to be a function of the variable x1; constant, linear or quadratic ‘absorbing functions’ are the typical choices (see, for instance, [16, 22, 52]). In this case, we allow for any arbitrary non-negative absorbing function. In order to solve the problem with PML layers, we distinguish two pressure fields: we denote by pFand pAthe restriction of the pressure field to the fluid domain and to the PML, respectively. This approach allows us to write explicitly the spurious reflection that the bounded PML domain produces in the fluid domain. Thus, the amplitudes of the pressure waves in the physical domain, pF, and in the PML, pA, are the solution of the following equations: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ΔpF+k2pF=0,0<x 1<a, 1 γ1 ∂ ∂x11 γ1 ∂pA ∂x1+∂2pA ∂x2 2 +k2pA=0,a<x 1<a ∗, pF(0,x 2)= e ik2x2, pF(a, x2)=pA(a, x2), ∂pF ∂x1 (a, x2)= 1 γ1(a) ∂pA ∂x1 (a, x2), pA(a∗,x 2)=0, with, as it was defined in Section 4.3.3, γ1(x1):=1,if 0 <x 1<a, 1+ i ωσ1(x1),if a≤x1<a ∗. As we have already shown in the previous subsection, if we introduce the complex change of variable )x1(x1)=x1 0 γ1(s)=x1+i ωx1 a σ1(s),x 1∈[a, a∗),(4.35)
4.4. Plane wave analysis of the PMLs 75 then ∂)x1 ∂x1 =γ1and ∂ ∂)x1 =1 γ1 ∂ ∂x1 , and hence, if we denote ˆpA()x1,x 2)=pA(x1,x 2), then we have 1 γ1 ∂ ∂x11 γ1 pA ∂x1+∂2pA ∂x2 2 +k2pA=0 ⇐⇒ ∂2ˆpA ∂)x2 1 +∂2ˆpA ∂x2 2 +k2ˆpA=0. Since we have recovered formally the Helmholtz equation, the solution of the PML problem can be written as superposition of plane waves: #pF(x1,x 2)=Ieik1x1+Re−ik1x1eik2x2,x 1∈(0,a), ˆpA()x1,x 2)=Teik1)x1+RAe−ik1)x1eik2x2,x 1∈[a, a∗), where Iis the amplitude of the incident wave, Tthat of the wave transmitted to the PML, and Rand RAare the amplitudes of the reflected waves in the physical domain and in the absorbing layer, respectively. By substituting (4.35) in the last equation, we can write the solution in the absorbing layer in the following equivalent form: pA(x1,x 2)=*Teik1x1e−cos θ c!x1 aσ1(s)+RAe−ik1x1ecos θ c!x1 aσ1(s)+eik2x2. Hence, from the boundary condition at x1= 0, we obtain I=1−R, and, from the transmission conditions at x1=a, R=RAand I=T. Notice that the latter implies that no spurious reflections arise at x1=a(which is the main feature of the PML technique); the terms involving Rand RAarise as a consequence of the waves reflected at x1=a∗. Finally, the homogeneous Dirichlet boundary condition at x1=a∗yields R=RA=e2ik1a∗ e2ik1a∗−e2 cos θ c!a∗ aσ1(s).(4.36) Summarizing, we have obtained the following analytical expression for the solution to the PML problem above: ⎧ ⎪ ⎨ ⎪ ⎩ pF(x1,x 2)=(1 −RA)e ik1x1+RAe−ik1x1eik2x2,x 1∈(0,a), pA(x1,x 2)=*(1 −RA)e ik1x1e−cos θ c!x1 aσ1(s)+RAe−ik1x1ecos θ c!x1 aσ1(s)+eik2x2, x1∈[a, a∗).
82 Chapter 5. An optimal PML in Cartesian coordinates are solution of the following equations (see, for instance, [52]): ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ΔpF+k2pF= 0 in ΩF, 1 γ1 ∂ ∂x 1 γ1 ∂pA ∂x +1 γ2 ∂ ∂y 1 γ2 ∂pA ∂y +k2pA= 0 in ΩA, ∂pF ∂n=gon Γ, pF=pAon ΓI, ∂pF ∂νx +∂pF ∂νy =1 γ1 ∂pA ∂νx +1 γ2 ∂pA ∂νy on ΓI, pA= 0 on ΓD, (5.2) where γ1(x1)=1,if |x1|<a, 1+ i ωσ1(|x1|),if a≤|x1|<a ∗, and γ2(x2)=1,if |x2|<b, 1+ i ωσ2(|y|),if b≤|x2|<b ∗. The main goal of this chapter is to determine how to choose the absorbing functions σ1 and σ2, so that pFbe an approximation as close as possible to the solution pof problem (5.1) in the physical domain. According to the results of the previous section, the natural candidates are unbounded functions σ1and σ2such that a∗ a σ1(s)=+∞and b∗ b σ2(s)=+∞. 5.3 Finite element discretization. In this section we describe a finite element method for the numerical solution of (5.2) and show that the resulting discrete problem is well posed only for certain unbounded absorbing functions. This will lead to additional constraints on σ1and σ2. We consider a partition in triangles of the physical domain ΩFand a partition in rectangles of the absorbing layer ΩA, matching on the common interface ΓIas shown in Figure 5.3. As usual, hdenotes the mesh-size. The reason why we use such hybrid meshes is that triangles are more adequate to fit the boundary of the obstacle, whereas rectangles will allow us to compute explicitly the integrals involving the absorbing function appearing in the elements in the layer. This is not strictly necessary, since these integrals can also be efficiently computed by means of standard quadrature rules as shown in the appendix. However, in this chapter, we will mainly consider exact integration to be able to assess the accuracy of the proposed PML independently of quadrature errors.
5.3. Finite element discretization. 83 ΓD ΓI ΩA ΩF Figure 5.3: Hybrid mesh on PML and physical domain. We will compute approximations ph Fand ph Aof the pressure amplitudes in the physical domain and in the absorbing layer, respectively, by using linear triangular finite elements for the former and bilinear rectangular finite elements for the latter. The degrees of freedom defining the finite element solution are the values of ph Fand ph Aat the vertices of the elements. Notice that because of the transmission condition pF=pAon ΓI, the values of ph Fand ph A must coincide at the vertices on the interface. Moreover, we impose the Dirichlet boundary condition ph A= 0 on ΓDto the finite element solution. Hence, ph Adoes not have degrees of freedom on the outer boundary. This fact will be essential for the resulting discrete problem to be well posed. Standard arguments in this finite element framework lead to the following discretization of problem (5.2)written in weak form: ΩF∇ph F·∇¯qhdx1dx2−ΩF k2ph F¯qhdx1dx2+ΩA γ2 γ1 ∂ph A ∂x1 ∂¯qh ∂x1 dx1dx2 +ΩA γ1 γ2 ∂ph A ∂x2 ∂¯qh ∂x2 dx1dx2−ΩA k2γ1γ2ph A¯qhdx1dx2=Γ g¯qh, for all function qh, continuous in ΩF∪ΩA, piecewise linear in ΩF, piecewise bilinear in ΩA, and vanishing on ΓD. Once the discrete problem is written in matrix form, it yields a system of linear equations whose unknowns are the nodal values of ph Fand ph A. The entries of the matrix are computed by assembling the element matrices; in particular, the following terms involve the unbounded absorbing functions: K γ2 γ1 ∂Ni ∂x1 ∂Nj ∂x1 dx1dx2,K γ1 γ2 ∂Ni ∂x2 ∂Nj ∂x2 dx1dx2,and K k2γ1γ2NiNjdx1dx2,(5.3) with Kbeing a rectangular element in ΩAand {Ni}the nodal finite element basis. For the discrete problem to be well posed, it is necessary that all the integrals above be finite, what is not trivial since they involve singular functions, whenever Kis a rectangle with one edge lying on the outer boundary ΓD. For instance, consider the element Kshown in Figure 5.4 (the forthcoming arguments and conclusions hold also true for all other elements with edges lying on ΓD). Notice that,
84 Chapter 5. An optimal PML in Cartesian coordinates ΓD b2 a∗ b1 a∗−h1 h2 h1 K P2 P1 Figure 5.4: Finite element with an edge lying on ΓD. since ph Avanishes at the vertices on ΓD, in this element we only need to compute the integrals (5.3) for the nodal functions N1and N2associated with the vertices denoted by P1and P2, respectively. These functions are given by N1(x1,x 2)=(x1−a∗)(x2−b2) h1h2 ,N 2(x1,x 2)=−(x1−a∗)(x2−b1) h1h2 , and their partial derivatives by ∂N1 ∂x1 (x1,x 2)=x2−b2 h1h2 ,∂N2 ∂x1 (x1,x 2)=−x2−b1 h1h2 , ∂N1 ∂x2 (x1,x 2)=x1−a∗ h1h2 ,∂N2 ∂x2 (x1,x 2)=−x1−a∗ h1h2 . Therefore, the integrals in (5.3) can be written as follows: K γ2 γ1 ∂Ni ∂x1 ∂Nj ∂x1 dx1dx2=±1 h2 1h2 2b2 b1 γ2(x2)(x2−bj)(x2−bi)dx2a∗ a∗−h1 dx1 γ1(x1),(5.4) K γ1 γ2 ∂Ni ∂x2 ∂Nj ∂x2 dx1dx2=±1 h2 1h2 2b2 b1 dx2 γ2(x2)a∗ a∗−h1 γ1(x1)(x1−a∗)2dx1,(5.5) K γ1γ2NiNjdx1dx2 =−1 h2 1h2 2b2 b1 γ2(x2)(x2−bj)(x2−bi)dx2a∗ a∗−h1 γ1(x1)(x1−a∗)2dx1.(5.6) We assume singularities of power type for the absorbing functions: σ1(x1)=O(a∗−x1)−αas x1→a∗,σ 2(x2)=O(b∗−x2)−αas x2→b∗.(5.7) Notice that the constraint that σ1and σ2have unbounded integrals holds true if and only if α≥1.
5.4. Determination of the absorbing function 85 From the definitions of γ1and γ2we have that γ1(x1)=O(a∗−|x1|)−αand γ2(x2)= O(b∗−|x2|)−α. Moreover, |γ1|≥1 and |γ2|≥1 and, hence, the integrals of 1/γ1(x1) and 1/γ2(x2) are always finite. For an element Kas that in Figure 5.4, γ2is bounded in the interval [b1,b 2] and, consequently, the integrals involving γ2(x2) are also finite. Finally, for the integral involving γ1(x1)wehave a∗ a∗−h γ1(x1)(x1−a∗)2dx1=a∗ a∗−hO(a∗−x1)2−αdx1<∞⇐⇒α<3. As a consequence of this analysis, we will restrict our choice of σ1and σ2satisfying (5.7) to exponents αsuch that 1 ≤α<3. 5.4 Determination of the absorbing function In this section we report the numerical experimentation performed to determine the most convenient unbounded absorbing functions. With this purpose, we have applied our PML method with different σ1and σ2to a scattering problem with known analytical solution and compared the accuracy of the numerical results. Consider problem (5.1) where the obstacle Ω is the unit circle centered at the origin. Given any inner point (x0 1,x 0 2) of this circle, it is well known that the function p(x1,x 2)= i 4H(1) 0k"(x1−x0 1)2+(x2−x0 2)2 satisfies the first and third equations of (5.1). Therefore, if we take g=∂p/∂n, then pis the unique solution of this problem. In our experiments we have taken x0 1=0.5m, x0 2= 0, and k=ω/c, with c= 340 m/s and different values of the frequency ω. For our computational domain we have taken a=b=2.0 m and a∗=b∗=2.25 m (see Figure 5.5). We have used uniform refinements of the mesh shown in Figure 5.5; the number Nof elements through the thickness of the PML is used to label each mesh. To measure the accuracy we have computed the relative error in the L2-norm in ΩF: Error = ΩFph F−p2dx1dx21/2 ΩF|p|2dx1dx21/2,(5.8) where ph Fis the numerical solution in the physical domain and pthe exact solution. According to the results of the previous section, it is enough to restrict the analysis to absorbing functions satisfying (5.7) with 1 ≤α<3. We have considered the integer powers: α= 1 and α= 2. In particular, we have tested functions of the following type, where βis a free parameter to be fitted:
86 Chapter 5. An optimal PML in Cartesian coordinates x y 1m x0=0.5m a∗=2.25 m a=2m b∗=2.25 m b=2m N=2 Figure 5.5: Domains and mesh in the scattering problem. •A. σ1(x1)= β a∗−x1 ,σ 2(x2)= β b∗−x2 ; •B. σ1(x1)= β (a∗−x1)2,σ 2(x2)= β (b∗−x2)2. Notice that, in both cases, σ1(a)= 0 and σ2(b)= 0. Hence, the corresponding coefficients γ1and γ2will be discontinuous. To avoid eventual side effects of these discontinuities in the coupling conditions on ΓI, we have also considered functions of the following type, which yield continuous γ1and γ2: •C. σ1(x1)= β a∗−x1−β a∗−a,σ 2(x2)= β b∗−x2−β b∗−b, •D. σ1(x1)= β (a∗−x1)2−β (a∗−a)2,σ 2(x2)= β (b∗−x2)2−β (b∗−b)2. In each case, we have fitted the parameter βso as to minimize the error. Figs. 5.6 to 5.11 show the results obtained with each type of absorbing functions and a range of values of β. We have used three meshes with refinement levels N= 2, 4, and 8, which have 464, 1720, and 6768 vertices, respectively. We report the results obtained with two frequencies: ω= 250 rad/s and ω= 750 rad/s. We report, in Table 5.1, the minimal relative errors and the optimal values of βdetermined for each type of absorbing function and each of the three meshes. It can be clearly seen from this table that the smallest errors are always attained for a function of Atype
5.4. Determination of the absorbing function 87 0 5 10 15 20 0 10 20 30 40 50 60 β/c Relative error (%) Type A Type B Type C Type D Figure 5.6: Relative errors for PML with different unbounded absorbing functions. Mesh: N=2;ω= 250 rad/s. 0 5 10 15 20 0 5 10 15 20 25 30 35 40 β/c Relative error (%) Type A Type B Type C Type D Figure 5.7: Relative errors for PML with different unbounded absorbing functions. Mesh: N=2;ω= 750 rad/s. 0 2 4 6 8 10 12 0 5 10 15 β/c Relative error (%) Type A Type B Type C Type D Figure 5.8: Relative errors for PML with different unbounded absorbing functions. Mesh: N=4;ω= 250 rad/s. 0 2 4 6 8 10 12 0 2 4 6 8 10 β/c Relative error (%) Type A Type B Type C Type D Figure 5.9: Relative errors for PML with different unbounded absorbing functions. Mesh: N=4;ω= 750 rad/s. with the parameter β≈c. To allow for comparison, we also include in the table the errors for this choice; namely, σ1(x1)= c a∗−x1 ,σ 2(x2)= c b∗−x2 .(5.9) Let us remark that, for each mesh, the CPU time needed to solve the problem is essentially the same for the four types of absorbing functions. The condition numbers of the system matrices remain basically of the same order of magnitude for all of the choices, too. As a definite conclusion of this experimentation, we propose to use the absorbing functions (5.9). Notice that an additional advantage of this proposal is that the resulting PML method does not need any parameter to be determined.
88 Chapter 5. An optimal PML in Cartesian coordinates 0 2 4 6 8 10 12 0 1 2 3 4 5 6 7 8 9 β/c Relative error (%) Type A Type B Type C Type D Figure 5.10: Relative errors for PML with different unbounded absorbing functions. Mesh: N=8;ω= 250 rad/s. 0 2 4 6 8 10 12 0 1 2 3 4 5 6 β/c Relative error (%) Type A Type B Type C Type D Figure 5.11: Relative errors for PML with different unbounded absorbing functions. Mesh: N=8;ω= 750 rad/s. Table 5.1: Minimal errors and optimal values of the parameter βfor PML with different unbounded absorbing functions. ω= 250 rad/sω= 750 rad/s Mesh Type βError(%) βError(%) A1.2c0.646 1.1c1.696 B2.2c2.160 1.8c2.305 N=2 C0.9c7.646 0.7c5.318 D14.4c14.995 11.8c9.739 (5.9) c0.763 c1.700 A1.0c0.131 1.1c0.437 B2.6c0.367 3.4c0.474 N=4 C1.1c2.113 0.8c1.411 D4.0c4.297 3.2c2.729 (5.9) c0.131 c0.447 A1.0c0.029 1.2c0.101 B2.8c0.070 2.6c0.111 N=8 C1.1c0.589 0.9c0.365 D7.6c0.957 6.8c0.602 (5.9) c0.029 c0.109 To assess the order of convergence of the proposed numerical method, we show in Figure 5.12 the error curves (log-log plots of errors versus mesh-size) for ω= 250 rad/s and ω= 750 rad/s. It can be seen from this figure that an order of convergence O(h2)is achieved. Let us recall that this is the optimal order for the used finite elements in L2-norm.
5.5. Comparison with classical absorbing functions 89 To end this section, we show in Figure 5.13 the real and imaginary parts of the solution computed with the proposed PML method for the mesh corresponding to N= 8 and ω= 750 rad/s. The solution is plotted in the physical domain and in the PML. 0.0156 0.0312 0.0625 0.125 10−3 10−2 10−1 100 101 h (mesh−size) Relative error (%) Relative error, ω=750 rad/s Relative error, ω=250 rad/s Order h2 Figure 5.12: Error curves of the PML method with absorbing functions (5.9). 5.5 Comparison with classical absorbing functions The aim of this section is to compare the proposed unbounded absorbing function (5.9) with the most standard classical choice: quadratic functions of the form σ1(x1)=σ∗(x1−a)2and σ2(x2)=σ∗(x2−b)2,(5.10) where σ∗is a parameter to be determined. For the comparison we have used the same numerical test as in the previous section. When these quadratic absorbing functions are used, the standard procedure to minimize the spurious reflections produced at the outer boundary of the PML consists of taking large values for σ∗. Notice that this agrees with the analysis in Section 4.4. However, larger values of σ∗lead to larger discretization errors. Therefore, σ∗cannot be chosen arbitrarily large because, otherwise, the discretization errors would be dominant, deteriorating the overall accuracy of the method. As shown in [51], for a given problem and a given mesh there is an optimal value of σ∗leading to minimal errors. Unfortunately, such optimal value depends strongly on the problem data as well as on the particular mesh. Thus, in practice, it is necessary to find in advance a reasonable value of σ∗. No theoretical procedure to tune this parameter is known
90 Chapter 5. An optimal PML in Cartesian coordinates Figure 5.13: Solution of the scattering problem computed by the PML method with absorbing functions (5.9). Mesh N=8,ω= 750 rad/s. to date. Some efforts have been done in [100], but the dependency of σ∗with respect to the mesh has not been avoided. Let us emphasize that a benefit of our proposed PML strategy is that it does not need of any parameter to be fitted. In Table 5.2 we compare the errors of the PML method with the unbounded absorbing functions (5.9) and with the quadratic absorbing functions (5.10). For the latter, we have used the optimal value of σ∗, which is also reported in the table. We also include in the table the condition number κof the system matrix for each discrete problem. A significant advantage of the proposed unbounded absorbing functions (5.9) can be clearly appreciated from this table. This is particularly remarkable for lowest frequencies, but the errors with the quadratic absorbing functions are larger in all cases, even though the optimal value of σ∗has been used. On the other hand, in spite of the singular character of the unbounded functions, the condition numbers of the resulting system matrices are essentially of the same order as those of the quadratic functions. On the other hand, Table 5.2 also shows that the optimal value of σ∗strongly depends on the problem data (the frequency ωin this case) and the mesh. The errors and the condition numbers would be significantly larger if any other value than the optimal σ∗were
5.5. Comparison with classical absorbing functions 91 Table 5.2: Comparison of PML methods with unbounded and quadratic absorbing functions. Unbounded (5.9) Quadratic (5.10) ω(rad/s) Mesh Error(%) κσ ∗Error(%) κ N= 2 0.763 6.7e+02 22.28 c11.644 4.7e+02 250 N= 4 0.131 5.1e+03 29.57 c3.675 5.0e+03 N= 8 0.029 4.1e+04 38.37 c1.134 4.6e+04 N= 2 1.700 1.1e+02 27.67 c7.602 1.1e+02 750 N= 4 0.447 7.0e+02 35.52 c2.291 9.4e+02 N= 8 0.109 5.6e+03 43.49 c0.698 8.2e+03 N= 2 6.958 2.7e+02 27.89 c11.620 2.9e+02 1250 N= 4 1.946 1.1e+03 36.94 c3.336 1.7e+03 N= 8 0.430 9.7e+03 45.70 c0.919 1.5e+03 used. This can be appreciated from Figure 5.14 and 5.15, where the relative error and the condition number are respectively plotted as functions of σ∗, for the mesh corresponding to N= 4 and ω= 750 rad/s. 0 50 100 150 20 0 0 10 20 30 40 Error for quadratic PML σ*/c Relative error Figure 5.14: Relative error for quadratic absorbing functions with varying σ∗. Mesh: N=4;ω= 750 rad/s. 0 50 100 150 20 0 0 1e03 2e03 3e03 4e03 5e03 Condition number for quadratic PML σ*/c Condition number Figure 5.15: Condition number for quadratic absorbing functions with varying σ∗. Mesh: N=4;ω= 750 rad/s. As a conclusion, the proposed PML method with unbounded absorbing function (5.9) clearly beats the classical choice of bounded absorbing functions. Moreover, it overcomes the problem of determining optimal parameters.
98 Chapter 5. An optimal PML in Cartesian coordinates Figure 5.24: Real part of the pressure field generated by an incident plane wave, k=2πm−1. Figure 5.25: Imaginary part of the pressure field generated by an incident plane wave, k=2πm−1. Figure 5.26: Real part of the pressure field generated by an incident plane wave, k=10πm−1. Figure 5.27: Imaginary part of the pressure field generated by an incident plane wave, k=10πm−1. The computation of the integrals involving γ2depend on the location of the element K. If −b≤b1<b 2≤b, then γ2= 1 and the integrals are trivial. If b≤b1<b 2≤b∗, then b2 b1 dx2 γ2(x2)=b2 b1 ω(b∗−x2) ω(b∗−x2)+ic dx2=h2+ic ωlog ω(b∗−b2)+ic ω(b∗−b1)+ic and b2 b1 γ2(x2)(x2−bj)(x2−bi)dx2=b2 b1 (x2−bj)(x2−bi)dx2+iω cb2 b1 (x2−bj)(x2−bi) b∗−x2 dx2. The first integral above is trivial, whereas for the second one straightforward computations
5.8. Computation of the element matrices 99 Figure 5.28: Real part of pressure field generated by a monopole, k= 2πm−1. Figure 5.29: Imaginary part of the pressure field generated by a monopole, k=2πm−1. Figure 5.30: Real part of the pressure field generated by a monopole, k=10πm−1. Figure 5.31: Imaginary part of the pressure field generated by a monopole, k=10πm−1. lead to b2 b1 (x2−bj)(x2−bi) b∗−x2 dx2=h2 2(2b∗−b2−b1)−h2(2b∗−bi−bj) +(b∗−bj)(b∗−bi) log b∗−b1 b∗−b2. Similar results are valid if −b∗≤b1<b 2≤−b. Alternatively, all these integrals can be computed using standard quadrature rules. In principle, these rules could lead to large truncation errors due to the singular character of the unbounded absorbing functions. However, our preliminary experiments show that the effect of numerical quadrature does not seem to modify significantly the accuracy of the proposed PML method. To show this we have solved the numerical test from Section 5.4 with the integrals
100 Chapter 5. An optimal PML in Cartesian coordinates 0 .1 m 0.5 m 0.25 m 2 m Figure 5.32: Domain and coarse mesh for the double-slit interference test. Figure 5.33: Wave field generated by the excitation of the waveguides. computed by Gauss-Legendre rules with 4 and 9 nodes. We report in Table 5.3 the relative errors of the solutions computed with each rule and with exact integration. It can be clearly seen that the errors in the numerical integration are negligible. Moreover, for 250 rad/s and 750 rad/s the 4-nodes rule shows a slightly better performance,
5.8. Computation of the element matrices 101 −ka 0 ka −0.4 −0.2 0 0.2 0.4 θ Real part of wave field r=2.0 r=1.5 r=1.0 Figure 5.34: Interference pattern at different distances rfrom the waveguides apertures. Table 5.3: Comparison of quadrature rules using PML with unbounded absorbing functions. Gauss-Legendre Exact Integration ω(rad/s) Mesh 4 nodes 9 nodes N= 2 0.763689 0.770405 0.763485 250 N= 4 0.130572 0.130413 0.130580 N= 8 0.028858 0.028755 0.028860 N= 2 1.699869 1.699611 1.699889 750 N= 4 0.446922 0.446910 0.446922 N= 8 0.109444 0.109467 0.109443 N= 2 6.957597 6.958152 6.958012 1250 N= 4 1.946417 1.946320 1.946313 N= 8 0.429920 0.429913 0.429912 which agrees with the fact that low-order integration schemes are preferable for singular integrands.
102 Chapter 5. An optimal PML in Cartesian coordinates
Chapter 6 An exact bounded PML in radial coordinates Contents 6.1 Introduction ..............................104 6.2 Scattering problem ..........................105 6.3 Statement of the PML equation ...................106 6.4 PML fundamental solution ......................108 6.5 PML integral representation formula ...............115 6.6 Addition theorem ...........................118 6.7 Existence and uniqueness of solutions for the PML equation . 120 6.8 Coupled fluid/PML problem .....................122 6.9 Discretization and numerical results ................125 6.10 Appendix ................................128 6.10.1 Technical results ........................... 128 6.10.2 Some classical results about the Hankel functions ......... 135 103
104 Chapter 6. An exact bounded PML in radial coordinates 6.1 Introduction We have shown in Chapter 5 that the Cartesian PML technique based on singular absorbing functions, leads to accurate numerical results. The aim of this chapter is to analyze mathematically the existence and uniqueness of solution of the corresponding coupled fluid/PML problem. We also prove that this choice leads to a theoretically exact bounded PML. More precisely, this kind of absorbing function on a circular annular layer allows recovering the exact solution of the time-harmonic scattering problem in the domain of interest, up to discretization errors, even though the thickness of the layer is finite. This is the reason why we call “exact” PML methods to those based on such absorbing functions. Standard PML techniques based on bounded absorbing function lead to partial differential equations in the PML with bounded coefficients. Thus, the theoretical procedure to prove the well-posedness of the coupled fluid/PML problems is based on the Freedholm alternative in standard Sobolev spaces. However, when a non-integrable absorbing function is used, the coefficients in the PML equation become unbounded, and the natural functional framework involves a non-standard weighted Sobolev space. In this case standard arguments cannot be applied due to the lack of a compactness result. As an alternative, we reproduce the classical steps used for the Helmholtz equation, taking advantage of the series representation of the solution in the PML domain. Thus we prove a result of existence and uniqueness for the coupled fluid/PML problem. The analysis of the theoretical error for other PML techniques is typically based on the construction of an analogous Dirichlet-to-Neumann operator using the solution in the PML. We also use this approach to prove that the solution in the fluid domain of the coupled fluid/PML problem is exactly equal to the solution of the scattering problem in an unbounded domain. The outline of this chapter is as follows: we recall in Section 6.2 how the classical timeharmonic scattering problem can be stated in a bounded domain by using a DtN operator. Section 6.3 is devoted to settling PML equations based on non-integrable absorbing functions on an annular domain surrounding the physical one. Once the fundamental solution for the PML is calculated in Section 6.4, a integral representation formula is proved in Section 6.5. We rewrite a classical “addition theorem” for the PML fundamental solution in Section 6.6. Using these tools, we prove a characterization theorem for the solution of the radial PML in Section 6.7 and derive a theorem of existence and uniqueness of solution for the PML problem. In Section 6.8, we use this result to recover the classical solution of the scattering problem by means of a coupled fluid/PML problem. We prove existence and uniqueness of solution for this coupled problem and write a weak formulation, as well. Finally, in Section 6.9, we report some numerical results obtained with a standard finite element method.
6.2. Scattering problem 105 6.2 Scattering problem Let Ω be a bounded two-dimensional open set with a Lipschitz boundary Γ. We aim to solve a scattering problem in the unbounded domain R2\Ω, which we assume connected (see Figure 6.1). As in Chapter 5, we are going to focus our attention on the time-harmonic scattering problem. More precisely, we consider the following Dirichlet boundary value problem for the Helmholtz equation, which models the wave propagation with frequency ω>0 and velocity of propagation c>0: Find p∈H1 loc(R2\Ω) such that −Δp−k2p=0 inR2\Ω,(6.1) p=fon Γ,(6.2) lim r→∞ √r∂p ∂r −ikp=0,(6.3) where r:= |x|is the radial polar coordinate for x∈R2,k:= ω/c is the wave number, and f∈H1 2(Γ) is the Dirichlet boundary data. Let us remark that we could analogously consider the corresponding Neumann boundary value problem. The existence of solution to both problems is well known in the literature (see for instance [104]). Ω ΩR SR ν R2\BR Γ Figure 6.1: Scatterer and artificial circular boundary. Let BR:= {x∈R2:|x|<R}be an open ball of radius Rsuch that Ω ⊂BR. Let SR:= {x∈R2:|x|=R}be its boundary and νits outward unit normal vector (see Figure 6.1). The DtN operator of the problem above is defined as follows: G:H 1 2(SR)−→ H−1 2(SR) g−→ ∂˜p ∂νSR (6.4)
106 Chapter 6. An exact bounded PML in radial coordinates where ˜p∈H1 loc(R2\BR) is the unique solution of −Δ˜p−k2˜p=0 inR2\BR, ˜p=gon SR, lim r→∞ √r∂˜p ∂r −ik˜p=0. Let us recall that this operator is explicitly given by the following series (see [84, 90]): Gg = ∞ n=−∞ gnk[H(1) n](kR) H(1) n(kR)einθ, where θis the angular polar coordinate of x,gn:= 1/(2πR)!SRg(x)e −inθ dS is the nth Fourier coefficient of g, and H(j) ndenotes the nth Hankel function of jth kind, j=1,2 (see for instance [103]). Problem (6.1)–(6.3) can be equivalently settled in the bounded domain ΩR:= x∈R2\Ω:|x|<R % by means of this DtN operator as follows: Find p∈H1(ΩR)such that −Δp−k2p= 0 in ΩR,(6.5) p=fon Γ,(6.6) ∂p ∂ν=Gp|SRon SR.(6.7) Clearly, if pis the solution of Problem (6.1)–(6.3) (see Figure 6.1), then p|ΩRis the unique solution of the problem above. 6.3 Statement of the PML equation Radial PML methods are based on simulating dissipation in an annular domain, D := {x∈R2:R<|x|<R }, surrounding the physical domain of interest (see Figure 6.2). This can be done by means of a complex-valued radial stretching proposed by Collino and Monk [52], which leads to the following partial differential equation written in polar coordinates: −1 r∂ ∂r ˆγ(r)r γ(r) ∂ˆp ∂r+γ(r) ˆγ(r)r ∂2ˆp ∂θ2−γ(r)ˆγ(r)k2ˆp= 0 in D,(6.8) where γ(r):=1+ i ωσ(r) and ˆγ(r):=1+ i ωr r R σ(s),(6.9)
6.3. Statement of the PML equation 107 R R Figure 6.2: Domain D. with the so-called absorbing function σ:[R, R)−→ [0,∞) being monotonically increasing and smooth. Although in practice σwill be chosen in C∞([R, R)), in all what follows it is enough to assume that σ∈C 2,1([R, R]) for all R∈(R, R) (recall that C2,1([R, R]) := {F∈C2([R, R]); F is Lipschitz-continuous in [R, R]}). Notice that we do not assume that σ(R) = 0. This function has been typically chosen bounded (see [52]). As an alternative, we propose to choose a non-integrable function, i.e., such that R R σ(s)=+∞; (6.10) for example, σ(r):=c/(R−r). Under the previous assumptions on σ, limr→R|γ(r)|= limr→R|ˆγ(r)|=+∞. Moreover, the coefficients of the differential equation (6.8) satisfy lim r→R ˆγ(r)r γ(r)= 0 and lim r→R γ(r) ˆγ(r)r=+∞.(6.11) Indeed, given >0, let A:= !R− Rσ(s). Because of (6.10), ∃r∈(R−, R) such that !r R−σ(s)≥A. Then r R σ(s)=R− R σ(s)+r R− σ(s)≤2r R− σ(s)≤2σ(r). Hence, lim r→R 1 σ(r)r R σ(s)=0,
114 Chapter 6. An exact bounded PML in radial coordinates where we have used (6.23) and Lemma 6.10.3 (which is proved in Section 6.10 below). Thus we conclude the theorem. To finish this section, we settle some decay properties of the fundamental solution Φ+. Lemma 6.4.7. For fixed x∈D, there holds uniformly in θy∈(−π,π]: lim ry→RˆγyγyΦ+(x,y)=0, lim ry→Rˆγy γy ∂Φ+(x,y) ∂ry =0, lim ry→Rγy ˆγy ∂Φ+(x,y) ∂θy =0. Proof. Let x∈D fixed. Using the asymptotic classical estimates for Hankel functions (see [103]) and Lemma 6.4.5, we can check that Φ+(x,·) and their derivatives go to zero exponentially and uniformly in all the directions, as ry→R. In fact, these three limits are satisfied for each Hankel function of first kind and order n. From item i) of Lemma 6.4.5, Im(d(x,y)) tends to infinity uniformly in all directions θy∈(−π,π] when rygoes to R. Hence, using the asymptotic behavior (6.82) and (6.83), we can check that the first kind Hankel functions and their derivatives go to zero exponentially while those of the second kind (and their derivatives) increase exponentially, in both cases uniformly in all the directions. Indeed, using the asymptotic behavior (6.82) and the limits (6.17)-(6.19), the outcoming fundamental solution satisfies, lim ry→RˆγyγyH(1) n(kd(x,y)) = lim ry→R2ˆγyγy πkd(x,y)ei(kd(x,y)−nπ 2−π 4)1+O1 |d(x,y)| =2 πkRlim ry→R√γyeikˆry1+O1 |d(x,y)|=0,(6.25) lim ry→Rˆγy γy ∂H(1) n(kd(x,y)) ∂ry = lim ry→Rˆγy γy ∂d(x,y) ∂ry ik2 πkd(x,y)ei(kd(x,y)−nπ 2−π 4)1+O1 |d(x,y)| =2k πRlim ry→R√γyeikˆry1+O1 |d(x,y)|=0,(6.26)
6.5. PML integral representation formula 115 and lim ry→Rγy ˆγy ∂H(1) n(kd(x,y)) ∂θy = lim ry→Rγy ˆγy ∂d(x,y) ∂θy ik2 πkd(x,y)ei(kd(x,y)−nπ 2−π 4)1+O1 |d(x,y)| =2k πRlim ry→R√γy ˆγy eikˆry1+O1 |d(x,y)|=0,(6.27) In the previous limits we have taken into account that lim ry→R√γyeikˆry= 0. In fact, lim ry→R√γyeikˆry=1 2ik lim ry→R ∂ ∂ry (e 2ikˆry)1 2 =1 2ik lim ry→R e2ikˆry ry−R1 2 =0, because of the definition of ˆr. Moreover, all the previous limits are uniformly in all the directions θy∈(−π,π] since θyis not involved in the definition of ˆry. 6.5 PML integral representation formula The aim of this section is to obtain an integral representation of the solutions of the PML equation (6.8). We search for smooth solutions which, furthermore, belong to the functional space V:=q∈H1 loc(D) : q2 V:= R Rπ −π ˆγ(r)r γ(r) ∂q ∂r 2 dθ dr +R Rπ −π γ(r) ˆγ(r)r ∂q ∂θ 2 dθ dr +R Rπ −π|ˆγ(r)γ(r)r||q|2dθ dr < +∞,. As a first step, we restrict our analysis to solutions of (6.8) in the space W:=V∩C1(D)∩C2(D), where D:= D ∪SR={x∈R2:R≤|x|<R }. Since the weights involved in the definition of V belong to L1 loc(R, R) and are positive, V is a Banach space when endowed with the norm ·V(see Kufner & S¨andig [75]) and, moreover, V ⊂H1 loc(D), so that q∈H1(K) even for compact sets Kintersecting SR. First, we prove two preliminary results.
116 Chapter 6. An exact bounded PML in radial coordinates Lemma 6.5.1. If q∈V, then lim ˜ R→RS˜ R|ˆγ||q|2dS =0. Proof. For q∈V, we define the complex-valued function Fgiven by F(r):=Sr ˆγ|q|2dS =π −π rˆγ(r)|q(r, θ)|2dθ. From the definition of V, it is immediate to check that Fand γF belong to L1(R, R). We also define G(r):=π −π ∂ ∂r rˆγ(r)|q(r, θ)|2dθ =γ(r)π −π|q(r, θ)|2dθ +2rˆγ(r)π −π Re ∂q ∂r(r, θ)¯q(r, θ)dθ. For q∈V, the first term in the above sum is integrable in (R, R). Regarding the second term, we have R R2rˆγ(r)π −π Re ∂q ∂r(r, θ)¯q(r, θ)dθdr ≤2'R Rπ −π ˆγ(r)r γ(r) ∂q ∂r 2 dθ dr(1 2R Rπ −π|ˆγ(r)γ(r)r||q|2dθ dr1 2 , which is again finite for q∈V. Thus G∈L1(R, R). Moreover, straightforward computations show that Gis the distributional derivative of F. Hence F∈W1,1(R, R) and, consequently, F∈C([R, R]) (see for instance Theorem VIII.2 in [44]). Now we can conclude the lemma by showing only that limr→RF(r) = 0, since Re(γ) and Im(γ) are non negative. We proceed by contradiction. Suppose limr→RF(r)=0;in such a case, since |γ|>1 and γis not integrable near r=R,!R R|γF|dr =∞, which would contradict the fact that γF ∈L1(R, R). Lemma 6.5.2. If ˆp∈Wis a solution of (6.8) and q∈V, then lim ˜ R→RS˜ R ˆγ γ ∂ˆp ∂rqdS=0. Proof. Let ˜ R∈(R, R). Since q∈V⊂H1 loc(D), if we multiply (6.8) by q∈V and integrate by parts in ˜ D:={x∈R2:R<|x|<˜ R}, we obtain S˜ R ˆγ γ ∂ˆp ∂rqdS=˜ D ˆγ γ ∂ˆp ∂r ∂q ∂r +˜ D γ ˆγr2 ∂ˆp ∂θ ∂q ∂θ −k2˜ D γˆγˆpq +SR 1 γ ∂ˆp ∂rqdS =˜ R RSrˆγ γ ∂ˆp ∂r ∂q ∂r +γ ˆγr2 ∂ˆp ∂θ ∂q ∂θ −k2γˆγˆpqdS dr +SR 1 γ ∂ˆp ∂rqdS. (6.28)
6.5. PML integral representation formula 117 Because of the definition of V, the expression between brackets above belongs to L1(R, R). Consequently, if we define H(˜ R):=S˜ R ˆγ γ ∂ˆp ∂rqdS, then, according to (6.28), from [44] again we have H∈C([R, R]). On the other hand, from the Cauchy-Schwarz inequality, we have R R|γ(r)H(r)|dr =R Rπ −π ˆγ(r)r∂ˆp ∂rqdθ dr ≤'R Rπ −π ˆγ(r)r γ(r) ∂ˆp ∂r 2 dθ dr(1 2R Rπ −π|ˆγ(r)γ(r)r||q|2dθ dr1 2 . which is finite for ˆp, q ∈V. Thus, γH ∈L1(R, R) and, hence, the same argument used to prove the previous lemma, allow us to conclude that limr→RH(r)=0. Now, the following step is to establish an integral representation of the solution of the PML equation. Theorem 6.5.3. If ˆp∈Wis a solution of (6.8), then the following integral representation formula holds true: ˆp(x)= 1 γ(R)SR∂Φ+(x,y) ∂ry ˆp(y)−∂ˆp ∂ry (y)Φ+(x,y)dSy,x∈D.(6.29) Proof. We fix an arbitrary x∈D and use the notation from Figure 6.3. As shown in the proof of Theorem 6.4.6, Φ+(x,·) satisfies (6.23). Hence, since Ayis diagonal, by using the Green’s second theorem and (6.8), we have ∂˜ DAygradyˆp(y)·nΦ+(x,y)−AygradyΦ+(x,y)·nˆp(y)dSy(6.30) =˜ Ddiv(Aygradyˆp(y))Φ+(x,y)−div(AygradyΦ+(x,y))ˆp(y)dy=0, where nis the outward unit normal vector to ˜ D. By using Lemma 6.10.3 below (see Appendix) and (6.30), we obtain ˆp(x) = lim →0S(x,) AygradyΦ+(x,y)·nˆp(y)dSy −S(x,) Aygradyˆp(y)·nΦ+(x,y)dSy =1 γ(R)SR∂Φ+(x,y) ∂ry ˆp(y)−∂ˆp ∂ry (y)Φ +(x,y)dSy −S˜ R ˆγy γy ∂Φ+(x,y) ∂ry ˆp(y)dSy+S˜ R ˆγy γy ∂ˆp ∂ry (y)Φ +(x,y)dSy.
118 Chapter 6. An exact bounded PML in radial coordinates To conclude the proof, it is enough to show that the last two integrals go to zero as ˜ R→R. For the first one we use Lemmas 6.4.7 and 6.5.1. For the second one, first we replace Φ+(x,y)byζ(y)Φ+(x,y) where ζis a smooth cutoff function vanishing in a neighborhood of xand taking the value 1 in {x∈R2:˜ R≤|x|<R }. Thus, the value of the integral does not change and ζ(y)Φ+(x,y)∈V, because of Lemma 6.4.7. Hence, the integral goes to zero as ˜ R→Ras a consequence of Lemma 6.5.2. 6.6 Addition theorem To characterize the solution of the PML equation, it is useful to write the fundamental solution as a series involving Bessel functions of first kind and order n, which we denote as usual by Jn. Theorem 6.6.1. Let x∈Dbe fixed. For all y∈Dsuch that |y|<|x|, there holds: Φ+(x,y)= i 4 ∞ n=−∞ H(1) n(kˆrx)J n(kˆry)e in(θx−θy).(6.31) This series and its term by term first derivatives with respect to ryare absolutely and uniformly convergent on compact subsets of the set {y∈R2:R≤|y|<|x|}. Proof. Let ˜ R∈[R, |x|). First, we define the following functions: φ(y):=⎧ ⎨ ⎩ i 4H(1) 0(k|x−y|),if 0 ≤|y|<R, Φ+(x,y),if R≤|y|<R , and, for each n∈N, pn(y):=Jn(kry)e inθy,if 0 ≤|y|<R, Jn(kˆry)e inθy,if R≤|y|≤˜ R. All these functions are continuous for |y|<˜ R, analytic for |y|<Rand C2for R≤|y|≤˜ R. Moreover, they satisfy lim r→R+ 1 γ ∂φ ∂r = lim r→R− ∂φ ∂r and lim r→R+ 1 γ ∂pn ∂r = lim r→R− ∂pn ∂r ,n∈N. Furthermore, straightforward computations allow us to show that φand pn,n∈N, all are solutions of the PML equation (6.8) in R≤|y|<˜ Rand solutions of the Helmholtz equation in 0 <|y|<R. Next, we proceed as in the proof of the addition theorem for the Helmholtz equation (see [54]), taking care of the fact that the Helmholtz equation is substituted by the PML equation for R≤|y|<˜ R. Thus we obtain S˜ R ˆγ γpn ∂φ ∂r −∂pn ∂r φdS =0.(6.32)
6.6. Addition theorem 119 On the other hand, straightforward computations allow us to show that, for all n∈N, if we define qn(y):=H (1) n(kˆry)e inθy,y∈D,(6.33) then qnare solutions of the PML equation (6.8) and belong to W. By applying the analogous of Theorem 6.5.3 for qninstead of ˆp, on the annular domain ˜ R<|y|<R instead of D, and taking into account that φ=Φ +(x,·) in this domain, we obtain S˜ R ˆγ γqn ∂φ ∂r −∂qn ∂r φdS =qn(x).(6.34) Next, multiplying equation (6.32) by H(1) n(kˆr(˜ R)) and (6.34) by Jn(kˆr(˜ R)), we have S˜ R ˆγ γH(1) n(kˆr(˜ R))pn ∂φ ∂r −H(1) n(kˆr(˜ R))∂pn ∂r φdS =0, S˜ R ˆγ γJn(kˆr(˜ R))qn ∂φ ∂r −Jn(kˆr(˜ R))∂qn ∂r φdS =J n(kˆr(˜ R))qn(x). If we subtract the first from the second equation, taking into account that H(1) n(kˆr(˜ R))pn= Jn(kˆr(˜ R))qnon S˜ R, we obtain Jn(kˆr(˜ R))qn(x)=H(1) n(kˆr(˜ R))kγ(˜ R)J n(kˆr(˜ R)) −Jn(kˆr(˜ R))kγ(˜ R)[H(1) n](kˆr(˜ R))ˆγ(˜ R)˜ R γ(˜ R)π −π φeinθ dθ =−2i ππ −π φeinθ dθ, (6.35) where we have used the explicit value of the Wronskian H(1) n(z)J n(z)−Jn(z)[H(1) n](z)= −2i/(πz) (see [2]). Since φ∈C(D), φ|S˜ Radmits a Fourier series, i.e., φ(y)= ∞ n=−∞ φne−inθy,y∈S˜ R,(6.36) where, from (6.35) and (6.33), φn:= 1 2ππ −π φeinθ dθ =i 4H(1) n(kˆrx)J n(kˆr(˜ R)) einθx.(6.37) Finally, we conclude (6.31) from (6.36) and (6.37), since ˆr(˜ R)=ˆry, for y∈S˜ R. The uniform convergence of the series (6.31) and its term by term first derivatives on compact subsets of {y∈R2:R≤|y|<|x|} is straightforward from the uniform convergence of the analogous series in the addition theorem for the fundamental solution of the Helmholtz equation (see [54]), and the fact that |ˆr(r)|is a monotonically increasing function.
120 Chapter 6. An exact bounded PML in radial coordinates 6.7 Existence and uniqueness of solutions for the PML equation Now we are able to characterize the smooth solutions of the PML equation (6.8): Theorem 6.7.1. If ˆp∈Wis a solution of (6.8), then there exists a sequence {an}such that, for all x∈D, ˆp(x)= ∞ n=−∞ anH(1) n(kˆrx)e inθx. Proof. Let ˆp∈W be a solution of (6.8). For fixed x∈D, if we apply the PML integral representation formula (6.29) and Theorem 6.6.1, then we have ˆp(x)= 1 γ(R)SR∂Φ+(x,y) ∂ry ˆp(y)−∂ˆp ∂ry (y)Φ +(x,y)dSy =1 γ(R)SR'i 4 ∞ n=−∞ kγ(R)H (1) n(kˆrx)J n(kR)e in(θx−θy)ˆp(y) −∂ˆp ∂ry (y)i 4 ∞ n=−∞ H(1) n(kˆrx)J n(kR)e in(θx−θy)(dSy = ∞ n=−∞ anH(1) n(kˆrx)e inθx, where an=i 4 1 γ(R)SRkγ(R)J n(kR)e −inθyˆp(y)−∂ˆp ∂ry (y)J n(kR)e −inθydSy. Now, we prove the existence and uniqueness of smooth solutions of the following problem for the PML equation with Dirichlet data g: Find ˆp∈Wsuch that −div(Agrad ˆp)−γˆγk2ˆp= 0 in D,(6.38) ˆp=gon SR.(6.39) Theorem 6.7.2. If g∈Hs(SR)with s>3/2, then there exists a unique solution ˆp∈Wof (6.38)-(6.39). Moreover, this solution is given by ˆp(x)= ∞ n=−∞ gn H(1) n(kR)H(1) n(kˆrx)e inθx,(6.40) where gn:= 1/(2πR)!SRg(x)e −inθxdS are the Fourier coefficients of g. Moreover the series and its term by term first derivatives converge uniformly on compact subsets of D.
6.7. Existence and uniqueness of solutions for the PML equation 121 Proof. First, we are going to prove that ˆpas defined by (6.40) belongs to W. We split the proof into three steps. The first one consists in proving that ˆp∈C 2(D). This step is essentially identical to what is known for the Helmholtz problem. In fact, taking into account classical estimates of the Hankel functions of first kind for large order (see [103]), for xin any compact subset of D and Nlarge enough, there holds ⎛ ⎝ |n|≥N gn H(1) n(kR)H(1) n(kˆrx)e inθx⎞ ⎠ 2 ≤ |n|≥N H(1) n(kˆrx) H(1) n(kR) 2 |n|≥N|gn|2 ≤Cg2 L2(SR) |n|≥NR |ˆrx|2|n| , where we recall that |ˆrx|>R(here and thereafter Cdenotes a generic constant, not necessarily the same at each occurrence). From this, we conclude the uniform and absolute convergence of the series (6.40) on compact subsets of D. Analogous procedures allow us to prove the uniform convergence of the corresponding series for the first and the second derivatives. Since each term in each series is continuous, we conclude that ˆp∈C2(D). The second step consists in proving that ˆp∈C(D). Since g∈Hs(SR) with s>1/2, for xin any compact subset of Dwe have for Nlarge enough ⎛ ⎝ |n|≥N gn H(1) n(kR)H(1) n(kˆrx)e inθx⎞ ⎠ 2 ≤ |n|≥N 1 ns H(1) n(kˆrx) H(1) n(kR) 2 |n|≥N n2s|gn|2 ≤C |n|≥N 1 ns 2 |n|≥N n2s|gn|2≤Cg2 Hs(SR) |n|≥N 1 n2s, the latter because of the decay behavior of the Fourier coefficients of functions in Hs(SR) (see [81]). This allows us to conclude that ˆp∈C(D). The same arguments as above applied to the term by term derivatives of the series allow us to show that, for g∈Hs(SR) with s>3/2, ∂ˆp/∂r and ∂ˆp/∂θ belong to C(D), too. From the previous steps, clearly ˆp∈H1 loc(D). Hence, since the weights in the norm of V are positive bounded functions in compact subsets of D, in order to prove that ˆp∈V we only need to prove that there exists δ>0 such that R R−δπ −π' ˆγ(r)r γ(r) ∂ˆp ∂r 2 + γ(r) ˆγ(r)r ∂ˆp ∂θ 2 +|ˆγ(r)γ(r)r||ˆp|2(dθ dr < +∞.(6.41) Since Im(ˆrx)→+∞,asrx→R, using standard estimates for H(1) 0and [H(1) 0], and for H(1) nand [H(1) n],n∈N, the latter uniform in n(see [46]), it is straightforward to prove that
122 Chapter 6. An exact bounded PML in radial coordinates the following limits hold uniformly in θx∈(−π,π] and n∈N: lim rx→Rˆγxγx H(1) n(kˆrx) H(1) n(kR)=0, lim rx→R 1 nˆγx γx ∂ ∂rx'H(1) n(kˆrx) H(1) n(kR)(=0, lim rx→R 1 nγx ˆγx H(1) n(kˆrx) H(1) n(kR)=0. Hence, for δsmall enough, we have R R−δπ −π|ˆγ(r)γ(r)r||ˆp|2dθ dr ≤Cδπ −π ∞ n=−∞ gneinθ 2 dθ =Cδg2 L2(SR), R R−δπ −π ˆγ(r)r γ(r) ∂ˆp ∂r 2 dθ dr ≤Cδπ −π ∞ n=−∞ ngneinθ 2 dθ ≤Cδg2 H1(SR), R R−δπ −π γ(r) ˆγ(r)r ∂ˆp ∂θ 2 dθ dr ≤Cδπ −π ∞ n=−∞ ngneinθ 2 dθ ≤Cδg2 H1(SR), which allow us to conclude that the integral (6.41) is finite. Therefore, from the three previous steps we deduce that ˆp∈W. Next, since ˆpis continuous on SR, evaluating (6.40) for x∈SRwe have ˆp(x)= ∞ n=−∞ gneinθ, so that (6.39) follows from the convergence of the Fourier series of g. On the other hand, straightforward computations allow us to show that each term in the series (6.40) is a solution of (6.38). Thus, ˆpis a solution too, because we have already shown the uniform convergence on compact subsets of D of that series and its term by term first and second derivatives. Finally, ˆpis the unique solution of (6.38)-(6.39) in W because of Theorem 6.7.1 and the uniqueness of the Fourier expansion of g. 6.8 Coupled fluid/PML problem Our next goal is to study the coupled fluid/PML problem and to prove that the solution of the classical scattering problem is recovered when the PML is used. Theorem 6.7.2 allows us to define a “Dirichlet-to-Neumann” PML operator, ˆ G:H 1 2(SR)→H−1 2(SR). First we define Gfor sufficiently smooth data as follows: ˆ G(g)= 1 γ(R) ∂ˆp ∂rSR ,g∈Hs(SR), with s>3/2,(6.42)
6.8. Coupled fluid/PML problem 123 where ˆpis the unique solution in W of (6.38)-(6.39). This definition can be extended to g∈H1 2(SR) by means of a density argument, because of the following result: Theorem 6.8.1. There exists a unique bounded linear operator ˆ G:H 1 2(SR)→H−1 2(SR) satisfying (6.42), which coincides with Gas defined by (6.4). Proof. If g∈Hs(SR), s>3/2, and Gis defined by (6.4), then ˆ Gg =Gg. Indeed, ˆ Gg =1 γ(R) ∂ˆp ∂rSR =1 γ(R) ∞ n=−∞ gnk H(1) n(kR) dˆr dr(R)[H(1) n](kˆrx)e inθ = ∞ n=−∞ gnk H(1) n(kR)[H(1) n](kR)e inθ =Gg. Consequently, the definition of ˆ Gextends uniquely to the whole space H1 2(SR) and ˆ G= G. Therefore, ˆ Gcan be equivalently used instead of Gin the definition of problem (6.5)- (6.7). Moreover we have the following result. Theorem 6.8.2. For f∈H1 2(Γ), there exists a unique solution (p, ˆp)∈H1(ΩR)×Vof the following problem: −Δp−k2p=0 in ΩR,(6.43) −div(Agrad ˆp)−γˆγk2ˆp=0 in D,(6.44) p=fon Γ,(6.45) ∂p ∂ν=Agrad ˆp·νin H−1 2(SR),(6.46) p=ˆpon SR.(6.47) Moreover, pcoincides with the solution of (6.5)-(6.7) and, hence, it coincides with the solution of the scattering problem (6.1)-(6.3) in ΩR. Proof. Let p∈H1(ΩR) be the solution of (6.5)-(6.7). Then pis the restriction to ΩRof the solution of (6.1)-(6.3). Hence p|SRis arbitrarily smooth. Thus, ˆ G(p|SR)=(1/γ(R))∂ˆp/∂r, with ˆp∈W being the solution of (6.38)-(6.39). Therefore (p, ˆp)∈H1(ΩR)×V is a solution of the coupled fluid/PML problem. To prove the uniqueness, it is enough to show that the solution (po,ˆpo) of problem (6.43)- (6.47) with f= 0 vanishes. By applying local regularity up to the boundary results for transmission problems (in particular Theorem 4.20 from [85]), we conclude that ˆpo∈C2(D). Notice that this regularity comes from the assumed smoothness on the absorbing function: σ∈C2,1([R, R]) ∀R∈(R, R). Hence ˆpo∈W and Agrad ˆpo·ν=1 γ(R) ∂ˆpo ∂r on SR.
130 Chapter 6. An exact bounded PML in radial coordinates Analogously, if we derive with respect to φy,wehave −ρysin φy=∂ξr(ρy,φ y) ∂φy cos θy−rysin θy ∂ξθ(ρy,φ y) ∂φy , ρycos φy=∂ξr(ρy,φ y) ∂φy sin θy+rycos θy ∂ξθ(ρy,φ y) ∂φy , whose solution is ∂ξr(ρy,φ y) ∂φy =ρy(−sin φycos θy+ cos φysin θy)=−ρysin(φy−θy), ∂ξθ(ρy,φ y) ∂φy =ρy ry (cos θycos φy−sin θysin φy)=ρy ry cos(φy−θy). Summarizing we have ∂ry ∂ρy = cos(φy−θy),∂θy ∂ρy =1 ry sin(φy−θy).(6.55) ∂ry ∂φy =−ρysin(φy−θy),∂θy ∂φy =ρy ry cos(φy−θy).(6.56) Now we are going to study the behavior of the complex distance d(·,·) between close points and near the circumference SR. The following lemma collects several limits that will be used in the proof of Lemma 6.10.3 below. Lemma 6.10.1. For fixed x∈Dand φy∈(−π,π], lim ρy→0 ∂ˆry ∂ρy =γxcos(φy−θx),(6.57) lim ρy→0 ∂d(x,y) ∂ρy ="γ2 xcos2(φy−θx)+ˆγ2 xsin2(φy−θx),(6.58) lim ρy→0 ˆry−ˆrxcos(θy−θx) ry−rxcos(θy−θx)=γx.(6.59) Proof. It is clear that, since d(x,·) and ˆr(|·|) are in C1(D\{x}), we can compute derivatives at any point of D different form x. In order to obtain the first equation, we have ∂ˆry ∂ρy =∂ˆry ∂ry ∂ry ∂ρy .(6.60) From the definition of ˆryit is immediate to check that ∂ˆry ∂ry=γy. Then, from (6.60) and (6.55) we obtain ∂ˆry ∂ρy =γycos(φy−θy).
6.10. Appendix 131 If ρy= 0, then x=yand (6.57) is obtained. For the second equation, firstly we are going to bound the limit. For x∈D and a fixed angle φy∈(−π,π], using the Taylor polinomial of first order for function F(ρy)=d(x,x+ρy(cos φy,sin φy)), which is smooth in [0,+∞), we obtain lim ρy→0 d(x,y) ρy = lim ρy→0 F(ρy) ρy =F(0) = lim ρy→0 ∂d(x,y) ∂ρy .(6.61) Hence, it is clear that F(0) ∈[C1,C 2] where the positive constants C1,C 2are given by lemma 6.4.4. Therefore, F(0) = lim ρy→0 ∂d(x,y) ∂ρy = lim ρy→0∂d(x,y) ∂ˆry ∂ˆry ∂ρy +∂d(x,y) ∂θy ∂θy ∂ρy = lim ρy→0 ˆry−ˆrxcos(θy−θx) d(x,y)γycos(φy−θy) + lim ρy→0 ˆryˆrxsin(θy−θx) d(x,y) sin(φy−θy) ry =γxcos(φy−θx) lim ρy→0 ˆry−ˆrxcos(θy−θx) d(x,y)+ˆrxˆγxsin(φy−θx) lim ρy→0 sin(θy−θx) d(x,y). In the last equality we have used that ˆrx rx=ˆγx. Now, taking into account that F(0) =0, we have F(0) =γxcos(φy−θx) lim ρy→0 ∂ˆry ∂ρy+ˆrx rysin(θy−θx) sin(φy−θy) ∂d(x,y) ∂ρy +ˆrxˆγxsin(φy−θx) lim ρy→0 cos(θy−θx)1 rysin(φy−θy) ∂d(x,y) ∂ρy =γxcos(φy−θx)γxcos(φy−θx) F(0) +ˆrxˆγxsin(φy−θx)1 rx sin(φy−θx) F(0) =F(0)−1γ2 xcos2(φy−θx)+ˆγ2 xsin2(φy−θx). Solving the above equation to obtain F(0), we conclude (6.58). For the last limit (6.59), we can use again an analogous argument involving Taylor polinomial of first order over each factor in the quotient and (6.57), to obtain lim ρy→0 ˆry−ˆrxcos(θy−θx) ry−rxcos(θy−θx)= lim ρy→0 ∂ˆry ∂ρy+ˆrxsin(θy−θx)∂ξθ ∂ρy ∂ry ∂ρy+rxsin(θy−θx)∂ξθ ∂ρy =γx.(6.62) The following integral will be also used below. Lemma 6.10.2. For a∈Cwith Re(a)=0, there holds 1 2ππ −π 1 acos2θ+a−1sin2θdθ = sign(Re(a)).(6.63)
132 Chapter 6. An exact bounded PML in radial coordinates Proof. Firstly it is easy to see that I=π −π 1 acos2θ+a−1sin2θdθ =2aπ 2 −π 2 1 cos2θ 1 a2+ tan2θdθ. Using the change of variable s= tan θ,wehave ds dθ =1 cos2θand then I=2a∞ −∞ 1 a2+s2ds =2 a∞ −∞ 1 1+a−2s2ds = lim R→∞ 2 aR −R 1 1+a−2s2ds. (6.64) In order to evaluate the previous improper integral, we apply the Residues’s theorem. For this purpose we have to calculate the residues of the function f(z)= 1 1+a−2z2= 1 21 1+ia−1z+1 1−ia−1zin the poles ±ia. In fact, Res f(z)=a 2iin z=ia, −a 2iin z=−ia. (6.65) −RR CR ia−1 Figure 6.8: Path in the complex plane (when Im(ia)>0). Now we have two alternatives to be considered: if Im(ia)>0 we obtain R −R 1 1+a−2s2ds =2πiResf(ia)−CR 1 1+a−2s2ds, (6.66) whereas if Im(ia)<0, R −R 1 1+a−2s2ds =2πiResf(−ia)−CR 1 1+a−2s2ds, (6.67)
6.10. Appendix 133 where Cr={reiθ :θ∈[0,π]}with r=Rand r=−R. When Rtends to infinity, the integral on CRgoes to zero and we conclude I= lim R→∞ 2 aR −R 1 1+a−2s2ds = sign(Im(ia))2π= sign(Re(a))2π. (6.68) Now, we are in a position to prove the next lemma, which has been used in Theorems 6.4.6 and 6.5.3. Lemma 6.10.3. For x∈Dfixed, if Φ±are the fundamental solutions given by (6.21) and (6.22), and ϕ∈C1(D), then lim →0S(x,) AygradyΦ±(x,y)·nϕ(y)dSy −S(x,) Aygradyϕ(y)·nΦ±(x,y)dSy=ϕ(x),(6.69) where S(x,)={y∈R2:|x−y|=}and nis its inward unit normal vector. Proof. We prove the lemma for Φ+. An analogous proof is valid for Φ−. First, we check that the limit of the second integral in (6.69) is zero. Since ϕ∈C1(D) and the coefficients of Ayare bounded for all ysuch that |x−y|≤, we only have to prove that lim →0S(x,) Φ+(x,y)dSy=0. The above limit is easy to check by using Lemma 6.4.4 and the following estimate: Φ+(x,y)−i 42i πlog kd(x,y) 2+2Cei π+1 ≤C|d(x,y)|2|log d(x,y)|, for |x−y|small enough, which in its turn follows from the asymptotic behavior of Hankel functions (see [103]). In the above expression, Ceis the Euler’s constant. Regarding the first integral in (6.69), since n=−eρon S(x,), from (6.21), (6.51) and (6.52), we have AygradyΦ+(x,y)·n=ˆγy γy ∂Φ+(x,y) ∂ry er+γy ˆry ∂Φ+(x,y) ∂θy eθ·(−eρ) =−ki 4[H(1) 0](kd(x,y)) ˆγy γy ∂d(x,y) ∂ry cos(φy−θy)+ γy ˆry ∂d(x,y) ∂θy sin(φy−θy) =: −ki 4[H(1) 0](kd(x,y))M(x,y),(6.70)
134 Chapter 6. An exact bounded PML in radial coordinates where M(x,y) denotes the expression between brackets above. By using the following elementary identities (see Figure 6.9): |x−y|cos(φy−θy)=ry−rxcos(θy−θx), |x−y|sin(φy−θy)=rxsin(θy−θx), we obtain M(x,y)=ˆγy ˆry−ˆrxcos(θy−θx) ry−rxcos(θy−θx)cos2(φy−θy)+ γyˆrx rx sin2(φy−θy)|x−y| d(x,y), which, in particular, together with Lemma 6.4.4 and the third limit from Lemma 6.10.1, show that M(x,y) is bounded for |x−y|small enough. θy−θx x O y θy φy |x−y| θy−φy Figure 6.9: Polar coordinates systems centered at the origin Oand at a point x∈R2 On the other hand, by using a classical estimate of [H(1) 0](z) (see [103]) and Lemma 6.4.4, we have [H(1) 0](kd(x,y)) −2i π 1 kd(x,y)≤C|x−y||log(kd(x,y))|,(6.71) for |x−y|small enough. Because of this, we proceed from (6.70) as follows: S(x,) AygradyΦ+(x,y)·nϕ(y)dSy=−π −π k i 4M(x,y)2i π 1 kd(x,y)ϕ(y)dφy −π −π k i 4M(x,y)[H(1) 0](kd(x,y)) −2i π 1 kd(x,y)ϕ(y)dφy.(6.72) Manipulating the second integral from (6.71), we obtain π −π k i 4M(x,y)[H(1) 0](kd(x,y)) −2i π 1 kd(x,y)ϕ(y)dφy ≤Cπ −π k 4|M(x,y)|2|log(kd(x,y))||ϕ(y)|dφy−→ 0as→0,(6.73)
6.10. Appendix 135 since M(x,y) is bounded for |x−y|≤and |d(x,y)|=O() uniformly in all directions (see Lemma 6.4.4). Thus, we only have to calculate the limit of the remaining integral in (6.72). For this purpose, we compute the following limit as ρy=|x−y|→0 for fixed φy∈(−π,π]: lim ρy→0 ρyM(x,y) d(x,y)=⎛ ⎜ ⎜ ⎝ˆγxcos(φy−θx) lim ρy→0 ∂ˆry ∂ρy +ˆrxsin(θy−θx)∂θy ∂ρy ∂d(x,y) ∂ρy +γxˆrxsin(φy−θx) lim ρy→0 cos(θy−θx)∂θy ∂ρy ∂d(x,y) ∂ρy ⎞ ⎟ ⎟ ⎠lim ρy→0∂d(x,y) ∂ρy−1 =γxˆγxlim ρy→0∂d(x,y) ∂ρy−2 =γxˆγx γ2 xcos2(φy−θx)+ˆγ2 xsin2(φy−θx), where we have used L’Hˆopital’s rule, Lemma 6.10.1 and (6.55). Therefore, since ρy=on S(x,), by using the above limit, the boundedness of M(x,y) and Lemma 6.4.4, we have from (6.72) and (6.73) lim →0S(x,) AygradyΦ+(x,y)·nϕ(y)dSy =ϕ(x)1 2ππ −π γxˆγx γ2 xcos2(φy−θx)+ˆγ2 xsin2(φy−θx)dφy = sign Re γx ˆγxϕ(x) =ϕ(x), because of Lemma 6.10.2 with a=γx/ˆγx, which can be shown that has a positive real part. 6.10.2 Some classical results about the Hankel functions In the two dimensional case, the elementary solutions of Helmholtz equation are the Hankel functions of first and second kind of order n, which are denoted by H(1) nand H(2) n, respectively. In this appendix, we focus our attention on some simple asymptotic properties of the Hankel functions. For a more detailed analysis we refer to Lebedev [80] or Watson [103]. Firstly, we can write the Hankel functions using Bessel functions as H(1) n(z):=J n(z)+iY n(z), H(2) n(z):=J n(z)−iY n(z),
136 Chapter 6. An exact bounded PML in radial coordinates where Jnand Ynare the Bessel functions of first and second kind of order n. The following series expansions for the Bessel functions are known: Jn(z)= ∞ k=0 (−1)k k!(n+k)! z 2n+2k,(6.74) Yn(z)=2 πJn(z) log z 2−1 π n−1 k=0 (n−k−1)! k!z 22k−n −1 π ∞ k=0 (−1)k k!(n+k)! z 2n+2k(ψ(k+1)+ψ(k+n+ 1)),(6.75) where ψ(k+1)=−Ce+1+···+1 k,k∈Nand Ceis the Euler’s constant. The sum of the series (6.74) and (6.75) are analytic functions of zin the complex plane cut along (−∞,0]. Lemma 6.10.4. The Hankel functions of order 0 verify the following asymptotic behavior when the modulus of the argument is small, i. e., when |z|→0 H(1) 0(z)=2i πlog z 2 2Cei π+1+O(|z|2|log z|),(6.76) dH(1) 0 dz (z)=2i π 1 z+O(|z||log z|).(6.77) Proof. For n= 0, using (6.74) and (6.75) it is immediate to conclude that H(1) 0(z)=J 0(z)+iY 0(z) =1+O(|z2|)+i2 π1+O(|z|2)log z 2+2Ce π+O(|z|2) =2i πlog z 2+2Cei π+1+O(|z|2|log z|), when |z|→0. On the other hand, taking into account (6.74) and (6.75) for n=1, dH(1) 0 dz (z)=−H(1) 1(z)=−J1(z)−iY 1(z) =O(|z|)−i2 πO(|z|) log z 2−1 π 2 z+O(|z|) =2i π 1 z+O(|z||log z|)+O(|z|)=2i π 1 z+O(|z||log z|), when |z|→0. Analogously, using the definition of H(2) n, (6.74) and (6.75) we obtain the results for the Hankel functions of second kind.
6.10. Appendix 137 Lemma 6.10.5. The following asymptotic behavior with respect to the order nis verified Jn(z)= zn 2nn!1+O1 n,(6.78) H(1) n(z)=2n(n−1)! iπzn1+O1 n,(6.79) when n→∞. Moreover, analogous asymtotic behaviour is obtained for their derivatives: dJn dz (z)= zn−1 2n(n−1)! 1+O1 n,(6.80) dH(1) n dz (z)= −2nn! iπzn+1 1+O1 n,(6.81) when n→∞. Proof. Firstly, from (6.74) we have Jn(z)= ∞ k=0 (−1)k k!(n+k)! z 2n+2k=zn 2nn! ∞ k=0 (−1)kn! k!(n+k)! z 22k =zn 2nn!1+O1 n=O1 n, when n→∞. In a similar way, Yn(z)=−(n−1)! π2 zn1+O1 n+O1 n+O1 n! =−(n−1)! π2 zn1+O1 n, when n→∞. Then, from the definition of H(1) n, we obtain the asymptotic behavior (6.79). Finally, we can obtain estimates (6.80) and (6.81) again taking into account the series expansions of the Bessel and Hankel functions. Remark 6.10.6. We must remark how the expressions (6.78)-(6.81) must be understood. For instance, (6.78) implies that there exists C>0and N∈Nsuch that for all n>N,if n>N then 1−C n≤ Jn(z) zn 2nn! ≤1+C n. We finish this section with a more technical and classical lemma proved, for example, in [103].
138 Chapter 6. An exact bounded PML in radial coordinates Lemma 6.10.7. For large argument, we have the following asymptotic behavior of the Hankel functions H(1) n(z)=2 πz ei(z−nπ 2−π 4)1+O1 |z|,(6.82) dH(1) n dz (z)=±i2 πz ei(z−nπ 2−π 4)1+O1 |z|,(6.83) when |z|→∞and |arg(z)|<π−δ, where δis an arbitrary small positive number.
Part III Computational applications on dissipative acoustics 139