D Lafarge

Keywords: porous media, rigid frame, viscothermal fluid, nonlocal theory, local theory, homogenization, Bloch wave, Kirchhoff equation, Kirchhoff-Langevin equation, temporal dispersion, spatial dispersion, circular tube, rigid cylinders, Helmholtz resonators

This work is dedicated to present and validate a new and generalized macroscopic nonlocal theory of sound propagation in rigid-framed porous media saturated with a viscothermal fluid. This theory allows to go beyond the limits of the classical local theory and within the limits of linear theory, to take not only temporal dispersion, but also spatial dispersion into account. In the framework of the new approach, a homogenization procedure is proposed to upscale the dynamics of sound propagation from Navier-Stokes-Fourier scale to the volume-average scale, through solving two independent microscopic action-response problems. Contrary to the classical method of homogenization, there is no length-constraint to be considered alongside of the development of the new method, thus, there is no frequency limit for the medium effective properties to be valid. In absence of solid matrix, this procedure leads to Kirchhoff-Langevin's dispersion equation for sound propagation in viscothermal fluids.

The new theory and upscaling procedure are validated in three cases corresponding to three different periodic microgeometries of the porous structure. Employing a semi-analytical method in the simple case of cylindrical circular tubes filled with a viscothermal fluid, it is found that the wavenumbers and impedances predicted by nonlocal theory match with those of the long-known Kirchhoff's exact solution, while the results by local theory (Zwikker and Kosten's) yield only the wavenumber of the least attenuated mode, in addition, with a small discrepancy compared to Kirchhoff's.

In the case where the porous medium is made of a 2D square network of cylindrical solid inclusions, the frequency-dependent phase velocities of the least attenuated mode are computed based on the local and nonlocal approaches, by using direct Finite Element numerical simulations. The phase velocity of the least attenuated Bloch wave computed through a completely different quasi-exact multiple scattering method taking into account the viscothermal effects, shows i ii

It has been my extreme good fortune to meet Denis Lafarge and to have him as my supervisor, and at the same time, a very good friend for the past three years. He gave me the opportunity to work on his very own idea which he had been thinking about and also working on himself for the past few years, and it has been an honor for me to have it as my PhD subject and complete the research and enjoy its outcome. In addition, working with Denis on this particular subject has been very rewarding to me considering sustained contemplation of the principles of fluid mechanics and thermodynamics. I here would like to deeply thank him for his tremendous support, and sharing with me a huge part of his own time for discussing on various subjects, in particular physics and music, as well as other subjects whether related or not to the research helping me out with his wide vision.

I would also like to greatly thank Yves Aurégan as my official advisor, although being very busy in academic life as the director of a large and high quality population laboratory, for sharing his time with me whenever required without hesitation and helping me out with challenges with his professional and intelligent academic and administrative advices at the right time.

I am grateful to Claude Boutin and Mario Silveirinha for having reviewed precisely this thesis and for their interests in this work. Many thanks also to Julius Kaplunov, Vincent Pagneux, and Camille Perrot for their participation in my thesis defence. Having them with interesting questions and comments helped me to have a deeper vision and review of the research. I am also grateful to Nicholas Fang and Andrew Norris for having red the thesis and for their helpful comments and interests on the work. v vi Acknowledgements I have been truly fortunate to discuss through a meeting with Frédéric Hecht and Olivier Pironneau, who helped me solving problems related to a bug in FreeFem++ . The remarkable success of the nonlocal theory presented in chapter 3 by the results of Finite Element computations, could not be achieved without this meeting. I also thank Fabien Chevillotte, for providing me with the first examples of some programmes in FreeFem++, which helped me to learn efficiently this tool. The validation of the theory in chapter 3, could not be possible without collaboration of Aroune Duclos; specially concerning the precise results achieved by multiple scattering method. His presence either professionally and as a friend as well has been very productive for me.

I appreciated very much constructive discussions I had on my work during three years with the members of the group of Acoustics and Mechanics of Porous Materials in the LAUM. I would like to thank Bruno Brouard, Olivier Dazel, Jean-Philippe Groby, and Sohbi Sahraoui.

I also thank Bernard Castagnède, Claude Depollier, Michel Henry, Pascal Ruello, Christian Inguere, Mohamed Tabellout and Jean-Pierre Dalmont, for offering me interesting teaching activities, and their helps to have access easily to teaching materials.

The LAUM provided me with an excellent environment to pursue my interests in research and teaching. I thank all people, scientific and administrative staff, technicians and PhD students, who contributed to establish an exceptionally pleasant and warm atmosphere in the LAUM. In this respect, my special thanks go to the actual directors of the laboratory, Joël Gilbert and Laurent Simon, and also previous directors, Yves Aurégan and Rachid El Guerjouma. I thank also Anne-Marie Brulé, Valérie Hermann and Julie Béhu for their efficient administrative assistance.

During my stay at Ecole Polytechnique Fédérale de Lausanne (EPFL), in Switzerland, I greatly benefited from support and encouragements of Andrew Barry, Marc Parlenge and Farhad Rachidi-Haeri. Thanks must particularly be expressed to Jian Zhao, for guiding me in initial steps of research towards the field of mechanical wave propagation by proposing me to work on an interesting subject, which motivated me to do a PhD in this field subsequently. I enjoyed very much working with him and I am grateful to Jian for his admirable support.

Martine Ben Amar and Armand Ajdari have had an important presence in my educational and research life for some years; supported and motivated me regularly, and gave me valuable and important academic advices at challenging moments from the time I was performing my Master degrees in Paris.

Introduction

In this work, we present and validate a new macroscopic nonlocal theory of sound propagation in an unbounded rigid-framed porous media saturated with a viscothermal fluid. The porous medium is either isotropic or having a wave guide symmetry axis, and macroscopically homogeneous. Despite the success of the local theory to predict the macroscopic behaviour of long-wavelength sound propagation in such a medium with relatively simple microgeometries, it remains far to describe correctly the coarse-grained dynamics of the medium when the wave length become smaller than, or even comparable to a characteristic length of the medium, or when the microgeometry of the porous medium become more complex. These can be observed, practically, either in high frequency regime, or when the microgeometry is formed in a way to exhibit the resonance behaviour at not necessarily high frequencies. It will be shown throughout this work that the origin of these difficulties is that the classical theory, which is based on the asymptotic approach of the so-called two-scale homogenization method, is not capable to describe the macroscopic non local effects, contrary to the new proposed non local theory.

In the framework of the non local theory allowing for spatial dispersion, a new homogenization method is suggested by an analogy with Maxwell's electromagnetic theory and the establishment of a thermodynamic identification, in order to upscale the dynamics of viscothermal fluids from pore level to macroscopic level. The spatial dispersion is incorporated in the new theory by considering that in Fourier space, the the two acoustic susceptibilities -effective density and bulk modulus-, depend not only on the frequency, as in local theory, but also on the wavenumber. We provide an upscaling procedure to obtain these two effective properties of the medium in terms of the spatial-averaged values of the microscopic fields which can be computed through solving two independent action-response problems at the pore level. The theory is validated in several cases with different microgeometries.

We present the non local theory in chapter 1. It is shown that the equations governing the longitudinal wave motions in a viscothermal fluid can be put in a Maxwellian nonlocal form. The procedure to derive the corresponding nonlocal density and compressibility which is introduced there, then serves as a guide, when the solid structure is present, to obtain the correct nonlocal theory upscaling procedure.

In chapter 2, a successful test of this theory is presented in the simple case of cylindrical circular tubes filled with a viscothermal fluid . It is found that the wavenumbers and impedances predicted by nonlocal theory match with those of the long-known Kirchhoff's exact solution. On the contrary, the results by local theory (Zwikker and Kosten's), yield only the wavenumber of the least attenuated mode, with a small discrepancy compared to Kirchhoff's. Zwikker and Kosten's local theory is derived in Appendix A, using the language of temporal and spatial dispersion. To evaluate the Zwikker and Kosten quantities in chapter 2, the formulae reported in Appendix A are used.

In chapter 3, the nonlocal theory is verified in the case where the microgeometry of the porous medium is nontrivial, in the form of an unbounded two-dimensional square lattice of rigid cylinders permeated by a viscothermal fluid. On the one hand, we will compare the complex frequency-dependent phase velocity associated with the least attenuated plane wave, predicted by the new theory using Finite Element Method simulations, with that of the corresponding least attenuated Bloch mode, obtained by the quasi-exact multiple scattering method, and show that the two are in remarkable agreement. The main microscopic equations to be solved numerically, in order to compute macroscopic properties of the porous medium according to local and nonlocal theory, are briefly reviewed.

The essential elements of calculation to obtain Bloch wavenumbers through the multiple scattering method are also presented in this chapter.

In chapter 4, the theory will be validated when the the geometry of the porous medium is in the form of a daisy chained Helmholtz resonators. Firstly, an analytical plane wave modelling is employed to obtain the the least attenuated Block mode of the medium. Secondly, using the upscaling procedures in nonlocal theory and using the plane wave modelling lead to two effective density and bulk modulus functions in Fourier space. Once we have these functions, we are able to find, in particular, the least attenuated mode propagating in the medium, through a dispersion equation coming from the macroscopic equations in nonlocal theory.

The chapters 1 and 2, and 3, can be considered as relatively independent texts and have been written in a article style. Chapter 1 is the text of an article which has been submitted recently to the journal Wave Motion. Therefore, throughout

Chapter 1

Nonlocal theory of sound propagation in homogeneous rigid-framed porous media

1
Elaborating a Maxwellian representation of longitudinal wave propagation in a viscothermal fluid, a new general nonlocal macroscopic theory of sound propagation in homogeneous porous media saturated with viscothermal fluid is proposed here in the case where the solid frame is rigid and isotropic. Allowing in the most general manner for temporal and spatial dispersion this theory, contrary to the conventional two-scale homogenization theory, is suitable to homogenize the so-called metamaterials. Moreover, for propagation in periodic media along a symmetry axis, the complete Bloch mode spectrum is expected to be achieved without any frequency constraint.

Introduction

What are the equations governing small-amplitude sound propagation in a fluidsaturated porous material, at the macroscopic level? Here, we propose an answer to this question in the form of a new nonlocal theory of sound propagation, applicable in the case where the macroscopically homogeneous and isotropic material is rigid-framed and permeated by a viscothermal fluid. It is clear that a macroscopically homogeneous medium is also necessarily unbounded. Macroscopic homogeneity and isotropy are assumed here in Russakoff's sense of volume averaging [1].

1 Nonlocal theory of sound propagation in porous media

In the framework of this theory, we are able to upscale successfully the dissipative fluid dynamics in a medium with periodic microstructure, as well. In this case, the coarse-grained propagation is considered along a single symmetry axis.

The solid phase need not be necessarily connex. Thus the upscaling procedure is also valid for an homogeneous, isotropic or periodic distribution of motionless rigid scatterers of arbitrary shape, embedded in the viscothermal fluid.

The main existing homogenization method that may be used to predict the macroscopic properties of sound propagation in the above media is the socalled two-scale asymptotic homogenization method [2,3,4,5,6] which assumes wide scale separation between characteristic macroscopic wavelengths and typical correlation lengths of the solid and fluid phases. Applied to most sound absorbers used in practice, it consistently leads to local models taking into account temporal dispersion. The temporal dispersion effects are manifested in frequencydependent density and bulk modulus, describing the effective properties of the medium, which can be well-approximated by simple formulae involving a few geometrical parameters of the pore space [7,8,9].

This homogenization method, however, does not appear to give the general solution of the problem. For instance, in geometries with structures in the form of Helmholtz's resonators, it does not predict the characteristic metamaterial behavior related to the Helmholtz's resonance [10], which is expressed in a resonant bulk modulus. Moreover, in ordinary periodic geometries, it does not predict, beyond the Rayleigh scattering regime, the high-frequency behavior of the Bloch mode spectrum and the presence of the band gaps. Recently, powerful high-frequencies extensions of the conventional homogenization method have been introduced [11,12,13]; high-frequency meaning here the absence of two-scale separation in the usual sense of homogenization. These extensions may be suitable to homogenize periodic media in the vicinity of 'cell resonances', however, a new sort of two-scale separation is required once again.

The present general physical solution to the aforementioned wave propagation problems do not require any explicit scale separation. It takes advantage of an analogy with electromagnetics and a thermodynamic identification. In three forthcoming papers, it will be shown to successfully predict: the Kirchhoff's radial mode wavenumbers and impedances in a circular tube filled with a viscothermal fluid [14], the Bloch's wavenumbers and impedances in periodic square arrays of rigid cylinders permeated by a viscothermal fluid [15], and the metamaterial behavior of a line of daisy-chained Helmholtz's resonators [10].

The new feature in the theory which explains its absence of limitation for propagation in macroscopically homogeneous materials, is that it accounts both for temporal dispersion and spatial dispersion, without making any perturbative simplification or introducing any explicit scale-separation condition. The only scale separation that subsists is an unavoidable one. Given the finite spatial width of Russakoff's test function employed to smooth-out the material irregularities, the material properties are necessarily subjected to remaining small point-topoint fluctuations. The description of the propagation of macroscopic waves with wavelengths comparable to the correlation length of these fluctuations will no longer be possible in terms of an equivalent homogeneous frequencyand-wavenumber-dependent medium. This problem does not exists for periodic geometries, since these fluctuations entirely disappear by averaging over periods.

In this case, waves of typical wavelengths smaller than the irreducible cell dimension can be treated by the present nonlocal macroscopic theory.

The paper is organised as follows. Sound propagation in a viscothermal fluid is first revisited in section 1.2, where it is shown that the equations governing the longitudinal motions are susceptible to be put in a Maxwellian nonlocal form.

After making the appropriate thermodynamic identification which implies that the acoustic counterpart of the electromagnetic H field is the thermodynamic pressure field, we conclude that the longitudinal wave motions derive from an equivalent nonlocal density and an equivalent nonlocal compressibility. The nonlocal density which plays the role of electric permittivity depends only on inertial and viscous effects, and the nonlocal bulk modulus which plays the role of magnetic permittivity depends only on elastic and thermal effects.

It is remarkable that, once this thermodynamic identification of the acoustic H field is made, the corresponding density and bulk modulus operators are directly related to the solutions of two independent action-response problems. On the one hand, the density operator reflects the nonlocal response of the fluid subjected to an external force. On the other hand, the bulk modulus operator reflects the nonlocal response of the fluid subjected to an external rate of heat supply.

These fundamental observations are next directly used in section 1.3, when the viscothermal fluid is permeating a macroscopically homogeneous porous structure. Generalizing the electromagnetic recasting of the equations, and the aforementioned thermodynamic identification as well, we conjecture that there is, now, a macroscopic density operator which reflects the nonlocal response of the permeating fluid subjected to an external force, and a macroscopic bulk modulus operator which reflects the nonlocal response of the permeating fluid subjected to an external rate of heat supply. This leads to stating definite recipes for determining the operators from the microgeometry.

Electromagnetic recasting of the acoustic equations in a viscothermal fluid

Within the approximations used in ordinary near-equilibrium fluid-mechanics, small amplitude wave motions in viscothermal fluids can be analysed in terms of two disconnected types of motion: shear motions with transverse velocity variations and no condensation, pressure, and temperature variations, and longitudinal motions with longitudinal velocity variations and nonvanishing condensation, pressure, and temperature variations [16].

In section 1.3 we consider macroscopic sound propagation in a homogeneous, thus unbounded, fluid-saturated porous medium. In practice the material is finite and the macroscopic perturbation comes from a source placed in the external free fluid. The source may generate shear waves and longitudinal waves. However, since we intend to describe the propagation of sound waves in the material, we are not concerned by the shear waves which, by nature, involve no pressure variations.2

Thus in this section, we are interested to investigate only the longitudinal motions. Within the ordinary Navier-Stokes-Fourier model of a viscothermal fluid, the linearised equations of longitudinal motions are written as [19,20,21]

ρ 0 ∂v ∂t = -∇p + ( 4η 3 + ζ)∇(∇ • v) (1.1a) ∂b ∂t + ∇ • v = 0 (1.1b) γχ 0 p = b + β 0 τ (1.1c) ρ 0 c p ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ (1.1d)
where the wave variables v, b, p, τ , are the fluid velocity, condensation, thermodynamic excess pressure, and excess temperature respectively, and the fluid constants ρ 0 , η, ζ, γ, χ 0 , β 0 , c p , T 0 , κ, represent the ambient density, first viscosity, second viscosity, ratio of heat coefficients c p /c v , adiabatic compressibility, thermal expansion coefficient, specific heat coefficient at constant pressure, ambient temperature, and thermal conduction coefficient, respectively.

Regarding these equations we observe that, if we make abstraction of our knowledge that the condensation b represents the quantity ρ ′ /ρ 0 , where ρ ′ is the excess density, we may view the equation (1.1b) as a sort of definition of a field b which is interestingly similar, in its philosophy, to the macroscopic electromagnetic equation-definition valid in any material medium

∂B ∂t + ∇ × E = 0 (1.2)
To express it in more detail, we note that an enlightening way of looking at the above equation (1.2) is to observe that, to any macroscopic electromagnetic field present in a material medium -in a given rest state of thermodynamic equilibrium -may be associated a 3-vector potential macroscopic field A. In fact, the medium defines its own privileged rest frame in which we work, therefore the time component of the electromagnetic potential can always be set to zero by the jauge invariance. What we call, in this rest frame, the electric and magnetic fields E and B, are nothing but the quantities which are related to this 3-potential by the equation-definitions E = -∂A/∂t and B = ∇ × A.

Eq.(1.2) then is a direct consequence of the definitions of E and B in terms of the electromagnetic potential. The physics is expressed in the additional electromagnetic field equation

∂D ∂t = ∇ × H (1.3)
and the constitutive relations

D = ǫ E (1.4a) H = μ-1 B (1.4b)
where ǫ and μ are two constitutive permittivity operators in which the physical properties of the medium are encoded. For sufficiently small amplitudes and in an homogeneous medium, these are linear difference-kernel operators [22].

Here, the particle displacement field a which is a 3-vector as well, can be viewed as the acoustic counterpart of the 3-potential A from which, may be derived the fields v and b, by the equation-definitions v = ∂a/∂t and b = -∇•a, respectively. Eq.(1.1b) then is a direct consequence of these definitions.

This suggests that we can complete the Eq.(1.1b) by expressing the acoustic equations analogous to (1.3) and (1.4)

∂d ∂t = -∇h (1.5)
1 Nonlocal theory of sound propagation in porous media and

d = ρv (1.6a) h = χ-1 b (1.6b)
where ρ and χ are two constitutive linear difference-kernel operators in which the fluid physical properties are encoded.

The operator χ-1 or the field h can always be arbitrarily chosen. The physics then is wholly expressed by the operator ρ. Indeed, this is the customary choice utilized in electromagnetic theory in presence of spatial dispersion [START_REF] Landau | Electrodynamics of continuous media[END_REF]22,23], which consists in setting by definition μ-1 = µ -1 0 . Here we could set by definition χ-1 = χ -1 0 .

However, in electromagnetism a different choice exists which seems to be more natural, and which is based on the fact that the Poynting vector S = E × H is interpreted as the electromagnetic part of energy current density [START_REF] Fournet | Electromagnétisme à partir des équations locales[END_REF]. Although, this point of view is usually considered in the case where there is no spatial dispersion, we believe that one should go beyond and generalize its validity in presence of spatial dispersion. Thus, a general nontrivial H field exists which is defined through S.

Here, we borrow this concept from electromagnetism and bring it to acoustics. As such, the acoustic Poynting vector s

s = hv (1.7)
is interpreted as the acoustic part of energy current density. Taking for this acoustic energy current density the expression-definition suggested by Schoch [START_REF] Schoch | Remarks on the concept of acoustic energy[END_REF] s ≡ pv, the above condition requires that the field h identically matches the thermodynamic excess pressure p

h ≡ p (1.8)
Accounting for the fact that the operators ρ and χ1 must be expressed via difference-kernels functions, allowing to respect time and space homogeneity, the nonlocal Maxwellian acoustic equations then take the form

∂b ∂t + ∇ • v = 0 (1.9) ∂d ∂t = -∇p (1.10)
1.2 Electromagnetic recasting of the acoustic equations with

d(t, x) = t -∞ dt ′ dx ′ ρ(t -t ′ , x -x ′ )v(t ′ , x ′ ) (1.11a) p(t, x) = t -∞ dt ′ dx ′ χ -1 (t -t ′ , x -x ′ )b(t ′ , x ′ ) (1.11b)
We now intend to determine these kernel functions. It turns out that the direct determination of the functions ρ(t, x) and χ -1 (t, x) is not possible due to some mathematical singularities. Their physical origin is related to the missing terms in the Newton-type and Fourier constitutive laws introduced in the fundamental momentum and energy balance laws leading to (1.1a) and (1.1d). Nevertheless, working in Fourier space it is easy to find the Fourier kernels ρ(ω, k) = ρ(ω, k) and χ -1 (ω, k) = χ -1 (ω, k) such that the wave physics described by Eqs.(1.1a-1.1d), namely

iωρ 0 v = ikp + ( 4η 3 + ζ)k 2 v (1.12a) -ωb + k • v = 0 (1.12b) γχ 0 p = b + β 0 τ (1.12c) iωρ 0 c p τ = iωβ 0 T 0 p + κk 2 τ (1.12d)
be exactly the same as that described by Eqs.(1.9-1.11), namely

k • v = ωb (1.13a) kp = ωd (1.13b) with d = ρ(ω, k)v (1.14a) p = χ -1 (ω, k)b (1.14b)
Comparing the two sets of equations (1.12) and (1.13-1.14), we obtain the following expressions for the Fourier density and bulk-modulus kernels

ρ(ω, k) = ρ(ω, k) = ρ 0 1 + 4η 3 + ζ ρ 0 k 2 -iω (1.15) χ -1 (ω, k) = χ -1 (ω, k) = χ -1 0 1 - γ -1 γ 1 - iω -iω + κ ρ 0 cv k 2 (1.16)
where the following general thermodynamic identity [START_REF] Epstein | Textbook of thermodynamics[END_REF] has been used

γ -1 = T 0 β 2 0 ρ 0 χ 0 c p (1.17)
Although, this will not be our concern in this paper, we note that the above Maxwellian recasting cannot express totally its power within the too-simplified Newton and Fourier constitutive laws utilized here. Its degenerate character may be seen through the above expressions in different related manners.

First, the above expressions depend on k only via the modulus k. With more complete constitutive laws not leading to degeneracies, no complete decoupling would exist between shear and longitudinal motions, and ρ(ω, k) would be a tensor ρ ij (ω, k) having the general form [START_REF] Landau | Electrodynamics of continuous media[END_REF]

ρ ij (ω, k) = ρ t (ω, k)(δ ij -k i k j /k 2 ) + ρ l (ω, k)k i k j /k 2 (1.18)
where ρ t and ρ l are transversal and longitudinal kernels which depend only on the magnitude of the wave vector (and on ω). Here, the transverse part disappears and the longitudinal part reduces to ρ l (ω, k)δ ij , because they are contracted just by longitudinal velocities, which finally leads to a scalar ρ verifying ρ(ω, k) = ρ(ω, k).

Another way to see the existing degeneracies is to note that, for the reasons relating to the causality, the above expression for ρ(ω, k) should satisfy Kramers-Kronig dispersion relations [START_REF] Landau | Electrodynamics of continuous media[END_REF][START_REF] Landau | Statistical Physics[END_REF]]

ℜ[ρ(ω, k)] -ρ 0 = 1 π P.V. ∞ -∞ ℑ[ρ(ξ, k)] ξ -ω dξ (1.19) ℑ[ρ(ω, k)] = - 1 π P.V. ∞ -∞ ℜ[ρ(ξ, k)] ξ -ω dξ + σ(k) ω (1.20)
and similarly for χ -1 (ω, k), but without the pole. For the density ρ, these relations are verified in a trivial way. The only nonzero term in the right-hand side of (1.19) and (1.20) is the pole relating to σ

(k) = (4η/3 + ζ) k 2 .
For the bulk modulus χ -1 , the degeneracy is not directly apparent in the frequency and wavenumber dependencies (1.16). It is expressed by the very fact that we are dealing with a scalar rather than a tensor: the vectors H and B are replaced by scalars h (or p) and b. In reality, H should not be viewed as a vector but an antisymmetric tensor of rank two with contravariant indexes and weight W = -1 [START_REF] Fournet | Electromagnétisme à partir des équations locales[END_REF]. The field h should not be viewed as a scalar but a symmetric tensor of rank two with contravariant indexes and weight W = -1, akin to a contravariant stress field. That is precisely the tensorial character of the quantity h, which would allow ρ to be a tensor. In the same manner, B should not be regarded as a vector but an antisymmetric tensor of rank two with covariant indexes and weight W = 0. The field b should not be regarded as a scalar but a symmetric tensor of rank two with covariant indexes and weight W = 0, representing a covariant strain field. The nondegenerate tensor operator χ-1 that could be written within a sufficiently general thermodynamic framework for the fluid mechanics equations, would therefore become a tensor operator ( χ-1 ) ijkl . As such, the shear and compressional motions would be merged and treated simultaneously, instead of being artificially disconnected.

To get a glimpse of some of the deep consequences of the Maxwellian recasting sketched above, we notice that in electromagnetics, when the fields vary sufficiently rapidly, the dissipation processes have no time to occur and the general operator relations D i = ǫij E j and H ij = (μ -1 ) ijkl B kl reduce to the ones implying that D and H are determined by E and B in the same space-time position. With material homogeneity, isotropy, and center symmetry taken into account, the only available tensor is δ ij (in Cartesian coordinates) [START_REF] Landau | Electrodynamics of continuous media[END_REF] and the only relations compatible with the characteristic antisymmetry in H and B are,

D i = ǫ 0 δ ij E j and H ij = µ -1 0 1
2 (δ ik δ jl -δ il δ jk )B kl , i.e., they are proportionality relations D = ǫ 0 E, and H = µ -1 0 B, with ǫ 0 and µ 0 , the physical constants named electric and magnetic permittivity, respectively.

Here in acoustics, similarly, we suppose that when the fields vary sufficiently rapidly, the dissipation processes have no time to take place, and the general operator relations d i = ρij v j and h ij = ( χ-1 ) ijkl b kl reduce to the ones implying that d and h are determined by v and b in the same space-time position. As before, the only available tensor is δ ij , and the only relations compatible with the characteristic symmetry in h and b are

d i = ρ 0 δ ij v j and h ij = χ -1 0 δ ij δ kl b kl + µ 0 (δ ik δ jl +δ il δ jk -2 3 δ ij δ kl )b kl , with ρ 0 , χ -1
0 , and µ 0 , physical constants interpreted as density, compressibility modulus, and rigidity modulus, respectively [START_REF] Landau | Theory of Elasticity[END_REF]. It may be checked that the latter identically vanishes within the degenerate recasting we discuss here, whereas the former two yield the ambient density and adiabatic bulk modulus.

Regarding the above considerations, it seems that Frenkel's idea in phonon theory of fluid thermodynamics [START_REF] Frenkel | Kinetic Theory of Liquids[END_REF], reported long ago, stating that fluids should behave like elastic solids at very short times, is automatically recovered. In this connection we note that, recently, a phonon theory of thermodynamics of liquids has been developed using Frenkel's theoretical framework of the phonon states in 1 Nonlocal theory of sound propagation in porous media a liquid, which has shown good agreement for the calculations of heat capacity coefficients of different liquids in a wide range of temperature and pressure [START_REF] Bolmatov | The phonon theory of liquid thermodynamics[END_REF].

Although, the present Maxwellian description shows to have degeneracies and, apparently, could be well-defined only in a larger fluid-mechanics thermodynamic framework, we shall not care about this issue, as our purpose is to use this theory to have a generalized description of sound propagation in fluid-saturated porous materials. In this case, describing the fluid behavior at the level of the ordinary Navier-Stokes-Fourier equations, will be entirely sufficient for the noise control applications of the theory.

The fact that the correct longitudinal Navier-Stokes-Fourier wave physics [16] is encoded in the degenerate expressions (1.15-1.16), can be easily checked by writing the dispersion equation corresponding to the system (1.13-1.14)

ρ(ω, k)χ(ω, k)ω 2 = k 2 (1.21)
Substituting in the above equation, the expressions (1.15-1.16), we get the known Kirchhoff-Langevin's dispersion equation of longitudinal waves [20,19] -

ω 2 + c 2 a -iω κ ρ 0 C V + 4η 3 + ζ ρ 0 k 2 - κ ρ 0 C V iω c 2 i -iω 4η 3 + ζ ρ 0 k 4 = 0 (1.22)
with c 2 a the adiabatic sound speed squared defined by c 2 a ≡ 1/ρ 0 χ 0 , and c 2 i the isothermal sound speed squared defined by c 2 i ≡ c 2 a /γ. Therefore, it is verified that the correct sound normal modes, and heat conduction normal mode [16], are encoded in the given expressions.

We arrive now at the fundamental two important elements which will be used in the next section. Let e be the direction along which we study the sound propagation, and x = x•e be the coordinate along this direction. It will be shown next that the two functions ρ(ω, k) and χ -1 (ω, k) are related to the solutions of two independent action-response problems obtained by putting, respectively, a fictitious harmonic pressure term P(t, x) = P 0 e -iωt+ike•x in the Navier-Stokes Eq.(1.1a), or the Fourier Eq.(1.1d).

Firstly, adding the potential bulk force f = -∇P = -ikPe to the right-hand side of Eq.(1.1a) belonging to the equation system (1.1a-1.1d), and writing the 1.2 Electromagnetic recasting of the acoustic equations fields as

v(t, x) = ve -iωt+ike•x (1.23a) b(t, x) = be -iωt+ike•x (1.23b) p(t, x) = pe -iωt+ike•x (1.23c) τ (t, x) = τ e -iωt+ike•x (1.23d)
we can easily get the response amplitudes v, b, p, τ . We then observe that the same expression as in (1.15) for ρ(ω, k) is obtained through the equation

-ρ(ω, k)iωv = -ike(p + P 0 ) (1.24)
In this problem, the response pressure p is added to the fictitious deriving pressure amplitude P 0 to represent a sort of total effective pressure field h.

This establishes a direct relation between the Fourier coefficient ρ(ω, k) (1.15) of the operator density, and the response of the fluid subjected to an external harmonic bulk potential force.

Secondly, putting the bulk rate of heat supply Q = β 0 T 0 ∂P/∂t = -iωβ 0 T 0 P in the right-hand side of Eq.(1.1d) belonging to the equation system (1.1a-1.1d), and writing the fields as before, we get the response amplitudes, v, b ′ , p, τwith a prime on the condensation for later convenience. We then observe that the same expression as in (1.16) for χ -1 (ω, k) is obtained through the equation

p + P 0 = χ -1 (ω, k)(b ′ + γχ 0 P 0 ) (1.25)
In this problem, the response pressure p is again added to the fictitious deriving pressure amplitude P 0 to represent a sort of total effective pressure field h, and the term γχ 0 P 0 is added to the response condensation b ′ , in order to represent a sort of total field b, related to h by the constitutive relation (1.6b). The interpretation of the corresponding decomposition of b in two terms is that b ′ is a nonisothermal part of the response, while γχ 0 P 0 is an isothermal response part; γχ 0 being the isothermal compressibility.

This establishes in turn a direct relation between the Fourier coefficient χ -1 (ω, k) (1.16) of the operator bulk modulus, and the response of the fluid subjected to an external harmonic bulk rate of heat supply.

These relations are employed directly in the next section to arrive at the wanted new macroscopic theory of sound propagation.

Generalization to fluid-saturated rigid-framed porous media

We consider, now, that the viscothermal fluid which is pervading the connected network of pores of a rigid-framed porous material, is set to vibrate following an incident small-amplitude sound wave. Assume that a characteristic length L exists, allowing to smooth out the irregularities of the porous structure. The aim is to describe the sound propagation in this complex medium.

As stated in general terms, this problem directly reminds the linear theory of electromagnetic wave propagation in matter, investigated by H.A. Lorentz in his theory of electrons [17]. In the same way as Lorentz could not hope to follow in its course each electron, here we do not intend to analyse the detail of the wave propagation at the pore level. Following Lorentz we remark that it is not the microlevel wavefield that can make itself felt in the experiments, that are carried out at the macroscopic level, but only the resultant effect produced by some macroscopic averaging. A macroscopic description of the sound propagation in the medium will be possible if we fix from the outset our attention not on the pore level irregularities, but only on some mean values. We proceed now to clarify it through some definitions.

Basic definitions

The homogeneous porous medium occupies the whole space and is composed of two regions: the void (pore) region V f which is a connex region permeated by the fluid, and a solid-phase region V s . The pore-wall region or solid-fluid interface is denoted by ∂V. The characteristic function of the pore region is defined by

I(x) = 1, x ∈ V f 0, x ∈ V s (1.26)
Given a field a(t, x) in the fluid, such as the velocity field, a macroscopic mean value A = a may be defined by volume-averaging in a sphere of typical radius L/2, as was done by Lorentz [17]. Since the purpose of the averaging is to get rid of the pore-level irregularities, the length L is taken sufficiently large to include a representative volume of the material. This Lorentz's averaging was subsequently refined by Russakoff [1] who replaced the integration in a sphere by a convolution with a smooth test function

f L of characteristic width L A(t, x) = a (t, x) = dx ′ I(x ′ )a(t, x ′ )f L (x -x ′ ) (1.27)
which is a much better smoothing procedure from a signal analysis standpoint, especially on account of the fact that the material homogeneity is in general only approximately realized over the distance L.

The presence of characteristic function I in (1.27) ensures that the integration is taken in the only fluid region. The test function f L is normalized over the whole space

dxf L (x) = 1 (1.28)
The so-called spatial averaging theorem [18] is written as

∇a = ∇ a + ∂V dx ′ a(t, x ′ )n(x ′ )f L (x -x ′ ) (1.29)
relating the average of the gradient of a microscopic field a to the gradient of the averaged field. The macroscopic homogeneity implies that the quantity

φ = I (1.30)
is a constant independent of the position. It represents the fluid volume fraction or the porosity of the medium. Thus for a macroscopically homogeneous medium, the commutation relation (1.29) results in

∂V dx ′ n(x ′ )f L (x -x ′ ) = 0 (1.31)
The constancy of φ or the vanishing of the integral (1.31) cannot exactly be realized because of the irregularities of the medium; some inherent indeterminacy is generally associated with the notion of spatial averaging. In practice, however, the integral (1.31) often decreases sufficiently to be neglected for all practical purposes, as soon as the averaging length L/2 reaches the value of some typical coarse graining radius. This is the situation implicitly considered here.

Spatial periodicity is an important particular case where exact homogeneity (1.31) is gained, but at the cost of an indefiniteness in the choice of the averaging length L. Indeed, considering wave propagation along a single symmetry axis, and denoting by P the irreducible period peculiar to this axis, we may take for the smoothing test function f L , any of the following functions f n , n = 1, 2, ...

f n (x) = f n (x)g(y, z), with dydz g(y, z) = 1, dx f n (x) = 1
1 Nonlocal theory of sound propagation in porous media and

f n (x) = 1/nP, |x| < nP/2 0, |x| > nP/2
These functions in turn define a series of different macroscopic averaged fields a n performed with coarse-graning averaging length L = nP . Now that the volume-average operation has been defined including the above points concerning the periodic case, given a small amplitude wave perturbation to the ambient equilibrium in the fluid, it is time to establish a Lorentz macroscopic theory of wave propagation of averaged wavefield quantities.

Pore-level equations

At the pore-level the linearized fluid-mechanics equations describing the fluid motion are written as

ρ 0 ∂v ∂t = -∇p + ( 4η 3 + ζ)∇(∇ • v) -η∇ × (∇ × v) (1.32a) ∂b ∂t + ∇ • v = 0 (1.32b) γχ 0 p = b + β 0 τ (1.32c) ρ 0 C P ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ (1.32d) in V f , and v = 0 (1.33a) τ = 0 (1.33b) on ∂V.
Comparing to the preceding equations in the free fluid (1.1a-1.1d), a rotational viscous term has been added in (1.32a) to account for the shear motions generated at the pore walls by virtue of the no-slip condition (1.33a). The condition (1.33b) comes from the fact that the solid frame which is much more thermally-inert than the fluid, is assumed to remain at the ambient temperature.

Macroscopic equations

Given the pore-level equations, we seek the macroscopic equations governing wave propagation of the averaged quantities

V ≡ v , and B ≡ b (1.34)
Since the velocity vanishes at the pore walls, the following direct commutation relation always holds true

∇ • v = ∇ • v = ∇ • V (1.35)
Thus, the Eq.(1.32b) is immediately translated at the macroscopic level

∂B ∂t + ∇ • V = 0 (1.36)
The electromagnetic analogy then suggests that the system of macroscopic equations can be carried through by introducing new Maxwellian fields H and D, and also operators ρ and χ-1 , such that

∂D ∂t = -∇H (1.37) with D = ρV (1.38a) H = χ-1 B (1.38b)
As we have seen in the previous section, such form of equations, with the scalar H and scalar ρ, is suitable to treat nonlocal propagation of longitudinal waves in an isotropic medium. Assuming scalar H and scalar ρ, we disregard the propagation of macroscopic shear waves. Another case that is adapted to this form concerns the propagation of longitudinal waves along a symmetry axis x set in the direction e. In the former case, accounting for the time homogeneity and the material macroscopic homogeneity, the nonlocal relations (1.38) are written

D(t, x) = t -∞ dt ′ dx ′ ρ(t -t ′ , x -x ′ )V (t ′ , x ′ ) (1.39) H(t, x) = t -∞ dt ′ dx ′ χ -1 (t -t ′ , x -x ′ )B(t ′ , x ′ ) (1.40)
1 Nonlocal theory of sound propagation in porous media with scalar kernels ρ(t, x) and χ -1 (t, x) whose Fourier amplitudes verify ρ(ω, k) = ρ(ω, k), and χ -1 (ω, k) = χ -1 (ω, k). In the latter case with a symmetry axis along e, we have D = De and V = V e, then the constitutive relations become

D(t, x) = t -∞ dt ′ dx ′ ρ(t -t ′ , x -x ′ )V (t ′ , x ′ ) (1.41) H(t, x) = t -∞ dt ′ dx ′ χ -1 (t -t ′ , x -x ′ )B(t ′ , x ′ ) (1.42)
with scalar kernels ρ(t, x) and χ -1 (t, x) whose Fourier amplitudes are ρ(ω, k), and χ -1 (ω, k).

In the above, the integrations over time t ′ in one hand, and over space coordinate x ′ or x ′ in another hand, convey the image of the so-called temporal dispersion effects and spatial dispersion effects, respectively. They can be served here as definitions for these effects.

We need now to identify the macroscopic field H. provided that S = pv . Thus H is identified as the field satisfying the relationdefinition

H v ≡ pv (1.44)
This identification of an effective macroscopic pressure field H different from the usual fluid-volume-averaged mean pressure p = 1 φ p , appears very natural. Generalized to the case of a bounded open material, it yields a field H that is continuous at the macroscopic interface of the material. This follows from the continuity of the normal component of the field V which is required by mass flow conservation and the continuity of the normal component of the acoustic part of the energy current density that may be supposed to hold true, provided no resistive surface layer exists at the boundary of the material. Thus, in the present electromagnetic-acoustic analogy, it may be noted that the continuity of the normal component of velocity V replaces the continuity of the tangential components of the field E, and the continuity of the scalar H replaces the continuity of the tangential components of H [START_REF] Landau | Electrodynamics of continuous media[END_REF].

Identification of constitutive operators

Generalizing the relations established for the longitudinal motions in viscothermal fluids to the present case, leads us to suggest now that the above Fourier coefficients ρ(ω, k) and χ -1 (ω, k) are directly related to the macroscopic response of the permeating fluid subjected to a harmonic fictitious pressure term P(t, x) = P 0 e -iωt+ike•x added to the pressure, either in the Navier-Stokes Eq.(1.32a), or the Fourier Eq.(1.32d).

Thus to determine the kernel ρ(ω, k) we first consider solving the action-response problem

ρ 0 ∂v ∂t = -∇p + ( 4η 3 + ζ)∇(∇ • v) -η∇ × (∇ × v) + f (1.45a) ∂b ∂t + ∇ • v = 0 (1.45b) γχ 0 p = b + β 0 τ (1.45c) ρ 0 c P ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ (1.45d) in V f , and v = 0 (1.46a) τ = 0 (1.46b)
on ∂V. The external force appears as before in the form of

f = -∇P = -ikeP 0 e -iωt+ike•x (1.47)
The unique solutions to the above system (1.45a-1.47), for the fields v, b, p, τ , take the form

v(t, x) = v(ω, k, x)e -iωt+ike•x (1.48a) b(t, x) = b(ω, k, x)e -iωt+ike•x (1.48b) p(t, x) = p(ω, k, x)e -iωt+ike•x (1.48c) τ (t, x) = τ (ω, k, x)e -iωt+ike•x (1.48d)
The response amplitudes v(ω, k, x), b(ω, k, x), p(ω, k, x), and τ (ω, k, x) are bounded functions which are uniquely determined by the microgeometry.

1 Nonlocal theory of sound propagation in porous media

The above problem, once solved, we can extract of the response pressure p(t, x) = p(ω, k, x)e -iωt+ike•x , its macroscopic part, denoted by

P(t, x) = P(ω, k)e -iωt+ike•x (1.49)
whose amplitude P(ω, k) is determined through the equation

P(ω, k) = p(ω, k, x)v(ω, k, x) .e v(ω, k, x) .e (1.50)
This expression comes from the relation-definition

P v = pv (1.51)
which has been inspired by (1.44). Then using the following relation

-iωρ(ω, k) v = -
ρ(ω, k) = k(P(ω, k) + P 0 ) ω v(ω, k, x) .e (1.53)
As this point, we see that the fields p(ω, k, x) and v(ω, k, x) are needed to be known in order to determine from microgeometry the effective density of the fluid-saturated porous medium. Hence, instead of solving (1.45a-1.46b) we have just to solve the following system of equations to get the amplitudes of the fields in (1.48)

-iωρ 0 v = -(∇ + ike)p + ( 4η 3 + ζ)(∇ + ike) (∇ • v + ike • v) (1.54a) -η(∇ + ike) × (∇ + ike) × v -ikeP 0 -iωb + ∇ • v + ike • v = 0 (1.54b) γχ 0 p = b + β 0 τ (1.54c) -ρ 0 c P iωτ = -β 0 T 0 iωp + κ(∇ + ike) • (∇ + ike)τ (1.54d) in V f , and v = 0 (1.55a) τ = 0 (1.55b)
on ∂V.

The procedure to determine the kernel χ -1 (ω, k) is quite similar. We consider again, initially solving the action-response problem with an excitation appearing this time in the energy balance equation, such that

ρ 0 ∂v ∂t = -∇p + ( 4η 3 + ζ)∇(∇ • v) -η∇ × (∇ × v) (1.56a) ∂b ′ ∂t + ∇ • v = 0 (1.56b) γχ 0 p = b ′ + β 0 τ (1.56c) ρ 0 c P ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ + Q (1.56d) in V f , and v = 0 (1.57a) τ = 0 (1.57b)
on ∂V, with the external heating

Q = β 0 T 0 ∂P ∂t = -iωβ 0 T 0 P 0 e -iωt+ike•x (1.58)
The solutions to the above problem take the same form as specified before through Eqs. (1.48). As previously done, we extract of the response pressure p(t, x) = p(ω, k, x)e -iωt+ike•x , its macroscopic part P(t, x) = P(ω, k)e -iωt+ike•x whose amplitude P(ω, k) is determined by the equation

P(ω, k) = p(ω, k, x)v(ω, k, x) .e v(ω
P(ω, k) + P 0 = χ -1 (ω, k) b ′ (ω, k, x) + φγχ 0 P 0 (1.60)
where the factor of porosity is inserted in the last term to account for the fact that P 0 = φP 0 . That gives rise to nonlocal Equivalent-fluid bulk modulus χ -1 (ω, k)

χ -1 (ω, k) = P(ω, k) + P 0 b ′ (ω, k, x) + φγχ 0 P 0 (1.61)
Obviously, we need to know the fields p(ω, k, x) and v(ω, k, x), and b ′ (ω, k, x) in order to determine the effective bulk modulus of the fluid-saturated porous medium. Substituting (1.48) in (1.56a-1.57b), we have subsequently the following 1 Nonlocal theory of sound propagation in porous media system whose solutions are the field amplitudes we seek

-ρ 0 iωv = -(∇ + ike)p + ( 4η 3 + ζ)(∇ + ike) (∇ • v + ike • v) (1.62a) -η(∇ + ike) × (∇ + ike) × v -iωb ′ + ∇ • v + ike • v = 0 (1.62b) γχ 0 p = b ′ + β 0 τ (1.62c) -ρ 0 c P iωτ = -β 0 T 0 iωp + κ(∇ + ike) • (∇ + ike)τ -iωβ 0 T 0 P 0 (1.62d) in V f ,

Periodic media

In the special case of periodic media the response amplitudes (1.48), in both density and bulk modulus relating action-response problems, are now periodic functions of x = x • e, where e is the symmetry axis along which the propagation is considered. To fix the solution, it is necessary to precise its supercell irreducible period L, which can be any integral multiple L = nP of the geometric irreducible period P of the medium along direction e. Thus we now write, e.g., for the velocity solution field

v n (t, x) = v n (ω, k, x)e -iωt+ike•x (1.64)
with the periodicity condition

v n (ω, k, x + nP e) = v n (ω, k, x) (1.65)
stated on the boundary of the supercell, being understood that nP is an irreducible period of the considered amplitude, and so on for the other fields p n , b n , and τ n .

Moving on to the determination procedures of the operators, there follows that a discrete infinite set of branches of an Equivalent-fluid density ρ(ω, k, n) are produced, such that

P n (ω, k) = p n (ω, k, x)v n (ω, k, x) n .e v n (ω, k, x) n .e (1.66)
and

ρ n (ω, k) = -ik(P n (ω, k) + P 0 ) -iω v(ω, k, x) n .e (1.67)
and likewise, a discrete infinite set of branches of an Equivalent-fluid bulk modulus χ -1 (ω, k, n) are produced, such that

χ -1 n (ω, k) = P n (ω, k) + P 0 b ′ n (ω, k, x) n + φχ 0 γP 0 (1.68)
where the symbol n represents the averaging over the supercell constituted of n geometric periods P in the propagation direction.

Characteristic wavenumbers and impedances of the normal mode solutions

The characteristic feature of the present macroscopic theory is that, without any simplifications, it allows for both temporal and spatial dispersion. Within the classical local equivalent-fluid theory which only accounts for temporal dispersion, for a given frequency ω, there is only one single normal mode that can propagate in the given positive x direction. With this single mode is associated a wavenumber q(ω) verifying the relation q 2 = ρ(ω)χ(ω)ω 2 , ℑ(q) > 0, where ρ(ω) and χ(ω) are the local density and compressibility functions. Here, since we fully take into account the spatial dispersion, several normal mode solutions might exist, with fields varying as e -iωt+iqe•x . At this time, each of these solutions should satisfy the following dispersion equation

ρ(ω, q)χ(ω, q)ω 2 = q 2 (1.69)
If we label l = 1, 2, ..., the different solutions q l (ω), ℑ(q l ) > 0, to the Eq.(1.69), the corresponding wave impedances Z l = H/V for propagation in the direction +x will be Z l (ω) = ρ(ω, q l (ω))/(ω 2 χ(ω, q l (ω))).

When the geometry is periodic, the index n will have to be added in these formulae. In our forthcoming papers, we show on different examples, that the present theory predicts the adequate wavenumbers q l (ω) and impedances Z l (ω). We note that the nonlocal functions ρ(ω, k) and χ(ω, k) will not be found to 1 Nonlocal theory of sound propagation in porous media systematically reduce to the local ones at long wavelengths. The asymptotic relations of the type ρ(ω) = lim k→0 ρ(ω, k), χ(ω) = lim k→0 χ(ω, k), which are in line with the perturbative philosophy of conventional homogenization theory and are also often assumed without proof for the electric permittivity [22], will not hold true in the vicinity of Helmholtz's resonances.

Conclusion

Following a general line of reasoning inspired by the electromagnetic theory, we have proposed a definite new procedure to perform nonlocal homogenization for sound propagation in rigid-framed homogeneous unbounded porous media saturated with a viscothermal fluid, either for general isotropic materials, or for propagation along a symmetry axis in periodic material. This procedure has been clarified first in a viscothermal fluid, leading to the Kirchhoff-Langevin's dispersion equation, and has been next generalized when a rigid porous matrix is embedded in this fluid. Contrary to the usual two-scale homogenization approach [2], or its recent high-frequency extensions [11,12,13], the proposed theory makes no perturbative simplification and will be valid as long as the medium can be considered macroscopically homogeneous in Lorentz-Russakoff's volume-average sense.

In practice, the randomness of the medium, which cannot entirely be smoothedout by the averaging procedure, will prevent constructing meaningful macroscopic description of the propagation, at wavelengths sufficiently small compared to typical correlation lengths involved in the medium properties fluctuations. But for periodic media, exact smoothing of microscopic variations can be achieved, and the nonlocal macroscopic theory will not be limited. Spatial variations of averaged fields, much smaller than the period, will be described as well by the macroscopic theory. There, a longstanding misconception might be expressed, once again, in connection with the notion of Lorentz-Russakoff's averaging. Contrary to what Lorentz mentioned in his theory of electrons, it is not necessary that the radius of the averaging sphere be so small that the state of the body, so far as it is accessible to our means of observation, may be considered as uniform throughout the sphere. Therefore, one does not need to assume that macroscopic wavelengths must be significantly greater than the length L. The ideal periodic case shows that no scale separation is required to establish a macroscopic theory.

We believe that the present nonlocal acoustic theory also indicates some misconceptions in electromagnetics regarding the nature of the field H in presence of spatial dispersion. Contrary to what is often asserted, a meaningful H field certainly exists in presence of spatial dispersion, even if we lack the adequate knowledge of thermodynamics necessary to define it. This problem is not encountered in acoustics which is a theory supported by near-equilibrium thermodynamics.

It is also suggested that the acoustic-electromagnetic analogy that we have used here, is a degenerate version of a much deeper one. The deeper version might be totally consistent with Frenkel's long overlooked idea [START_REF] Frenkel | Kinetic Theory of Liquids[END_REF][START_REF] Bolmatov | The phonon theory of liquid thermodynamics[END_REF], stating that a fluid behaves like a solid at very short times. The fluid capability to support solid-like shear waves at very short times is automatically implied by the tensor symmetry of the equations and the electromagnetic analogy.

Introduction

In a recent paper [1], by using Kirchhoff-Langevin's description of compressional wave propagation in a fluid, and an electromagnetic analogy as a powerful heuristic guide, we introduced two new upscaling procedures allowing to compute two 'acoustic permittivities' from microstructure of porous media. These describe in a nonlocal 'Maxwellian' manner, the phenomenon of macroscopic sound propagation in rigid-framed homogeneous porous materials saturated with a viscothermal fluid. The first is a macroscopic effective density playing the role of macroscopic electric permittivity, the second a macroscopic effective compressibility playing the role of macroscopic magnetic permittivity. By 'macroscopic' we mean that the theory is not concerned with the values of the acoustic fields at every microscopic spatial position, but only, with their 'macroscopic' effective values, obtained by averaging in a way to be precised. In the proposed theory, an energetic 'Umov-Poynting' definition of the macroscopic 34 2 Nonlocal theory of sound propagation in a circular pore pressure is introduced, and the permittivities are fully nonlocal, i.e. they are frequency and also wavenumber dependent.

The physical motivation of the new nonlocal theory is the recognition that spatial dispersion effects are not well described by the existing theories. Indeed, at the zero'th order of the asymptotic two-scale homogenization theory [2] and in all existing macroscopic models such as [3] and [4,5], spatial dispersion effects are entirely absent. In the full asymptotic two-scale homogenization theory [2,6], or its recent variants [7,8,9], some spatial dispersion effects are present, but in a limited manner.

The limited possible uses of the latter and other asymptotic methods, were in recent years highlighted by their inability to cope simultaneously with all geometries and frequencies, and in particular, to describe the whole dynamics of metamaterial structures with Helmholtz resonators [10] on one hand, and the complete Bloch mode spectrum in periodic structures on the other hand. With the proposed new 'Maxwellian' approach, these limitations disappear and no restrictions on periodic geometries or frequencies subsist, even if we have so far formulated the solution only for the case of propagation along a symmetry axis of an unbounded material. Therefore, it is important to note that the new theory is suitable to predict the exact properties of the so-called metamaterials [10], from microstructure.

As the new nonlocal theory is intended to provide the true physico-mathematical solution of the macroscopic wave propagation problem, it needs to be mathematically checked in unequivocal precise manner. A proper checking is especially desirable, also, because a detailed verification of the general ideas of the theory cannot be performed in electromagnetics, where the counterpart of the upscaling procedures proposed cannot be formulated yet.

In electromagnetics, the nonlocal Maxwell macroscopic theory comparable to the present acoustic one, is elusive. We believe that its very completion is not possible yet, because the necessary thermodynamics of electromagnetic fields in matter is missing.

In this paper we concentrate on a simple example for which an exact mathematical verification of the general physical considerations used to construct the nonlocal theory is possible. This example concerns sound propagation in straight ducts.

For ducts of circular cross-section the exact solution, accounting in the framework of near-equilibrium ordinary fluid mechanics for the effects of viscous losses and thermal conduction in the fluid, is known since G. Kirchhoff [11]. Kirchhoff's investigation had been in the thermodynamic framework of the ideal gas theory. Langevin [12,13] later showed that Kirchhoff's solution applied more generally, to a viscothermal fluid obeying an arbitrary equation of state. This available long known solution of the viscothermal wave propagation problem, offers the possibility to directly check in simple but nontrivial manner the nonlocal theory's physical considerations.

This paper is organized as follows. In section 2.2 we present Kirchhoff-Langevin's solution of the problem of small-amplitude sound propagation in a tube of circular cross-section filled with a viscothermal fluid. We recall how this solution allows to compute at given real angular frequency ω, the complex wavenumbers k (m,n) (ω) of the axisymmetric normal modes m = 0, n = 0, 1, 2, ..., where m and n are azimuthal and radial mode indexes. The corresponding wavenumbers k l (ω) = k (0,l-1) (ω), l = 1, 2, ..., are obtained as the complex roots of the transcendent Kirchhoff-Langevin's dispersion Eq.(2.47).

Anticipating on the Maxwellian theory's definition of a macroscopic pressure field H by means of the fundamental thermodynamic equation-definition

pv = H v (2.1)
-which we call the 'Umov-Poynting' definition since pv is interpreted as the acoustic part of macroscopic energy current density -where p is excess thermodynamic pressure, v velocity, and is the averaging operation over a cross-section, we then introduce at a given real angular frequency ω, the following complex impedance factors Z l (ω), l = 1, 2, ...

H = Z l v • e x (2.2)
where e x represents the unit vector along the x-axis.

These frequency-dependent Kirchhoff-Langevin's wavenumbers k l (ω) and impedances Z l (ω) enable us to evaluate two frequency-dependent Kirchhoff-Langevin's permittivities, namely the densities ρ l (ω) and bulk modulii χ -1 l (ω) associated with the different radial modes

n = l -1, l = 1, 2, ... ρ l (ω) = k l Z l /ω, χ -1 l (ω) = ωZ l /k l (2.3)
In section 2.3 we recall the principles of the proposed macroscopic Maxwellian nonlocal theory. There, two permittivities are introduced which are fundamentally nonlocal. They are expressed through two nonlocal density and bulk modulus operators, also described in Fourier space by frequency-dependent and wavenumber-dependent complex amplitudes ρ(ω, k) and χ -1 (ω, k). These func-tions ρ(ω, k) and χ -1 (ω, k) are independently computable thanks to the two conjectured upscaling procedures.

It will be shown that the Kirchhoff-Langevin's complex wavenumbers and impedances k l and Z l still make sense in the framework of the macroscopic theory.

The wavenumbers k l are the solutions of the 'Maxwell-Kirchhoff' dispersion equation

ρ(ω, k)χ(ω, k)ω 2 = k 2 (2.4)
and the impedances Z l then can be computed by

Z l = ρ(ω, k l )χ -1 (ω, k l ) (2.5)
They may be referred as Maxwell-Kirchhoff's wavenumbers and impedances.

The coincidence between Kirchhoff-Langevin's and Maxwell-Kirchhoff's complex wavenumbers and impedances serves as a test of the exactness of the two upscaling procedures of the new theory. Equivalently this can be expressed through the coincidence of Kirchhoff-Langevin's and Maxwell-Kirchhoff's densities and bulk modulii

ρ l (ω) = ρ(ω, k l ) χ -1 l (ω) = χ -1 (ω, k l ) (2.6)
A successful numerical checking of the above performed in section 2.4 will clearly indicate that the two nonlocal upscaling procedures described in [1] are exact, irrespective of the microgeometry. In forthcoming papers, similar successful numerical verifications will be made when the geometry is nontrivial but sufficiently simple to allow again for relatively easy detailed solutions. This will complete the numerical demonstration of the exactness and generality of the nonlocal upscaling procedures conjectured in [1].

We notice that Kirchhoff-Langevin's theory has been scarcely used in practice.

For the least attenuated plane mode which is most often the only one of importance, it gives the results indistinguishable from those of the simple approximate theory developed much later by Zwikker and Kosten [14]. Extended to arbitrary geometries this simpler treatment is nothing but that of the zero'th order homogenization theory. The propagation models ordinarily used in acoustic studies of porous media [4,5,15], are all developed in this local-theory framework.

In electromagnetism, it would correspond to the widely used simplification which consists in assuming that the permittivities have no dependencies on the wavenumber. This simplification, however, is clashing with the complete wave nature of the problem which implies to some extent spatial dispersion as well as temporal dispersion. Even if Zwikker and Kosten's theory works very well for the least-attenuated plane wave mode, it never provides the complete physical solution of the macroscopic (i.e. cross-section averaged) wave propagation problem. In particular, it lacks predicting the existence of the higher order modes, and also, the change of nature of the propagation in very wide tubes, illustrated by the change of nature of the fundamental mode, which tends to become surface-wave like.

Zwikker and Kosten's local theory is derived in Appendix A, using the language of temporal and spatial dispersion. It will be clear that the present complete 'Maxwell-Kirchhoff' description of sound propagation in a straight circular duct, is nothing but a Zwikker and Kosten's treatment generalized to include spatial dispersion.

Kirchhoff-Langevin's theory of sound propagation in a tube of circular cross-section

Kirchhoff's original investigations on the effects of viscosity and heat conduction on sound propagation in free air, and also, air inside a hollow solid tube, were conducted by treating the air as an ideal gas [11]. Langevin completed much later Kirchhoff's theory by considering air having the second viscosity and a general equation of state [12]. We commence by recalling this complete Kirchhoff-Langevin's theory, which is usually not presented without simplifications in the acoustic literature.

Linearized equations, in the Navier-Stokes-Fourier model

We are given a homogeneous viscothermal fluid which obeys an arbitrary caloric equation of state

ε = ε(s, υ) (2.7)
with ε the specific internal energy per unit mass, s the specific entropy, and υ = 1/ρ the specific volume. The thermodynamic pressure p, absolute temperature T , and specific heats at constant pressure and constant volume c p and c v are defined by

p ≡ - ∂ε ∂υ s , T ≡ ∂ε ∂s υ , c p ≡ T ∂s ∂T p , c v ≡ T ∂s ∂T υ (2.8)
2 Nonlocal theory of sound propagation in a circular pore Eqs.(2.7) and (2.8) and elimination of s results in a thermal equation of state p = p(T, υ). The fluid thermal expansion coefficient is defined by

β ≡ 1 υ ∂υ ∂T p (2.9)
Let us define two fixed reference sound speed values c a and c i , resp. adiabatic and isothermal, by

c 2 a ≡ ∂p ∂ρ s , c 2 i ≡ ∂p ∂ρ T (2.10)
An application of general thermodynamic methods [16] show that the quantities introduced in (2.8), (2.9), and (2.10), are not independent. They are related by the following thermodynamic identities

γ -1 = T β 2 c 2 a c p , c 2 a = γc 2 i (2.11)
where γ ≡ c p /c v is the specific heat ratio.

With v, the Euler's fluid velocity, σ the stress tensor, and q the heat flux, the equations of mass conservation, momentum conservation, and energy conservation, are expressed by

∂ρ ∂t + ∇ • (ρv) = 0 (2.12a) ∂(ρv i ) ∂t + ∂ j (ρv i v j -σ ij ) = 0 (2.12b) ∂(ρε) ∂t + ∇ • (ρεv + q) = σ ij ∂ j v i (2.12c)
TEqs.(2.7) and (2.8) results in the following relation

dε = -pd(1/ρ) + T ds (2.13)
As such, Eq.(2.12c) may be put in the following equivalent form of a balance law for entropy

∂(ρs) ∂t + ∇ • ρsv + 1 T q = 1 T (σ ij + pδ ij ) ∂ j v i + q • ∇ 1 T (2.14)
In the right-hand side, one sees the density of local entropy sources in the fluid, which is required to be positive by the second law. Now substituting in Eqs.(2.12a), (2.12b) and (2.12c), the total fields

ρ(t, x) = ρ 0 + ρ ′ (t, x) (2.15a) σ ij (t, x) = -pδ ij + σ ′ ij = -(P 0 + p ′ (t, x))δ ij + σ ′ ij (t, x) (2.15b) s(t, x) = s 0 + s ′ (t, x) (2.15c) T (t, x) = T 0 + τ (t, x) (2.15d)
where ρ 0 , P 0 , s 0 , T 0 represent the constant thermodynamic equilibrium values, and p ′ the thermodynamic excess pressure p ′ = p(T, 1/ρ) -p(T 0 , 1/ρ 0 ) = p(T, 1/ρ) -P 0 , the following linearized equations governing the small amplitudes perturbations are immediately obtained

∂ρ ′ ∂t + ρ 0 ∇ • v = 0 (2.16) ρ 0 ∂v i ∂t = -∂ i p ′ + ∂ j σ ′ ij (2.17) ρ 0 ∂s ′ ∂t = - 1 T 0 ∇ • q (2.18)
To close the system of equations, there remain to precise the constitutive laws which give the deviatoric stresses σ ′ ij and the heat flux q in terms of other variables. In the Navier-Stokes-Fourier theory used in this paper, it is assumed that the σ ′ ij are purely viscous; the Maxwell's stress terms [17] appearing in nonisothermal fluids and required by the kinetic theory of gases are neglected. In addition, molecular relaxation phenomena will not be considered. If necessary they may be incorporated as done in e.g. [18]. Within these simplifications, the constitutive equations are written, in the following Newton-Stokes form

σ ′ ij = η ∂ j v i + ∂ i v j - 2 3 δ ij ∇.v + ζδ ij ∇.v (2.19)
and Fourier form

q i = -κ∂ i T (2.20)
The coefficients of thermal conductivity κ, first and second viscosity η and ζ, are constants to be evaluated in the ambient state (P 0 , T 0 ). The second law of thermodynamics results in the inequalities

κ ≥ 0, η ≥ 0, ζ ≥ 0 (2.21)
The precise values of these fluid constitutive parameters are difficult to obtain by molecular theories; they are generally best found by experiment.

Inserting the constitutive equations (2.19-2.20) in the 5 balance equations (2.12a), (2.12b) and (2.14), we find for the system of 5 linearized conservation equations of mass, momentum, and energy

∂ρ ′ ∂t + ρ 0 ∇ • v = 0 (2.22a) ρ 0 ∂v i ∂t = -∂ i p ′ + ∂ j η ∂ j v i + ∂ i v j - 2 3 δ ij ∇.v + ζδ ij ∇.v (2.22b) ρ 0 ∂s ′ ∂t = κ T 0 ∇ 2 τ (2.22c)
In what follows, we keep using the 6 variables velocity, condensation b ≡ ρ ′ /ρ 0 , excess pressure, and excess temperature. To obtain a closed system of equations on these 6 variables we employ the equation of state ρ = ρ(p, T ), whose linearized version gives

b = ρ ′ /ρ 0 = γ ρc 2 a 0 p ′ -β 0 τ (2.
χ 0 ≡ 1/ρc 2 a 0
(2.25)

our complete system of six linearized viscothermal equations on the six variables v, b, p ′ , τ , is finally written as follows

ρ 0 ∂v ∂t = -∇p ′ + η∇ 2 v + ζ + η 3 ∇(∇ • v) (2.26a) ∂b ∂t + ∇ • v = 0 (2.26b) γχ 0 p ′ = b + β 0 τ (2.26c) ρ 0 c p ∂τ ∂t = β 0 T 0 ∂p ′ ∂t + κ∇ 2 τ (2.

Propagation in the circular tube

The solutions of the viscothermal fluid equations were given by Kirchhoff, for the case of axisymmetric wave propagation in a circular tube filled with ideal gas, taking into account viscous losses and thermal exchanges [11]. Kirchhoff made the simplification that the solid walls remain at ambient temperature due to the large heat capacity and conduction coefficient of the solid compared to the fluid. This is in general a well-verified simplification which is used in the following as well. Here, we derive the axisymmetric 'Kirchhoff-Langevin's' solutions of the above equations (2.26), more general than Kirchhoff's as they are written for a fluid having arbitrary equation of state (thus β 0 is different from the ideal gas value 1/T 0 ) and nonzero value of the second viscosity.

We notice that the reason to study only the axisymmetric solutions, is that we later intend to use the macroscopic nonlocal theory to derive anew the macroscopic characteristics of these solutions -wavenumbers and impedances. But this macroscopic theory, by definition, concerns the response to a 'macroscopic stirring' whose variations over the transverse cross-section of the pores are smoothed out and must not be considered. Therefore in circular pores, by symmetry, the fields meaningful to consider in the framework of the macroscopic theory are microscopically axisymmetric. Following Rayleigh's presentation of Kirchhoff's theory [19], we first substitute the state Eq.(2.26c) in Fourier's equation (2.26d) and use the thermodynamic identity (2.11) and definition (2.25) to obtain the following alternative form of (2.26d)

∂τ ∂t = γ -1 β 0 ∂b ∂t + κ ρ 0 c v ∇ 2 τ (2.27)
The equations are simplified, using the variable

τ ′ = β 0 τ γ -1 (2.28)
With this variable, Eqs.(2.26c) and (2.27) become

p ρ 0 = c 2 i b + (c 2 a -c 2 i )τ ′ (2.29)
and

∂τ ′ ∂t = ∂b ∂t + κ ρ 0 c v ∇ 2 τ ′ (2.30)
Then, assuming that the variables v, b, p, τ ′ , are varying with time as e -iωt , the equations (2.26) yield

-ρ 0 iωv = -∇p + η∇ 2 v + ζ + η 3 iω∇b (2.31a) -iωb + ∇ • v = 0 (2.31b) p ρ 0 = c 2 i b + (c 2 a -c 2 i )τ ′ (2.31c) -iωτ ′ = -iωb + κ ρ 0 c v ∇ 2 θ ′ (2.31d)
Eliminating the pressure and condensation, give rise to the following velocitytemperature equations

-iωv - η ρ 0 ∇ 2 v = -∇X (2.32a) X = c 2 a - η 3 + ζ ρ 0 iω τ ′ + κ ρ 0 c v iω c 2 i - η 3 + ζ ρ 0 iω ∇ 2 τ ′ (2.32b) ∇ • v -iωτ ′ - κ ρ 0 c v ∇ 2 τ ′ = 0 (2.32c)
Elimination of the velocity by taking the divergence of (2.32a) and using (2.31b) results in the temperature equation

ω 2 τ ′ + c 2 a -iω κ ρ 0 c v + 4η 3 + ζ ρ 0 ∇ 2 τ ′ + κ ρ 0 c v iω c 2 i -iω 4η 3 + ζ ρ 0 ∇ 2 ∇ 2 τ ′ = 0 (2.33)
Let λ 1 and λ 2 be the two, small and large, solutions of the associated Kirchhoff-Langevin's characteristic equation

ω 2 + c 2 a -iω κ ρ 0 c v + 4η 3 + ζ ρ 0 λ + κ ρ 0 c v iω c 2 i -iω 4η 3 + ζ ρ 0 λ 2 = 0 (2.34)
The small solution -mainly real -describes propagating acoustic waves with small bulk absorption, the large solution -purely imaginary -highly damped diffusive entropic waves.

The field τ ′ solution to (2.33) will have the form

τ ′ = A 1 ϕ 1 + A 2 ϕ 2 (2.35)
with functions ϕ 1 and ϕ 1 verifying

∇ 2 ϕ 1 = λ 1 ϕ 1 , ∇ 2 ϕ 2 = λ 2 ϕ 2 (2.36)
The velocity v will write

v = v ′ + B 1 ∇ϕ 1 + B 2 ∇ϕ 2 (2.37)
with v ′ the vortical part, such that

∇ 2 v ′ = -iωρ 0 η v ′ , ∇ • v ′ = 0 (2.38)
The relation between coefficients B and A follows from (2.32c)

B 1,2 = κ ρ 0 c v + iω λ 1,2 A 1,2 (2.39)
In summary, the viscothermal fields are generally decomposed as

v = v ′ + κ ρ 0 c v + iω λ 1 A 1 ∇ϕ 1 + κ ρ 0 c v + iω λ 2 A 2 ∇ϕ 2 (2.40a) b = 1 + κ ρ 0 c v iω λ 1 A 1 ϕ 1 + 1 + κ ρ 0 c v iω λ 2 A 2 ϕ 2 (2.40b) p ρ 0 = c 2 a + c 2 i κ ρ 0 c v iω λ 1 A 1 ϕ 1 + c 2 a + c 2 i κ ρ 0 c v iω λ 2 A 2 ϕ 2 (2.40c) τ ′ = A 1 ϕ 1 + A 2 ϕ 2 (2.40d)
Concerning the application to axisymmetric fields propagating in the right-going x direction in a tube of circular cross-section, we wish to determine normal modes as functions of x, proportional to e +ik l x , where the k l s, ℑ(k l ) > 0, l = 1, 2, ..., are complex constants to be specified. In what follows for convenience the index l, labeling the different axisymmetric modes solutions, will be omitted.

For these modes, the operator ∇ can be replaced by ike x + ∂ ∂r e r (e r representing the radial unit vector) and the operator ∇ 2 by ∂ 2 ∂r 2 + ∂ r∂r -k 2 , and the different fields a(t, x) by their amplitudes a such that a(t, x) = a(r)e -iωt+ikx . There follows that the corresponding ϕ 1 and ϕ 2 will be described by Bessel functions

ϕ 1,2 = J 0 r -λ 1,2 -k 2 (2.41)
Writing the vortical velocity v ′ in the form v ′ = u ′ e x + q ′ e r , with axial and radial amplitudes u ′ and q ′ independent of azimuthal angle ϑ, it is easy to see that Eqs.(2.38) imply u ′ is the solution to

∂ 2 ∂r 2 + ∂ r∂r u ′ = -iω ν + k 2 u ′ (2.42)
where ν ≡ η/ρ 0 is the kinematic viscosity and q ′ is determined by the relation

q ′ = -ik -iω ν + k 2 ∂u ′ ∂r (2.43)
As a result, u ′ and q ′ will be written as

u ′ = Aϕ, q ′ = A -ik -iω ν + k 2 ∂ϕ ∂r (2.44)
where ϕ is the Bessel function

ϕ = J 0 r iω ν -k 2 (2.45)
Finally, writing the total velocity in the form of v = ue x + qe r , u and q are obtained as

u = Aϕ + ik κ ρ 0 c v + iω λ 1 A 1 ϕ 1 + ik κ ρ 0 c v + iω λ 2 A 2 ϕ 2 (2.46a) q = A -ik -iω ν + k 2 ∂ϕ ∂r + κ ρ 0 c v + iω λ 1 A 1 ∂ϕ 1 ∂r + κ ρ 0 c v + iω λ 2 A 2 ∂ϕ 2 ∂r (2.46b)
The tube being assumed sufficiently inert thermally to remain at ambient temperature, no temperature-jump occurs on the tube wall r = R. On the other hand no-slip condition is applied on the wall. Thus, (2.46a-2.46b) and (2.40d) vanish on the fluid-solid boundaries. These three homogeneous equations have non vanishing solutions only if their determinant is zero, which consequently yields the following Kirchhoff's dispersion equation

κ ρ 0 c v + iω λ 1 1 ϕ 1w ∂ϕ 1 ∂r w - κ ρ 0 c v + iω λ 2 1 ϕ 2w ∂ϕ 2 ∂r w - k 2 -iω ν + k 2 iω λ 1 - iω λ 2 1 ϕ w ∂ϕ ∂r w = 0 (2.47)
referred here as Kirchhoff-Langevin equation to remind its validity for a general viscothermal fluid. The index w indicates that the functions and derivatives are evaluated at the tube wall r w = R.

Eq.(2.47) has a series of discrete complex solutions k l , ℑ(k l ) > 0, l = 1, 2, ..., which can be sorted by convention in ascending manner relating the values of ℑ(k l ) such that ℑ(k 1 ) < ℑ(k 2 ) < ... . To determine the solutions k l , the Newton-Raphson root-finding method may be employed, with initial values k 0l taken as if the dissipation effects are neglected. In this ideal lossless case, Eq.(2.47) will be considered in the limit η, ζ, κ → 0, and therefore k 0l will be the solutions to

J 1 R ω 2 c 2 a -k 2 = 0 (2.48) that means k 2 0l = ω 2 c 2 a - x 2 l R 2 (2.49)
where x l ≥ 0 are the successive zeros of the function J 1 (x). As J 1 (0) = 0 i.e.

x 1 = 0, the first solution is always k 01 = ω/c a , which corresponds to the plane wave mode. At a given frequency, there are one (plane wave mode k 01 ) or more real positive solutions describing right-going propagating waves, and an infinite discrete set of purely imaginary solutions with ℑ(k 0l ) > 0, describing evanescent waves which attenuate as e -ℑ(k 0l )x along positive x-axis.

In general, a few Newton-Raphson iterations suffice to make these starting losseless purely real or purely imaginary solutions k 0l converge towards the complete complex solutions k l . In this process, the solution with positive imaginary part is retained, as we consider waves propagating in the direction +x which can be created by a source placed on the left. The condition that the imaginary part is positive, automatically fixes the sign of the real part. We note that, for a given solution l and on account of the two independent conditions expressing the vanishing of τ ′ and u at the tube wall

r w = R A 1 ϕ 1w + A 2 ϕ 2w = 0 (2.50a) Aϕ w + ik κ ρ 0 c v + iω λ 1 A 1 ϕ 1w + ik κ ρ 0 c v + iω λ 2 A 2 ϕ 2w = 0 (2.50b)
the solution may be fixed in terms of only one arbitrary amplitude (related to the arbitrary sound pressure level). We denote it by A and write

A = Aik iω λ 1 - iω λ 2 ϕ 1w ϕ 2w , A 1 = -Aϕ w ϕ 2w , A 2 = Aϕ w ϕ 1w (2.51)
In summary, with this notation the different fields write as follows

u ikA = iω λ 1 - iω λ 2 ϕ 1w ϕ 2w ϕ - κ ρ 0 c v + iω λ 1 ϕ w ϕ 2w ϕ 1 + κ ρ 0 c v + iω λ 2 ϕ w ϕ 1w ϕ 2 (2.52a) q A = iω λ 1 - iω λ 2 ϕ 1w ϕ 2w k 2 -iω ν + k 2 ∂ϕ ∂r - κ ρ 0 c v + iω λ 1 ϕ w ϕ 2w ∂ϕ 1 ∂r + κ ρ 0 c v + iω λ 2 ϕ w ϕ 1w ∂ϕ 2 ∂r (2.52b) b A = -1 + κ ρ 0 c v iω λ 1 ϕ w ϕ 2w ϕ 1 + 1 + κ ρ 0 c v iω λ 2 ϕ w ϕ 1w ϕ 2 (2.52c) p -Aρ 0 = c 2 a + c 2 i κ ρ 0 c v iω λ 1 ϕ w ϕ 2w ϕ 1 + c 2 a + c 2 i κ ρ 0 c v iω λ 2 ϕ w ϕ 1w ϕ 2 (2.52d) τ ′ A = -ϕ w ϕ 2w ϕ 1 + ϕ w ϕ 1w ϕ 2 (2.52e)
The key step to determine the above field patterns is to specify the wavenumber k l . Suitable averaging of the above fields will then allow to compute the quantities making sense in the forthcoming macroscopic theory.

Let us denote by a bracket f the average of a field f , performed over the crosssection of the tube

f = 1 πR 2 R 0 dr 2π 0 rdϑf (2.53)
For later use in the macroscopic theory, we introduce the following notion of macroscopic mean pressure

H pv = H v (2.54)
Thermodynamic-relating reasons to define in such a way the macroscopic pressure were given in [1]. The notation H comes from the affinity of this concept with that of Maxwell magnetic field H. In accordance with this general definition, a characteristic complex impedance factor Z l can be defined for a given mode, by setting

H = Z l u , such that Z l = pu u 2 (2.55)
where the averages are computed knowing the wavenumber k l .

It may be noted that, by knowing the complex wavenumber k l (ω) and impedance Z l (ω) for a given radial mode solution l, an equivalent-fluid complex density ρ l (ω) and bulk modulus χ -1 l (ω) might be defined through setting relations having the usual form (see e.g. [?])

k l = ω/c l = ω √ ρ l χ l and Z l = ρ l c l = ρ l χ -1 l , i.e. ρ l (ω) = k l ω Z l = k l ω pu u 2 , χ -1 l (ω) = ω k l Z l = ω k l pu u 2 (2.56)
Indeed, for a given mode l, this means that we may write classical macroscopic equivalent-fluid equations of motion

ρ l (ω) ∂ ∂t u = - ∂ ∂x H, χ l (ω) ∂ ∂t H = - ∂ ∂x u
To verify this, consider one wave e -iωt+ik l x in the +x direction, we have

-ρ l (ω)iω u = -ik l H, -χ l (ω)iωH = -ik l u (2.57)
By introducing the definitions (2.55) the above gives the relations (2.56).

There are known formulae giving the average of Bessel functions ϕ, ϕ 1 and ϕ 2 and their products, in terms of other Bessel functions. Thus, when k l is known, the above Kirchhoff-Langevin's impedance factors Z l (2.55), and densities and bulk modulii ρ l (ω) and χ -1 l (ω) (2.56) are all known in closed form.

General nonlocal theory applied to sound propagation in a tube with circular cross-section

Applying the general nonlocal theory of sound propagation presented in [1], to the case of sound propagation in a tube of circular cross-section, the macroscopic averaging operation in the sense of [1] evidently becomes the cross-section average (2.53). In fact, the fields relating to the harmonic action-response problems [1] have the form a(ω, k, r)e -iωt+ikx , with amplitudes a(ω, k, r) independent of x, and thereby, needed to be averaged only over a cross-section.

Linearized macroscopic equations in Maxwellian nonlocal theory

Then introducing the macroscopic variables V = V • e x = u and B = b , where u is the axial velocity, the nonlocal theory (chapter 1) predicts that wave

D(t, x) = t -∞ dt ′ dx ′ ρ(t -t ′ , x -x ′ )V (t ′ , x ′ ) (2.60) H(t, x) = t -∞ dt ′ dx ′ χ -1 (t -t ′ , x -x ′ )B(t ′ , x ′ ) (2.61)
where ρ(t, x) and χ -1 (t, x) are constitutive kernel functions independent of temporal and spatial variations of the fields. They are determined only by the fluid constants and the microgeometry of the porous medium, i.e. here, the tube radius R. The integrations over t ′ determine the so-called temporal dispersion effects and the integrations over x ′ determine the so-called spatial dispersion effects [20,21]. The upscaling recipes, seen in chapter 1, lead to specify the Fourier coefficients ρ(ω, k) and χ -1 (ω, k) of these constitutive functions, and will be described in the next two sections.

Determination of the nonlocal density ρ(ω, k)

To compute ρ(ω, k) we first consider the response of the fluid subjected to the action of an external stirring force-per-unit-volume f , which derives from a fictitious harmonic pressure waveform inserted in Navier-Stokes equation. Thus we consider solving the action-response problem (see Eqs.(1.45-1.47))

∂v ∂t + 1 ρ 0 ∇p = ν∇ 2 v - ( η 3 + ζ) ρ 0 ∇ ∂b ∂t + 1 ρ 0 f (2.62a) ∂b ∂t + ∇ • v = 0 (2.62b) γχ 0 p = b + β 0 τ (2.62c) ∂τ ∂t = β 0 T 0 ρ 0 c p ∂p ∂t + κ ρ 0 c p ∇ 2 τ (2.62d)
for r < R, and

v = 0 (2.63a) τ = 0 (2.63b)
at r = R, with the deriving force given by f = -∇P = -ike x P 0 e -iωt+ikx (2.64)

With calculations entirely similar to those which have been done before, the equations (2.32) become

-iωv -ν∇ 2 v = -∇X + 1 ρ 0 f (2.65a) X = c 2 a - η 3 + ζ ρ 0 iω τ ′ + κ ρ 0 c v iω c 2 i - η 3 + ζ ρ 0 iω ∇ 2 τ ′ (2.65b) ∇ • v -iωτ ′ - κ ρ 0 c v ∇ 2 τ ′ = 0 (2.65c)
and the temperature equation (2.33) becomes

-ω 2 τ ′ -c 2 a -iω κ ρ 0 c v + 4η 3 -ζ ρ 0 ∇ 2 τ ′ - κ ρ 0 c v iω c 2 i -iω 4η 3 + ζ ρ 0 ∇ 2 ∇ 2 τ ′ + k 2 ρ 0 P = 0 (2.66)
A particular solution of the above equation is

τ ′ p = C k 2 ρ 0 P (2.67) with C = ω 2 -c 2 a -iω κ ρ 0 c v + 4η 3 + ζ ρ 0 k 2 + κ ρ 0 c v iω c 2 i -iω 4η 3 + ζ ρ 0 k 4 -1
(2.68)

The general solution of Eq.(2.66) is this particular solution added to the general solution (2.35) of the homogeneous equation (2.33)

τ ′ = τ ′ 0 + τ ′ p = A 1 ϕ 1 + A 2 ϕ 2 + C k 2 ρ 0 P (2.69)

Application of the general nonlocal theory 51

Similarly, the general expression of the velocity will be

v = ue x + qe r = v 0 + v p (2.70)
where v 0 is written as in (2.40a), and v p = u p e x is the particular solution with u p determined by (2.65c), and τ ′ = τ ′ p , i.e., iku p = (iω -κk 2 /ρ 0 c v )τ ′ p . Only the x component u will be required to compute ρ(ω, k). However, the radial component q needs also to be written as it is involved in the boundary conditions, by means of which the amplitudes A, A 1 , A 2 are finally fixed. Both components write accordingly

u = Aϕ + ik κ ρ 0 c v + iω λ 1 A 1 ϕ 1 + ik κ ρ 0 c v + iω λ 2 A 2 ϕ 2 + -iω + κ ρ 0 c v k 2 C ik ρ 0 P (2.71a) q = A -ik -iω ν + k 2 ∂ϕ ∂r + κ ρ 0 c v + iω λ 1 A 1 ∂ϕ 1 ∂r + κ ρ 0 c v + iω λ 2 A 2 ∂ϕ 2 ∂r (2.71b)
Now, we seek the excess pressure solution as the last required quantity. It has the general form

p = p 0 + p p (2.72)
where p 0 is given by (2.40c), and the particular solution p p is determined by (2.26c) with τ ′ = τ ′ p , b = b p , and iωb p = iku p . Thus

p ρ 0 = c 2 a + c 2 i κ ρ 0 c v iω λ 1 A 1 ϕ 1 + c 2 a + c 2 i κ ρ 0 c v iω λ 2 A 2 ϕ 2 + c 2 a -c 2 i κ ρ 0 c v iω k 2 C k 2 ρ 0 P (2.73)
The boundary conditions imply that the three quantities excess temperature and the two components of velocity (2.69), (2.71a) and (2.71b), should vanish at the tube wall r = r w = R. This yields a linear system whose solution uniquely determines the three response amplitudes A, A A direct verification of this conjectured upscaling procedure would be to see whether or not, when the amplitudes A and P 0 are adjusted so that the gradient -∇H = -ik l e x H in section 2.2.2 is the same as the gradient -∇(P + P) = -ik l e x (P + P) in this section, the averages of the velocities v appearing respectively in Eqs.(2.26) and (2.62) turn out to be exactly the same. This macroscopic coincidence of the two mean velocities obtained in two different problems, is highly nontrivial even in the present simplest case of straight duct geometry. In principle, the two problems have fundamentally different nature; one is an eigenvalue problem and the other is an action-response problem. In the two problems, the corresponding sets of microscopic field patterns are not the same; but after averaging, the two mean velocities divided by the corresponding two gradients are conjectured to be exactly the same.

A comparison between the first Eq.(2.57) and Eq.(2.75) shows that this matching is expressed in explicit equivalent form by the following equation

ρ l (ω) = ρ(ω, k l ) (2.76)
It is in this last convenient form, which must be valid for all different modes, that the validity of the upscaling procedure for ρ(ω, k) will be directly checked.

We may name 'Maxwell-Kirchhoff's' the nonlocal density function (2.75). There are known formulae giving the average of Bessel functions ϕ, ϕ 1 and ϕ 2 and their products, in terms of other Bessel functions, which allow to write the nonlocal density function (2.75) in closed form. This function is a complicate ratio of sums containing many terms, each involving the product of 6 Bessel functions, multiplied by factors involving ω, k 2 , λ 1 , λ 2 . To save time, instead of seeking its most compact final expression, we have made a direct Matlab programming of it with ω and k as input arguments.

Determination of the nonlocal bulk modulus χ -1 (ω, k)

The same type of calculations seen before to obtain ρ(ω, k), can be performed to directly compute χ -1 (ω, k). Here, we need first consider the response of the fluid subjected to the action of an external stirring rate of heat supply per unit volume and unit time Q, which derives from a fictitious harmonic pressure waveform inserted in Fourier equation. Thus we consider solving the actionresponse problem (see Eqs.(1.56-1.58))

∂v ∂t + 1 ρ 0 ∇p = ν∇ 2 v - ( η 3 + ζ) ρ 0 ∇ ∂b ∂t (2.77a) ∂b ∂t + ∇ • v = 0 (2.77b) γχ 0 p = b + β 0 τ (2.77c) ∂τ ∂t = β 0 T 0 ρ 0 c p ∂p ∂t + κ ρ 0 c p ∇ 2 τ + 1 ρ 0 c p Q (2.77d)
for r < R, and

v = 0 (2.78a) τ = 0 (2.78b)
at r = R, with the stirring rate of heat supply given by

Q = β 0 T 0 ∂P ∂t = -iωβ 0 T 0 P 0 e -iωt+ikx (2.79)
Through similar calculations as seen before, Eqs.(2.32) now become

-iωv -ν∇ 2 v = -∇X (2.80a) X = c 2 a - η 3 + ζ ρ 0 iω τ ′ + κ ρ 0 c v iω c 2 i - η 3 + ζ ρ 0 iω ∇ 2 τ ′ (2.80b) ∇ • v -iωτ ′ - κ ρ 0 c v ∇ 2 τ ′ = γχ 0 ∂P ∂t (2.80c)
and the temperature equation (2.33) becomes

-ω 2 τ ′ -c 2 a -iω κ ρ 0 c v + 4η 3 -ζ ρ 0 ∇ 2 τ ′ - κ ρ 0 c v iω c 2 i -iω 4η 3 + ζ ρ 0 ∇ 2 ∇ 2 τ ′ -(iω -νk 2 )γχ 0 iωP = 0 (2.81)
A particular solution to the above equation will be

τ ′ p = -C(iω -νk 2 )γχ 0 iωP (2.82)
with the same constant C as before.

The general excess temperature solution is this particular solution added to the general solution (2.35) of the homogeneous equation (2.33)

τ ′ = τ ′ 0 + τ ′ p = A 1 ϕ 1 + A 2 ϕ 2 -C(iω -νk 2 )γχ 0 iωP (2.83)
Likewise, the general velocity solution is in the form

v = ue x + qe r = v 0 + v p (2.84)
where v 0 writes as in (2.40a), and v p = u p e x is the particular solution with u p determined by (2.80c) and τ ′ = τ ′ p , i.e., iku p = (iω -κk 2 /ρ 0 c v )τ ′ p -iωγχ 0 P. We obtain for the two components of the velocity

u = Aϕ + ik κ ρ 0 c v + iω λ 1 A 1 ϕ 1 + ik κ ρ 0 c v + iω λ 2 A 2 ϕ 2 + -iω + κ ρ 0 c v k 2 C(iω -νk 2 ) -1 γχ 0 ω k P (2.85a) q = A -ik -iω ν + k 2 ∂ϕ ∂r + κ ρ 0 c v + iω λ 1 A 1 ∂ϕ 1 ∂r + κ ρ 0 c v + iω λ 2 A 2 ∂ϕ 2 ∂r (2.85b)
The general pressure solution, similarly is written as

p = p 0 + p p (2.86)
where p 0 writes as in (2.40c) and p p is determined by (2.26c) with τ ′ = τ ′ p and b = b p , iωb p = iku p . We will have

p ρ 0 = c 2 a + c 2 i κ ρ 0 c v iω λ 1 A 1 ϕ 1 + c 2 a + c 2 i κ ρ 0 c v iω λ 2 A 2 ϕ 2 -C(iω -νk 2 )iω c 2 a -c 2 i κ ρ 0 c v iω k 2 γχ 0 P -c 2 i γχ 0 P (2.87)
The general condensation solution will now also be required. It may be written from (2.31c) and the expressions (2.83) and (2.87), which yield

b = 1 + κ ρ 0 c v iω λ 1 A 1 ϕ 1 + 1 + κ ρ 0 c v iω λ 2 A 2 ϕ 2 -C(iω -νk 2 )iω 1 - κ ρ 0 c v iω k 2 γχ 0 P -γχ 0 P (2.88)
As before, the boundary conditions imply that the three quantities excess temperature and the two components of velocity (2.83), (2.85a) and (2.85b) should vanish at the tube wall r = r w = R. This yields a linear system whose solution uniquely determines the three response amplitudes A, A

χ -1 (ω, k) = P(ω, k) + P 0 b(ω, k, r) + γχ 0 P 0 (2.90)
A direct verification of this, would be to see whether or not, when the amplitudes A and P 0 are adjusted so that the amplitude H in section 2.2 is the same as the amplitude P + P 0 in this section, then, the average of the condensation b appearing in Eqs.(2.26) and the average of the condensation b + γχ 0 P, where b is the quantity appearing in (2.77), turn out to be also exactly the same. Again, this macroscopic accordance is by no means trivial, even in the present simplest case of straight duct geometry. A comparison between the second equation of (2.57) and Eq.(2.90) shows that this matching is also expressed by the following equation which must be valid for all different modes

χ -1 l (ω) = χ -1 (ω, k l ) (2.91)
It is in this last convenient form, that the validity of the upscaling procedure for χ -1 (ω, k) will be directly checked.

As for the nonlocal density, we may name 'Maxwell-Kirchhoff's' the nonlocal bulk modulus function (2.90). As before, this Maxwell-Kirchhoff's nonlocal bulk modulus is known in closed form as a complicate ratio of sums containing many terms, each involving the product of 6 Bessels. Again to save time, a direct programming of the function χ(ω, k) was made, with ω and k the input arguments.

Dispersion equation, wavenumbers, and impedances

The above nonlocal theory predicts that, at a given frequency, normal mode solutions with averaged fields varying as e -iωt+ikx will exist, for k solution to the dispersion equation

ρ(ω, k)χ(ω, k)ω 2 = k 2 (2.92)
Indeed, this equation comes from Eqs.(2.58-2.61), making use of macroscopic fields having the dependencies e -iωt+ikx . Hereafter we refer to this dispersion equation as to the nonlocal Maxwell-Kirchhoff's dispersion equation.

For the proposed nonlocal theory to be correct, Maxwell-Kirchhoff's dispersion equation (2.92) must be mathematically equivalent to the original Kirchhoff-Langevin's dispersion Eq.(2.47). In particular, both equations must have the same set of solutions k l at given ω. Moreover, the macroscopic impedances of the corresponding modes must also be the same. Recall that for any field freely propagating in the tube, the nonlocal theory definition of the macroscopic field H(t, x) originates from the fundamental 'Umov-Poynting' identification

p(t, x, r)u(t, x, r) = H(t, x) u(t, x, r) (2.93)
Since we have employed this expression to define 'Kirchhoff-Langevin's' macroscopic impedance factors Z l (2.55), Kirchhoff-Langevin's and Maxwell-Kirchhoff wavenumbers and macroscopic impedances match automaticaly, provided the following aforementioned relations

ρ l (ω) = k l ω Z l = ρ(ω, k l ), χ -1 l (ω) = ω k l Z l = χ -1 (ω, k l ) (2.94)
The calculations to be performed in order to show the mathematical equivalence between (2.47) and (2.92) appear very tedious, because of the large number of terms to be collected and rearranged to express the mean term pu .

In what follows, to check the validity of the theory, the easier way of direct numerical checking of relations such as (2.94) has been employed. The problems relating to the precision in the Matlab computations would limit the number of modes that can be valuably checked. The results which will be presented next, however, clearly validate the theory.

Check on the nonlocal theory

In order to validate the theory, we detail three significant different cases, representative of three main types of duct wave regimes, referred to the 'narrow' tube, 'wide' tube, and 'very wide' tube, respectively in acoustic literature [22].

In section 2.4.1 we will consider the case of low frequencies or 'narrow' tubes.

In this low frequency range, the viscous skin depth δ = (2ν/ω) 1/2 and thermal skin depth -having the same order for air -are greater than R. Calculations are performed at frequency f = 100Hz for a tube of radius R = 10 -4 m. This is the narrow tube configuration considered in [23]. The value R = 10 -4 m is typical for the pore size dimensions found in ordinary porous materials used in noise control applications [15] -such as pore-size parameter Λ of dynamically connected pores in [3]. With a viscous skin depth equal to two times the radius, the fundamental plane-wave like mode is mostly diffusive and the higher order modes are highly attenuated.

Certainly due to insufficient accuracy of Matlab Bessel's functions, only the first mode appears to be numerically very well characterized. However, it provides a first check of the theory: Kirchhoff-Langevin's and Maxwell-Kirchhoff's results for this mode are found to be identically the same, up to the numerical accuracy.

In section 2.4.2 we consider the case of high frequencies or 'wide' tubes. In this frequency range, the viscous skin depth and thermal skin depth become significantly smaller than R. The calculations are done at frequency f = 10kHz for a tube of radius R = 10 -3 m. This is the 'wide' tube configuration considered in [23]. With the radius now being more than 50 times the viscous skin depth, the fundamental plane-wave like mode is a well-propagating mode, whose macroscopic characteristics may be followed with great numerical accuracy by the Matlab nonlocal theory computations. Several higher order modes are also successfully described by the present numerical computations, whether they are below or above the cutoff frequency. Again, the results provide unequivocal validation of the proposed nonlocal theory.

Finally, in section 2.4.3, taking a tube radius of 1cm, and a frequency f = 500kHz, we consider the case of 'very wide' tubes. Here, the fundamental least attenuated mode is no longer plane-like as predicted by Zwikker and Kosten theory. Sound energy tends to concentrate near the walls. Once again Kirchhoff-Langevin's results are exactly reproduced by the nonlocal Maxwell-Kirchhoff's theory and provide unequivocal validation of the theory. In all foregoing calculations the parameters of the air are set to the values shown in Table 2.1. Given the radius R and the frequency f , we proceed as follows to evaluate the different quantities. To evaluate the Zwikker and Kosten quantities, the formulae reported in Appendix A are used. To evaluate the Kirchhoff-Langevin's and Maxwell-Kirchhoff's wavenumbers k KL and k M K , a Newton-Raphson scheme is used to solve, Eq.(2.47) and Eq.(2.92), respectively. In the Kirchhoff case, as we dispose of the explicit analytical expression of the function F (ω, k) = 0, we have an explicit analytical expression for its derivative with respect to k. In the Maxwell case, the expressions are too lengthy, thereby we make use of a numerical derivative. This is carefully done by evaluating the function F at two close values of the wavenumber k(1 ± ǫ d /2 e iθ ), with ǫ d a very small and adjustable parameter (e.g. ǫ d = 10 -9 ), and then averaging over several random orientations θ between 0 and 2π. The results which are reported below have been shown to be insensitive to the variation of ǫ d . Our stopping condition of the Newton scheme is that the relative error between two successive evaluations of k KL or k M K should be less than a very small fixed value ǫ s (e.g. ǫ s = 10 -12 ). The Newton scheme is more stable and converges in fewer iterations for Kirchhoff-Langevin's equation (2.47) than for Maxwell-Kirchhoff's equation (2.92).

ρ 0 T 0 c 0 η ζ κ χ 0 c p γ (kg/m 3 ) (K) (m/s) (kg ms -1 ) (kg ms -1 ) (W m -1 K -1 ) (P a -1 ) (J kg -1 K -1 )
Given a wavenumber k = k KL = k M K , the field patterns in the Kirchhoff-Langevin source-free propagation problem, and Maxwell-Kirchhoff actionresponse problems, can be explicitly written.

For the Kirchhoff-Langevin's quantities, the impedance Z KL has been computed first, using the explicit expression (2.55). The density ρ KL and bulk modulus χ -1 KL were then obtained using the relations (2.56).

For the Maxwell-Kirchhoff quantities, the density and bulk modulus were computed first, using the direct programming of the relations (2.74-2.75) and

(2.89-2.90), and as input argument k, the mode wavenumber, either coming from the Newton solution of the Kirchhoff-Langevin's equation, or the Newton solution of the Maxwell-Kirchhoff's equation. As the two wavenumbers may not coincide exactly due to finite precision and existing inaccuracies in Matlab Bessel functions, this does not produce exactly the same values of density and bulk modulus. Corresponding end-values will be distinguished by using indexes M K KL or M K M K respectively. Corresponding impedances could then be computed from (2.5), and corresponding wavenumbers, by (2.4). If the theory is correct, we expect to see consistency between these different quantities, when no numerical problems arise; what will be shown below.

Narrow tubes: R = 10 -4 m, f = 100Hz

In order to distinguish between Zwikker and Kosten's, Kirchhoff

i χ -1 M K KL 9.962016409534367 × 10 4 -1.043531769003892 × 10 3 i χ -1 M K M K 9.962016409534391 × 10 4 -1.043531769003666 × 10 3 i ℜ(∆χ -1 /χ -1 ) < 10 -14 ℑ(∆χ -1 /χ -1 ) < 10 -14
Kirchhoff-Langevin's value k KL is relatively insensitive to the starting value. Convergence to the given solution is obtained by starting from the lossless case solution k = ω/c a , the Zwikker-Kosten solution, or the value k = 6 + 2i taken in [23]. On the contrary, Maxwell-Kirchhoff's value k M K is sensitive to the starting value. Not reported here, a meaningless unattenuated value k M K is found, when using as starting value the lossless-relating solution k = ω/c a .

The errors ∆f /f indicate the relative differences computed between Maxwell-Kirchhoff's values and Kirchhoff-Langevin's reference values. Their small values show that they are numerically insignificant. There is complete matching of wavenumbers, impedances, densities, and bulk modulus, for the first least attenuated mode. For the first higher order mode, the imaginary part of the wavenumber is already in the order 5. × 10 4 i. This leads to a large imaginary part in the complex arguments of the Bessel functions ϕ and ϕ 1,2 . The resulting loss of precision prevents making precise checks with Matlab.

2.4.2 Wide tubes: R = 10 -3 m, f = 10kHz

For the least attenuated plane wave mode, the values obtained of the wavenumbers, impedances, densities and bulk modulii are:

The deviations are not significant owing to the calculation precision. This is again a clear validation of the nonlocal upscaling procedures. The cutoff frequency of the first higher order axisymmetric mode is a little above 10kHz. While this mode is still very significantly attenuated, its macroscopic characteristics k, Z,

k ZW 1.877218171102030 × 10 2 + 3.047328259173055i k KL 1.877217761268940 × 10 2 + 3.050105788888088i k M K 1.877217761268940 × 10 2 + 3.050105788888080i k M K KL 1.877217761268940 × 10 2 + 3.050105788888080i k M K M K 1.
χ -1 ZW 1.379578791782648 × 10 5 -1.406078512972857 × 10 3 i χ -1 KL 1.379578235612869 × 10 5 -1.407920077798555 × 10 3 i χ -1 M K KL 1.379578235612869 × 10 5 -1.407920077798545 × 10 3 i χ -1 M K M K 1.379578235612869 × 10 5 -1.407920077798562 × 10 3 i ℜ(∆χ -1 /χ -1 ) < 10 -15 ℑ(∆χ -1 /χ -1 ) < 10 -15
nevertheless, are obtained with a precision which, once again, shows the exactness of the upscaling:

It may be noted that the negative real part of the bulk modulus is the type of behaviour described for metamaterials [10], the negative real part of wavenumber also being present and associated with negative group velocity.

× i χ -1 M K KL -3.235050549225333 × 10 8 + 2.885427853735547 × i χ -1 M K M K -3.235050516110349 × 10 8 + 2.885427853625859 × i ℜ(∆χ -1 /χ -1 )
< 10 -10 ℑ(∆χ -1 /χ -1 ) < 10 -10 2.4.3 Very wide tubes: R = 10 -2 m, f = 500kHz

Recall that in this new regime of the wave propagation the least attenuated mode is no longer a plane mode. It tends to concentrate near the walls. The values obtained of the wavenumbers, impedances, densities and bulk modulii are They show the exactness of the upscaling procedure in this regime, as well. Since spatial nonlocality plays an essential role here, there are considerable differences between Zwikker and Kosten values and the exact ones. Recall that the local Zwikker and Kosten theory assimilates the field H with mean pressure p and the mean pressure with the pressure itself. But here, the pressure is no longer a constant over the section, so that the local approach is largely in error.

i χ -1 M K M K 8.279327305835144 × 10 4 -5.269923491003189 × 10 5 i ℜ(∆χ -1 /χ -1 ) < 10 -12 ℑ(∆χ -1 /χ -1 ) < 10 -12

Conclusion

The exact matching between the macroscopic translation of Kirchhoff-Langevin's results and the results obtained on the basis of the new nonlocal theory proposed in [1], provides a clear validation of the nonlocal-relating upscaling procedures.

The important concept, which leads us in [1] to these exact homogenization upscaling procedures, is the 'Umov-Poynting-Heaviside' concept of 'acoustic part of energy current density'. Using an analogy with electromagnetics, we had to identify exactly this acoustic part with the quantity s = pv, where p is the thermodynamic pressure. This served in turn as a basis to define a macroscopic pressure field, through the macroscopic relation definition pv = H v . And this thermodynamic definition can be recognized finally a key to make use of the solutions of two simple action-response problems, in the appropriate way, which lead to the independent computation of the two nonlocal acoustical susceptibilities ρ and χ.

In forthcoming papers it will be shown that the proposed nonlocal Maxwellian theory, providing exact homogenization procedures, is valid also in the case of nontrivial geometries.

We believe that the present theory highlights the unsatisfactory thermophysical state of affairs in macroscopic electromagnetic theory, where a comparable definition of the macroscopic magnetic field H cannot yet be proposed, because of lacking thermodynamic variables allowing to express the concept of 'electromagnetic part of energy current density'. For this reason, no electromagnetic analogue of the present acoustic upscaling procedures can yet be proposed.

Chapter 3

Nonlocal theory of sound propagation in porous media; case of two-dimensional arrays of rigid cylinders

Introduction

A new nonlocal theory of sound propagation in rigid-framed porous media saturated with a viscothermal fluid has been recently proposed [1], which is considered to provide for the first time an exact homogenization procedure. By this theory we can treat, in particular, the media which are microscopically periodic, and macroscopically homogeneous, and the propagation is along a symmetry waveguide axis. A first successful test of this theory has been made in the simple case of cylindrical circular tubes filled with a viscothermal fluid (see chapter 2). It was found that the wavenumbers and impedances predicted coincide with those of the long-known Kirchhoff's full solution [2]. Here, we want to verify the validity of this new nonlocal theory in the case where the microgeometry of the porous medium is nontrivial, in the form of an unbounded two-dimensional square lattice of rigid cylinders permeated by a viscothermal fluid (see Fig. 3.1). This geometry allows a direct quasi-analytical calculation of the medium properties by a multiple-scattering approach taking into account viscous and thermal effects [3]. If, as guessed, the theory is exact, a matching will be observed between the multiple scattering predictions and the new theory predictions, independently of the frequency range considered.

To verify this, on the one hand, we will compare the complex frequency-dependent phase velocity associated to the least attenuated plane wave, predicted by a FEM implementation the new theory, with that of the corresponding least attenuated 68 3 Case of two-dimensional arrays of rigid cylinders Bloch mode, obtained by the quasi-exact multiple scattering method, and show that the two are in remarkable agreement. On the other hand, comparison with the single mode complex phase velocity obtained by the previously existing local theory, will show the domain of validity of the local description and its limit in terms of frequency.

By local theory we refer to space locality. Nonlocality in time, or temporal dispersion, has been already taken into account for wave propagation in porous media [4,5,6,7]. That means, in Fourier space the effective density and bulk modulus depend on the frequency ω. In other terms, the field dynamics at one location retains a memory of the field values at this location but is not affected by the neighboring values. The local description is usually based on retaining only the leading order terms in the two-scale homogenization method [8, 9, 7] using an asymptotic approach in terms of a characteristic length of the medium, the period L in periodic media, supposed to be much smaller than the wavelength λ [10].

The nonlocal theory we propose takes not only temporal dispersion but also spatial dispersion into account. Here, the medium is assumed unbounded and homogeneous, so that spatial dispersion refers to the dependence of the permittivities -effective density and bulk modulus -on the wavenumber k [11].

The materials susceptible to show the nonlocal behaviour may be classified into two main groups regarding their microgeometry. The first comprises the materials who exhibit this behaviour in sufficiently high frequency regime. The second one concerns materials with microgeometry constituting the resonators, which exhibit the spatial dispersion phenomena even at not very high frequencies; the resonance phenomena act as a source generating nonlocal behaviour. In this chapter we investigate the first type of these geometries, and will see the second one in a forthcoming chapter where the geometry of daisy chained Helmholtz resonators will be treated.

The nonlocal theory we use here takes advantage of an analogy with electromagnetics to give a coarse-grained description of dynamics of small amplitude perturbations in the porous media; expressing the macroscopic governing equations in a Maxwellian form. The homogenization method employed in the present nonlocal theory results in the remarkable point that considerations on length-scale constraints inherent in local theory, originated from asymptotic approach, do not exist any more. For the microgeometry considered in this chapter, the latter is explicitly shown by the fact that the whole dynamics is described through the nonlocal Maxwellian approach; the normal mode related phase velocity is precisely predicted by this approach in a large frequency range. The chapter is organized as follows. In section 3.2 microscopic equations of a small perturbation in viscothermal fluids including balance laws and constitutive relations are expressed. Here, microscopic scale refers to the scale in which Navier-Stokes and Fourier equations are valid. In section 3.3 we will see how the macroscopic fields are defined through a spatial averaging in the present periodic media. In section 3.4 we review briefly the macroscopic governing equations and constitutive relations in local theory. The procedure to determine effective frequency-dependent density and bulk modulus through two actionresponse problems, is presented as well. Once these two effective properties are known, we can get directly the phase velocity of the single mode propagating and attenuating in the medium. In section 3.5 the nonlocal theory [1] is briefly presented for this case of periodic media and propagation according to a symmetry axis, here considered the positive x-axis with the unit vector e x . It will be shown how the wavenumber and frequency dependent effective density and bulk modulus are determined by solving two distinct systems of microscopic equations, coming from two independent action-response problems. Once these two effective functions are known, we can have access after solving a dispersion equation, to the phase velocity of the possible Bloch wavemodes propagating and attenuating in the medium. In section 3.6 we introduce the multiple scattering method including viscothermal effects, allowing to obtain the spectrum of Bloch wavenumbers for the geometry shown in Fig. 3.1. The aforementioned microscopic systems of equations leading to determine the effective properties of the medium in the framework of local and nonlocal theories are solved by Finite Element Method (FEM) using FreeFem++ [12]. The phase velocities coming from these two theories are compared to those obtained by the quasi-exact multiple scattering method in section 3.7, where we will observe clearly the power of the new nonlocal theory.

Microscopic equations

At the microscopic scale, the linear equations governing the dynamics of smallamplitude disturbances in a homogeneous viscothermal fluid come from balance equations of mass, momentum and energy, the constitutive relations of Navier-Stokes and Fourier, and the state equation of the fluid. These governing equations describe the small deviations of thermodynamic pressure p, density ρ, temperature T , velocity v and entropy s, from their rest state p 0 , ρ 0 , T 0 , v 0 = 0 and s 0 , up to the terms of first order. The two constitutive relations are written as

σ ′ ij = 2η e ij - 1 3 (∇.v)δ ij + ζ(∇.v)δ ij (3.1a) q = -κ∇T (3.1b)
The first one, is a linear relation between the shear stress σ ′ ij and strain rate, where e ij = 1 2 (∂ i v j + ∂ j v i ) is the symmetric part of the strain rate tensor, δ is the Kronecker symbol, and η and ζ are the first and second viscosity of the fluid. The second one, is the heat conduction Fourier's law, with q the heat flow, and κ coefficient of thermal conductivity.

Using these constitutive relations, the conservation equations of mass, momentum and energy in the bulk fluid V f for a fluid particle give

∂b ∂t + ∇ • v = 0 (3.2a) ρ 0 ∂v ∂t = -∇p + η∇ 2 v + ζ + 1 3 η ∇ (∇ • v) (3.2b) ρ 0 c p ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ (3.2c)
where b ≡ ρ ′ /ρ 0 with ρ ′ the density deviation, τ is the excess temperature, β 0 ≡ ρ 0 [∂(1/ρ)/∂T ] p and c p ≡ T 0 (∂s/∂T ) p represents the coefficient of thermal expansion and the specific heat at constant pressure, which are evaluated at the fluid rest sate. For convenience in future writing, we denote also by p the pressure deviation. When we expand the thermodynamic equations of state, ρ = ρ(p, s) and T = T (p, s) near the rest state up to the first term [13], then by omitting s in these equations and making use of the thermodynamic identities (∂ρ/∂s) p = -ρ 0 β 0 /c p , (∂T /∂p) s = β 0 T 0 /ρ 0 c p , c 2 0 ≡ (∂p/∂ρ) s representing the adiabatic sound speed squared, we conclude the following state equation in

V f γχ 0 p = b + β 0 τ (3.3)
where χ 0 ≡ ρ -1 0 (∂ρ/∂p) s is the coefficient of adiabatic compressibility at rest state, γ ≡ c p /c v the relative specific heats at constant pressure and constant volume, involved in the thermodynamic identity γ -1 = β 2 0 T 0 /ρ 0 c p .

In the solid phase region V s energy balance equation is reduced to

ρ s c s p ∂τ s ∂t = κ s ∇ 2 τ s (3.4)
where ρ s is the constant solid density, τ s solid excess temperature, and κ s solid coefficient of thermal conductivity.

On the solid-fluid interface ∂V, we have the conditions of the continuity of the temperature τ = τ s and heat flow κ∇τ = κ s ∇τ s . We admit in the following that the specific heat and coefficient of heat conductivity of the solid phase are sufficiently large to allow that the latter conditions combined with Eq.(3.4) are reduced to a single boundary condition τ = 0. Taking into account the no-slip condition on the fluid-solid interface, the boundary conditions for the velocity and excess temperature on ∂V are finally written as

v = 0 (3.5a) τ = 0 (3.5b)
The equations (3.2) and (3.3) with boundary conditions (3.5) establish a closed system with the field variables v, b, p and τ .

Averaging

We need first to define our macroscopic fields in order to describe their dynamics, specializing to the case of periodically structured porous media. We use here the spatial averaging method following Lorentz [14] and refined by Russakoff [15]. Let I be the fluid indicator function

I(r) = 1, r ∈ V f 0, r ∈ V s (3.6)
Given a microscopic field a(t, r) in the fluid region, its macroscopic value is defined by the following integration over the whole space

a (t, r) = dr ′ I(r ′ )a(t, r ′ )f (r -r ′ ) (3.7)
where f is a test function with a typical width L h which is defined as the homogenization length. In periodic media, L h is usually taken to be equal to a period, function f then being constant inside and zero outside. However, this homogenization length can include more than one spatial period which results in a different value for a [1]. Here we choose to take the average just over one irreducible spatial period according to the direction of propagation, which gives the effective properties of the porous medium depending of this choice. The presence of the characteristic function in the above, ensures that the integrand be non-zero only in the fluid region. In addition, the test function is chosen to be normalized over the whole space

drf (r) = 1 (3.8)
The macroscopic homogeneity of the medium implies that the volume fraction of the fluid, i.e., the porosity φ = I is constant all over the medium. The so-called spatial averaging theorem [16] is written as

∇a = ∇ a + ∂V dr ′ a(t, r ′ )n(r ′ )f (r -r ′ ) (3.9)
relating the average of the gradient of a microscopic field a to the gradient of the averaged field, where n(r ′ ) is the outward normal from the fluid, on the solid-fluid interface ∂V.

These definitions are employed next to formulate the macroscopic local and nonlocal theory.

Local theory

As the nonlocal theory is presented in the form of Maxwell acoustic equations, we present the macroscopic equations of local theory in a Maxwellian form as well, in order to compare these two theories and see more clearly the difference between them.

The macroscopic velocity and condensation are defined as

V = v (3.10a) B = b (3.10b)
Since the velocity vanishes on the pore walls, the following direct commutation relation always holds

∇ • v = ∇ • v = ∇ • V (3.11)
Thus, the averaged form of Eq.(3.2a) become

∂B ∂t + ∇ • V = 0 (3.12)
The electromagnetic analogy then suggests that the system of macroscopic equations can be carried through by introducing new Maxwellian fields H and D, and also operators ρ and χ-1 , such that Assuming that H and ρ are scalars, we disregard the propagation of macroscopic shear waves. Here, the propagation of longitudinal waves is considered along the symmetry axis e x , we have D = De x and V = V e x , then the constitutive local relations are written as

D(t, x) = t -∞ dt ′ ρ(t -t ′ )V (t ′ , x) (3.16) H(t, x) = t -∞ dt ′ χ -1 (t -t ′ )B(t ′ , x) (3.17)
We see in the above relations that temporal dispersion is taken into account, i.e., the fields D and H at a given place x and time t depend on the history of the fields V and B at the same place. The time invariance of the problem results in the particular t and t ′ time-difference dependence of the density and bulk modulus kernels.

In the local theory it turns out that the abstract 'Maxwell' macroscopic field H is the mean pressure p . Indeed, we shall always define this field H by writing:

pv (t, x) = H(t, x) v (t, x) (3.18)
As p is the thermodynamic excess pressure, and pv is interpreted as the acoustic part of the energy current density [17], this equality (3.18) may be viewed as a thermodynamic definition. The vector S = HV plays the role of an acoustic macroscopic 'Poynting' vector. Now, in the local theory, because the motion is almost divergence-free at the pore scale, the microscopic pressure gradients are always on the order of the macroscopic pressure gradients, and, because scale separation is assumed, the pressure can be viewed in first approximation as a slowly variable quantity equal to the mean pressure. Thus in (3.18), p can be replaced by p and extracted from the average; this leads to identifying H = p .

The Fourier transform of the constitutive relations (3.16) and (3.17) are written as

D(ω, x) = ρ(ω)V (ω, x) (3.19) H(ω, x) = χ -1 (ω)B(ω, x) (3.20)
We proceed next to review briefly how one can have access to the effective functions ρ(ω) and χ -1 (ω) from microscopic fields.

Determination of constitutive operators

The procedure to obtain effective properties of the medium in local theory logically derives from the only assumption that, because of scale separation, the motion may be viewed as divergence-free in first approximation, at the pore scale. There is however no complete generality in this assumption. It is a sort of simplification of the true wave problem which may be in error in geometries with resonators.

The two-scale homogenization using asymptotic analysis is usually employed to justify the procedure. This technique provides a powerful mathematical formalization of the above tacit physical assumption.

Two characteristic lengths are introduced: the wavelength λ and characteristic length of the unit cell L. There, it is supposed that the wavelength is much bigger than the characteristic unit cell length: ε ≡ L/λ ≪ 1. The microscopic fields are expanded involving the parameter ε, leading further to derive two independent sets of equations by which we can compute effective density and bulk modulus.

As such, the Fourier kernels density and bulk modulus are obtained via two independent action-response problems.

When a harmonic bulk force f (t) = f 0 e -iωt e x , with constant f 0 , is applied on the fluid, we will have the following action-response problem involving the amplitudes of the fields

∇.v = 0 (3.21a) -iωv = ∇p + η∇ 2 v + f 0 (3.21b) in V f v = 0 (3.22)
on ∂V, where the fields are the amplitudes of the solutions

v(t, r) = v(ω, r)e -iωt (3.23a) p(t, r) = p(ω, r)e -iωt (3.23b)
This problem is one of the two independent problems, obtained at leading order, by the aforementioned homogenization method. It is suitable to determine the density ρ(ω) in the local theory because, the above-mentioned tacit physical assumption is encapsulated in (3.21b), and, coherent with this simplification, the neglect of spatial dispersion is apparent in the fact that f is taken as a spatial constant.

Considering a periodic square cell with the length L, containing a single cylinder (Fig. From the solution field v the equivalent-fluid density is obtained as

ρ(ω) = - f 0 iωV (ω) (3.24)
Applying an excitation in the form of a stirring heating Q(t) = Q0 e -iωt , where Q0 is a constant, leads to the following action-response problem [7] for the amplitude of the excess temperature field τ (t, r) = τ (ω, r)e -iωt

-iωρ 0 c p τ = κ∇ 2 τ + Q0 (3.25) in V f τ = 0 (3.26) on ∂V.
This problem is the second one obtained at leading order, by the mentioned homogenization. It turns out suitable to determine the compressibility χ(ω) in the local theory because, the tacit physical assumption that the pressure field is a slowly variable quantity that may be viewed in first approximation as equal to the mean pressure, is encapsulated in (3.25) in the very fact that Q0 is taken as a spatial constant. This is directly the consequence of the neglect of spatial dispersion.

Considering as before the periodic square cell with the length L, containing a single cylinder (Fig. 3.1) and bounded in x ∈ [0, L], y ∈ [-L/2, L/2], there is a unique amplitude field τ (ω, r) solution to Eqs. (3.25-3.26), which is periodic with the period L such that it gives the same value on the cell boundaries x = 0 and x = L, ∀y; and also on the cell boundaries y = -L/2 and y = L/2, ∀x.

From the solution field τ , a factor ρ ′ analogous to the previous ρ is obtained as

ρ ′ (ω) = - Q0 iωT (ω) (3.27)
where T = τ .

In the framework of nonlocal theory, the following direct relation exists between the two functions ρ ′ and χ:

χ(ω) = χ 0 γ -(γ -1) ρ 0 c p ρ ′ (ω) (3.28)

Nonlocal theory 77

The reasoning to obtain this relation has been given in Appendix A; indeed, the fact that we were considering cylindrical circular tubes was not explicitly used.

The bulk modulus is thus written as

χ -1 (ω) = χ -1 0 γ + (γ -1) iωρ 0 c p τ (ω, r) Q0 -1
(3.29)

Phase velocity

Once the density ρ(ω) and bulk modulus χ -1 (ω) (or compressibility χ(ω)) are determined, we can obtain the constant of the medium for each frequency. For a given frequency ω, there is only one single normal mode that can propagate in the given positive x direction. With this single mode is associated a wavenumber q(ω) verifying the relation

ρ(ω)χ(ω)ω 2 = q 2 (3.30)
such that ℑ(q) > 0. The complex phase velocity c(ω) associated with this frequency ω is immediately written as

c(ω) = ρ -1 (ω)χ -1 (ω) (3.31)

Nonlocal theory

We intend here to write the macroscopic equations in a Maxwellian form allowing for both temporal and spatial dispersion. The field equations in nonlocal theory will be the same as in local theory. As before, the macroscopic condensation and velocity are defined as

V = v (3.32a) B = b (3.32b) The relation ∇ • v = ∇ • v = ∇ • V is
as well valid because of the boundary condition on the velocity on the pore walls. Thus, the averaged form of Eq.(3.2a) become

∂B ∂t + ∇ • V = 0 (3.33)
Here also, the electromagnetic analogy then suggests that the system of macroscopic equations can be carried through by introducing new Maxwellian fields H and D, and also operators ρ and χ-1 , such that The operators ρ and χ-1 are density and bulk modulus, describing the effective properties of the medium. They are uniquely fixed in principle, through the condition that H is to be identified through the acoustic part of energy current density S = pv , by setting [1] pv = H v (3.37)

The propagation of longitudinal waves is considered along the symmetry axis e x , we have D = De x and V = V e x , then the nonlocal constitutive relations, this time, are written as

D(t, x) = t -∞ dt ′ dx ′ ρ(t -t ′ , x -x ′ )V (t ′ , x ′ ) (3.38) H(t, x) = t -∞ dt ′ dx ′ χ -1 (t -t ′ , x -x ′ )B(t ′ , x ′ ) (3.39)
We see in the above relations that not only temporal dispersion but also spatial dispersion is taken into account, i.e., the fields D and H at a given time t depend on the fields V and H at all previous time and all points of the space. The time invariance and macroscopic homogeneity of the problem result in the dependence of the kernels on the differences t -t ′ and x -x ′ . The Fourier transform of the constitutive relations (3.38) and (3.39) are written as

D(ω, k) = ρ(ω, k)V (ω, k) (3.40) H(ω, k) = χ -1 (ω, k)B(ω, k) (3.41)
We proceed next to review briefly how one can have access to the effective functions ρ(ω, k) and χ -1 (ω, k) from microscopic fields.

Determination of constitutive operators

The above Fourier coefficients ρ(ω, k) and χ -1 (ω, k) are directly related to the macroscopic response of the permeating fluid subjected to a harmonic fictitious pressure term P(t, x) = P 0 e -iωt+ikx added to the pressure, either in the Navier-Stokes Eq.(3.2b), or the Fourier Eq.(3.2c).

Thus to determine the kernel ρ(ω, k) we first consider solving the action-response problem

∂b ∂t + ∇ • v = 0 (3.42a) ρ 0 ∂v ∂t = -∇p + η∇ 2 v + ζ + 1 3 η ∇ (∇ • v) + f (3.42b) ρ 0 c p ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ (3.42c) γχ 0 p = b + β 0 τ (3.42d) in V f , and v = 0 (3.43a) τ = 0 (3.43b) on ∂V.
The stirring force appears in the form of

f = -∇P = -ike x P 0 e -iωt+ikx (3.44)
The unique solutions to the above system (3.42a-3.44), for the fields v, b, p, τ , take the form

v(t, r) = v(ω, k, r)e -iωt+ikx (3.45a) b(t, r) = b(ω, k, r)e -iωt+ikx (3.45b) p(t, r) = p(ω, k, r)e -iωt+ikx (3.45c) τ (t, r) = τ (ω, k, r)e -iωt+ikx (3.45d)
The response amplitudes v(ω, k, r), b(ω, k, r), p(ω, k, r), and τ (ω, k, r) are bounded functions which are uniquely determined by the microgeometry. For the special case of periodic geometry, the solution is not unique, however. The unique formulation of the action response problem also includes a specification of the periodic cell. The condition that the response amplitudes are bounded functions is replaced by the condition that they are periodic over the chosen cells. This indeterminacy was absent in the case of the local theory, because the solutions were independent of the choice of the periodic cell. Here, in all foregoing calculations, the irreducible cell as illustrated in Figs 1 and 2 will always be chosen as the period.

The above problem, once solved, we use the fundamental 'Umov-Poynting' relation (3.37) to write

P v = pv (3.46)
where P = P(ω, k)e -iωt+ikx is the macroscopic part of the pressure response p(t, r) = p(ω, k, r)e -iωt+ikx , whose amplitude is determined by

P(ω, k) = p(ω, k, r)v(ω, k, r) .e x V (ω, k) (3.47)
Then using the Fourier transform of Eq.(3.34), applying Eq.(3.40), and admitting that the two parts P and P 0 just add to form the field H, viz.

-iωρ(ω, k)V (ω, k) = -ik(P(ω, k) + P 0 ) (3.48)
gives rise immediately to nonlocal Equivalent-fluid density ρ(ω, k)

ρ(ω, k) = k(P(ω, k) + P 0 ) ωV (ω, k) (3.49)
The strong motivation for this conjectured expression (3.49), is its simplicity and the fact that it is explicitly verified in absence of solid [1]. It has been exactly verified in cylindical circular tubes (chapter 2). Thus it must provide the exact upscaling procedure.

At this point, we see that the fields p(ω, k, r) and v(ω, k, r) are needed to be known in order to determine from microgeometry the effective density of the fluidsaturated porous medium. Hence, instead of solving (3.42a-3.43b) it is sufficient to solve the following system of equations to get the amplitudes of the fields in

(3.45) -iωb + ∇.v + ikv x = 0 (3.50a) -iωρ 0 v = -∇p -ikpe x + η∇ 2 v + 2ikη ∂v ∂x (3.50b) -ηk 2 v + ζ + 1 3 η ∇(∇.v) + ik ζ + 1 3 η (∇.v)e x +ik ζ + 1 3 η ∇v x -ζ + 1 3 η k 2 v x e x -ike x P 0 -iωρ 0 c p τ = -iωβ 0 T 0 p + κ∇ 2 τ + 2ikκ ∂τ ∂x -k 2 κτ (3.50c) γχ 0 p = b + β 0 τ (3.50d) in V f , and v = 0 (3.51a) τ = 0 (3.51b) on ∂V.
As in local case, considering a periodic square cell with the length L, containing a single cylinder (Fig. 3.1) and bounded in x ∈ [0, L], y ∈ [-L/2, L/2], the amplitude fields v(ω, k, r), b(ω, k, r), p(ω, k, r), and τ (ω, k, r) are periodic with the period L, such that they give the same values on the cell boundaries x = 0 and x = L, ∀y; and also on the cell boundaries y = -L/2 and y = L/2, ∀x.

The procedure to determine the kernel χ -1 (ω, k) is quite similar but slightly less direct, as was already the case in local theory. We now deal with the field B, and the way the latter is thought to be connected with the fields appearing in the new action-response problem, requires a little reflection.

We consider again, initially solving the action-response problem with an harmonic fictitious term P(t, x) = P 0 e -iωt+ikx added to the pressure, but appearing this time in the energy balance equation

∂b ′ ∂t + ∇ • v = 0 (3.52a) ρ 0 ∂v ∂t = -∇p + η∇ 2 v + ζ + 1 3 η ∇ (∇ • v) (3.52b) ρ 0 c P ∂τ ∂t = β 0 T 0 ∂p ∂t + κ∇ 2 τ + Q (3.52c) γχ 0 p = b ′ + β 0 τ (3.52d) in V f , and v = 0 (3.53a) τ = 0 (3.53b)
on ∂V, with the stirring heating

Q = β 0 T 0 ∂P ∂t = -iωβ 0 T 0 P 0 e -iωt+ikx (3.54)
The solutions to the above problem take the same form as specified before through Eqs.(3.45) and the same comments can be made regarding the special case of periodic microstructure. This problem, once solved, we again use the fundamental 'Umov-Poynting' relation-definition (3.37) to write

P v = pv (3.55) and thus P(ω, k) = p(ω, k, r)v(ω, k, r) .e x V (ω, k) (3.56)
for the amplitude of the macroscopic part P = P(ω, k)e -iωt+ikx of the pressure response p(t, r) = p(ω, k, r)e -iωt+ikx .

Then using the Fourier transform of Eq.(3.39), and admitting as before that the two parts P and P 0 just add to form the field H, viz.

P(ω, k) + P 0 = χ -1 (ω, k)B(ω, k) (3.57)
it remains to identify the field B.

This identification is obtained from the thermodynamic understanding that the field b ′ in the action-response problem (3.52-3.53), as it is fixed in particular by the Laplacian term expressing thermal conduction in Eq.(3.52c), is to be viewed as determining a nonisothermal response part in the macroscopic condensation field B. The latter is thus seen as the direct sum of two contributions: B = b = b ′ + b ′′ , with b ′ a nonisothermal response part determined by the above action-response problem, and b ′′ a complementary isothermal response part, by definition given by the isothermal relation b ′′ = γχ 0 P 0 . On account of the fact that P 0 = φP 0 , this results in the following relation [1] P(ω, k)

+ P 0 = χ -1 (ω, k) b ′ (ω, k, r) + φγχ 0 P 0 (3.58)
That gives the nonlocal equivalent-fluid bulk modulus χ -1 (ω, k)

χ -1 (ω, k) = P(ω, k) + P 0 b ′ (ω, k, r) + φγχ 0 P 0 (3.59)
Again, the strong motivation for the above simple ansatz (3.59) is its simplicity, and the fact that it is explicitly verified in absence of solid [1] and in cylindrical circular tubes (chapter 2). Thus we expect that it provides the exact upscaling procedure.

At this point, to determine the effective bulk modulus of the fluid-saturated porous medium, we have to look for the amplitude fields p(ω, k, r) and v(ω, k, r), and b ′ (ω, k, r). These fields can be obtained by substituting (3.45) in (3.52a-3.53b) which gives the following system

-iωb ′ + ∇.v + ikv x = 0 (3.60a) -iωρ 0 v = -∇p -ikpe x + η∇ 2 v + 2ikη ∂v ∂x (3.60b)
-ηk 2 v + ζ + 1 3 η ∇(∇.v) + ik ζ + 1 3 η (∇.v)e x +ik ζ + 1 3 η ∇v x -ζ + 1 3 η k 2 v x e x -iωρ 0 c p τ = -iωβ 0 T 0 p + κ∇ 2 τ + 2ikκ ∂τ ∂x -k 2 κτ -iωβ 0 T 0 P 0 (3.60c) γχ 0 p = b ′ + β 0 τ (3.60d) in V f , and v = 0 (3.61a) τ = 0 (3.61b) on ∂V.
Similarly, here also, considering a periodic square cell with the length L, containing a single cylinder (Fig. 3.1) and bounded in x ∈ [0, L], y ∈ [-L/2, L/2], the amplitude fields v(ω, k, r), b(ω, k, r), p(ω, k, r), and τ (ω, k, r) are periodic with the period L, such that they give the same values on the cell boundaries x = 0 and x = L, ∀y; and also on the cell boundaries y = -L/2 and y = L/2, ∀x.

Phase velocities

Contrary to the case of local theory, here, since we fully take into account spatial dispersion, several normal mode solutions might exist, with fields varying as e -iωt+iqx . Each solution should satisfy the following dispersion equation ρ(ω, q)χ(ω, q)ω 2 = q 2 (3.62)

With each frequency ω, several wavenumbers q l (ω), ℑ(q l ) > 0, l = 1, 2, ..., may be associated. The complex phase velocity corresponding to a given solution is written as

c l (ω) = ω q l (ω) (3.63)

Multiple scattering method

For the simple geometry represented in Fig. (3.1) a relatively simple calculation of the possible wavenumbers k l is feasible by a multiple scattering approach [3]. In this calculation, which is presented in some more detail in what follows, we adopt a description of the fluid motion in terms of three velocity potentials, the acoustic potential φ a , entropic potential φ e and vorticity potential ψ:

v = ∇ (φ a + φ e ) + ∇ × ψ (3.64)
The vorticity potential ψ has just one component, which is directed along the z-axis and is denoted by φ v . In harmonic regime, three independent Helmholtz equations

∇ 2 + (k α ) 2 φ α = 0, α = a, e, v (3.65)
are satisfied in V f , where (k α ) 2 , α = a, e, v, are the wavenumbers associated to acoustic, thermal and viscous waves, respectively. The former two (k a ) 2 and (k e ) 2 are the opposite of the small and large solutions λ 1 and λ 2 of Kirchhoff-Langevin's dispersion equation (see Eq.(2.34) in chapter 2); the latter is (k v ) 2 ≡ iω/ν.

Using Eqs.(2.40a), (2.40d), and (2.28) in chapter 2, it is easy to express the excess temperature in terms of potentials

β 0 τ γ -1 = ( κ ρ 0 c v + iω λ a ) -1 φ a + ( κ ρ 0 c v + iω λ e ) -1 φ e (3.66)
The boundary conditions at the solid-fluid interface for the potentials come from the fact that the displacement and excess temperature fields on ∂V are such that

u = 0 (3.67a) τ = 0 (3.67b)
These boundary conditions establish the coupling between the potentials, which results in the fact that a wave carried by one potential is scattered on the three types of waves on interacting with the solid.

In this chapter, because we investigate only the least attenuated mode and are not concerned with the terms of minor importance that concern the intrinsic bulk fluid propagation, the bulk attenuation will be neglected for the acoustic mode. This simplification will allow direct use of known results on Schlömilch series.

Thus we set as a simplification -λ a = (k a ) 2 ≡ (ω/c 0 ) 2 . Moreover, we also neglect the higher order terms governing the attenuation of entropic waves, and set as another simplification -λ e = (k e ) 2 ≡ iωρ 0 c p /κ. Consistent with the first simplification, we need not account for the thermal conductivity term in the first parenthesis in (3.66). After straightforward calculation using the thermodynamic identity γ -1 = T 0 β 2 0 c 2 0 /c p , the following relation is obtained

τ = T 0 β 0 c p iωφ a + ρ 0 c p β 0 κ φ e (3.68)
Considering one row of infinite number of cylinders, as is shown in Fig. 3.2, we expand the potentials in terms of right and left going plane waves

φ α 0 (r) = ∞ n=-∞ A +α 0n e ik α n .r + A -α 0n e -ik α n .r (3.69a)
φ α L (r) = ∞ n=-∞ A +α Ln e ik α n .(r-Lex) + A -α Ln e -ik α n .(r-Lex) (3.69b)
The ingoing or outgoing meaning of the four types of amplitudes A is apparent on the figure. The index α refers to the type a, e, or v of potential field. It is clear that the periodicity of the potential fields with respect to y coordinates implies that for each n the y component of the wavevectors k α n should be k α ny = 2πn/L, thus (k α ) 2 = (k α nx ) 2 + (2πn/L) 2 . Another symmetry consideration of the problem is that we are interested only with the solutions leading to a fluid motion symmetric around each cylinder. In this, we restrict to the motions that can be created by a 'macroscopic stirring'. This is analogous to our restriction in chapter 2, for the same reason, to axisymmetric motions. This restriction implies that the fields φ a and φ e are even functions, and φ v an odd function with respect to y coordinates. Thus, with regard to the terms in (3.69), after combining the up and down components n and -n there will appear cos(2πny/L) y-dependency for acoustic and entropic potentials, and sin(2πny/L) y-dependency for vorticity potential. To account explicitly for this symmetry in the notation, instead of (3.69) we employ the following condensed form of the potentials

φ α 0 (r) = ∞ n=0 C α n (y) A +α 0n e ik α nx x + A -α 0n e -ik α nx x (3.70a) φ α L (r) = ∞ n=0 C α n (y) A +α Ln e ik α nx (x-L) + A -α Ln e -ik α nx (x-L) (3.70b)
where

C α n (y) =    cos 2πny L , α = a, e sin 2πny L , α = v (3.71)
We note also, that to each n, α, and ω, we may associate a characteristic incidence angle θ α n , such that

k α sin(θ α n ) = 2πn L , k α cos(θ α n ) = k α nx (3.72)
For the acoustic type α = a, this angle will be real when the frequency is such that 2πn/(k α L) < 1, and complex, equal to π/2 -iξ at higher frequencies, with ξ > 0 ensuring that ℑ(k α nx ) > 0. For the entropic and vorticity types, this angle will be complex, chosen such that ℑ(k α nx ) > 0.

The first step in the calculation is to obtain the reflection and transmission properties of the row of cylinders, by the following scattering matrix, which relates the outgoing waves to the ingoing ones

  A - 0 A + L   =   T R R T     A - L A + 0   (3.73)
where

A - 0 =        A -a 0 A -e 0 A -v 0        , A + L =        A +a L A +e L A +v L        (3.74) O x y 0 L L 2 -L 2 A - 0 A + 0 A + L A - L Figure 3.2:
One row of infinite number of rigid cylinders.

and so on for the vectors A + 0 and A - L . Each of the vectors A +α L , A +α 0 , A -α L , and A -α 0 contains the whole ensemble of plane wave amplitudes with α = a, e, v, each of which had been indexed by n. The reflection and transmission matrices R and T thus have elements of the type R αβ pn and T αβ pn , where the indexes on the right refer to ingoing waves and the index on the left to outgoing ones. The presence of the different elements results from the interactions and transformations of the different kinds of potentials into one another, due to the boundary conditions (3.67). To compute R and T and thus construct the scattering matrix, the analysis of the scattering problem is divided in different elementary parts, combined in the end.

Consider a given ingoing potential φ β n,in (r) = C β n (y)e ik β nx x coming from the left on the row, with n arbitrary and β which is set to be replaced by either a, e or v. This ingoing potential of type β, n, by the scattering effects and through the boundary conditions, creates outgoing potentials of all of the three types a, e and v, and indexes p. These outgoing potentials can be interpreted either as reflected or transmitted fields by the row

φ αβ n,R (r) = ∞ p=0 R αβ pn C α p (y)e -ik α px x , α = a, e, v, x ≃ 0 (3.75a) φ αβ n,T (r) = ∞ p=0 T αβ pn C α p (y)e ik α px (x-L) , α = a, e, v, x ≃ L (3.75b)
where φ αβ n,R and φ αβ n,T refer to the reflected and transmitted potentials of type α, which have been created by the ingoing β, n-field. The coefficients T αβ pn and R αβ pn are the above-mentioned elements of the matrices R and T .

On the other hand, the same reflection and transmission fields can be expressed in terms of the potential fields scattered by all cylinders belonging to the row. The scattered field from the j-th cylinder at the position r j can be expanded on the basis of Hankel functions of the first kind

φ αβ n,scat (r, r j ) = ∞ m=-∞ B αβ mn i m H m (k α |r -r j |)e imθ r-r j (3.76)
where θ r-r j is the azimuthal angle of the vector rr j relative to e x , and B αβ mn is the unknown weighting coefficient, associated with the scattered or Hankel divergent wave H m . The row of cylinders being infinite, this coefficient is independent of the cylinder under consideration. Thus, the scattered field by the row is written by summing the above expression over all cylinders

φ αβ n,scat (r) = ∞ j=-∞ ∞ m=-∞ B αβ mn i m H m (k α |r -r j |)e imθ r-r j . (3.77)
We can now identify the reflected and transmitted fields through the following relations

φ αβ n,R (r) = φ αβ n,scat (r), x ≃ 0 (3.78a) φ αβ n,T (r) = φ αβ n,scat (r) + δ αβ φ β n,in (r), x ≃ L (3.78b)
where δ represents the Kronecker symbol. It is clear that with the help of the addition theorem of Bessel functions these relations determine in principle the reflection and transmission coefficients R αβ pn and T αβ pn in terms of the coefficients B αβ mn . To determine the coefficients B αβ mn , we proceed as follows.

With the same ingoing potential φ β n,in (r) = C β n (y)e ik β nx x coming from the left on the row, we can analyze the field incident on an arbitrary cylinder situated at r l , in terms of the convergent Bessel waves. The three types α = a, e, v of incident potentials will be present because of the row scattering, and each can be expanded in the basis of Bessel functions centred at r l as

φ αβ n,inc (r, r l ) = ∞ m=-∞ C αβ mn i m J m (k α |r -r l |)e imθ r-r l (3.79)
Again, the unknown weighting coefficient C αβ mn , associated with the incident or convergent Bessel wave J m , is independent of the cylinder under consideration.

Alternatively, this incident field on cylinder r l can be regarded as the sum of the ingoing field, and the fields scattered by cylinders j = l. According to this viewpoint, the incident field on cylinder r l is given by

φ αβ n,inc (r, r l ) = j =l φ αβ n,scat (r, r j ) + δ αβ φ β n,in (r, r l ) (3.80)
For each ingoing β, n-potential, we now apply the boundary conditions (3.67) to the total field resulting from the superposition of incident (3.79) and scattered (3.76) fields around one arbitrary cylinder. This allow to relate the coefficients B αβ mn to the C αβ mn as follows

B αβ mn = γ=a,e,v D γβ m C γβ mn , α = a, e, v (3.81)
through the coefficient D γβ m which is defined by

       k a H ′ m (k a R) k e H ′ m (k e R) m R H m (k v R) m R H m (k a R) m R H m (k e R) k v H ′ m (k v R) T 0 β 0 cp iωH m (k a R) ρ 0 cp β 0 κ H m (k e R) 0               D aβ m D eβ m D vβ m        = K β (3.82)
for β = a, e, v, where the vectors K β are defined as where,

K a ≡        -k a J ′ m (k a R) -m R J m (k a R) -T 0 β 0 cp iωJ m (k a R)        , K e ≡        -k e J ′ m (k e R) -m R J m (k e R) - ρ 0 cp β 0 κ J m (k e R)        , and K v ≡        -m R J m (k v R) -k v J ′ m (k v R) 0       
G α m (r, r l ) =        j =l i m H m (k α |r -r j |) cos(mθ r-r j ), α = a, e j =l i m H m (k α |r -r j |) sin(mθ r-r j ), α = v (3.84)
which can be rewritten using the addition theorem of Bessel functions as

G α m (r, r l ) =                      ∞ q=0 σ |m-q| (k α L) + σ m+q (k α L) i q J q (k α |r -r l |) cos(qθ r-r l )(1 - 1 2 δ q0 ), α = a, e ∞ q=0 σ |m-q| (k α L) -σ m+q (k α L) i q J q (k α |r -r l |) sin(qθ r-r l ), α = v (3.85)
where σ m (k α L) is a function representing the series of Hankel functions of the first kind defined by the relation

σ m (k α L) = σ -m (k α L) = (1 + (-1) m ) ∞ j=1 H m (jk α L), α = a, e, v (3.86)
which vanishes for the odd values of m. Notice that for the reason of better convergence, it is necessary to expand theses series in terms of the Schlömilch series [18]. In this expansion it is assumed that the wavenumber k a is real, which is the case as we neglect as explained before the bulk attenuation in the fluid.

Finally, the relation (3.80) is expressed in the basis of Bessel functions with respect to the coordinates centered at r l . First, the known formula e ik.(r-r l ) = ∞ m=-∞ i m J m (k|rr l |)e im(θ k -θ r-r l ) for the expansion of the plane waves on the basis of Bessels functions is used to express the ingoing part

φ β n,in (r, r l ) =                      ∞ m=0 (2 -δ m0 )e ik β nx L/2 cos(mθ β n ) i m J m (k β |r -r l |) cos(mθ r-r l ), β = a, e ∞ m=0 -2i(1 -δ m0 )e ik β nx L/2 sin(mθ β n ) i m J m (k β |r -r l |) sin(mθ r-r l ), β = v (3.87)
Next, the addition theorem of Bessel functions is used to express the scattered part.

We obtain in this manner three series of algebraic equations each of which corresponds to either an ingoing acoustic, entropic or vorticity field, to determine the coefficients C αβ mn and, as a result B αβ mn :

C αβ m =                      δ αβ (2 -δ m0 )e ik β nx L/2 cos(mθ β n ) + ∞ q=0 γ=a,e,v D αγ q C γβ q F α mq (σ), α = a, e, v; β = a, v -2iδ αβ e ik β nx L/2 sin(mθ β n ) + ∞ q=0 γ=a,e,v
D αγ q C γβ q F α mq (σ), α = a, e, v;

β = v (3.88)
where, Now, we consider an infinite number of rows separated by the distance L. We make use of the concept of scattering matrix which we have just studied for an arbitrary row, and apply the Bloch condition for this case of periodic medium. We have

F α mq (σ) = (1 - 1 2 δ m0 )(δ αa + δ αe ) [σ m+q (k α L) + σ m-q (k α L)] (3.89) +δ αv [σ m+q (k α L) + σ m-q (k α L)] , α = a,
A + L A - L = e ik B L A + 0 A - 0 (3.92)
where k B denotes the Bloch wavenumber to be determined where 0 and I are the are the zero and identity matrices, respectively. Since at this stage the reflection and transmission matrices R and T are known, we are able to solve the above eigenvalue problem numerically and get the macroscopic Bloch wavenumbers of the medium. As in the nonlocal theory, here also, with each frequency ω there might be associated several Bloch wavenumbers k B,l , l = 1, 2, 3, ..., and their corresponding complex phase velocities will be

c l (ω) = ω k B,l (3.94)

Results

In this section, we will present the results concerning the phase velocity of the single mode obtained through the local theory, the phase velocity of the least attenuated mode obtained using the nonlocal theory, and the phase velocity of the least attenuated Bloch mode coming from the multiple scattering method. These results are shown in a large frequency range for the three different values

Results

93

of porosity φ = 0.99, 0.9 and 0.7. For a constant period length L = 10µm, the corresponding cylinder radius R = L (1 -φ)/π for the three values of porosity becomes R = 0.56µm, 1.78µm, and 3.1µm, respectively. Fluid properties for all computations are presented in Table 3.1. We note that we have taken ζ = 0 for the second coefficient of viscosity, involved in nonlocal theory computations, in order to be closer to multiple scattering calculation in which the dissipation due to compression/dilation motions is not taken into account. As we have seen in section 3.4 the phase velocity in local theory is given by Eq.(3.31) which requires the Fourier effective density ρ(ω) and bulk modulus χ -1 (ω) of the medium. These two functions are computed through two independent sets of equations (3.21a-3.22) and (3.25-3.26).

ρ 0 T 0 c 0 η ζ κ χ 0 c p γ (kg/m 3 ) (K) m/s (kg ms -1 ) (kg ms -1 ) (W m -1 K -1 ) (P a -1 ) (J kg -1 K -1 ) 1.2 293 343 1.8 × 10 -5 0 2.6 × 10 -2 7.1 × 10 -6 1005 1.4
The equations, for local and nonlocal theory, are solved in a periodic square cell including a single cylinder (Fig. 3.2), by Finite Element Method using FreeFem++ as a free software to solve partial differential equations numerically, based on Finite Element Method which can be used for coupled systems. This software has the possibility to generate mesh automatically and is capable of a mesh adaptation, handling the general boundary conditions, to include, now, periodic boundary conditions which is required to solve the present sets of equations. However, the presence of a bug, relating to the handling of boundary conditions, without receiving any error message, created a very long delay to obtain the correct results. This bug has been corrected, finally, by intervention of the developer of the software. FreeFem++ provides us with a powerful tool when the solution of the problem varies locally and sharply, creating a new mesh adapted to the Hessian of the solution. The weak form of the equations to be solved is firstly needed in order to implement the FEM simulations through the software.

To obtain the phase velocity of the least attenuated wave according to nonlocal theory, first we have to obtain the Fourier kernels ρ(ω, k) and χ -1 (ω, k), via solving separately two sets of equations (3.50a-3.51b) and (3.60a-3.61b) by FEM using FreeFem++. For a given k which is involved in excitations (3.44) and (3.54) we solve the two systems of equations for a large frequency range. A priori, we can continue this procedure for several k to approximate finally the complex functions ρ(ω, k) and χ -1 (ω, k) and then using the dispersion relation (3.62) to find the natural constants of the medium q and subsequently the corresponding phase velocities by (3.63). However, these functions do not seem to have a simple form, especially at sufficiently high frequencies where the spatial dispersion effects

Results

95

become important. Thus, to find the complex normal modes of the medium, we have proceeded through Newton-Raphson method. Initial values q 0 are chosen for each given frequency to compute the complex functions ρ(ω, q 0 ) and χ -1 (ω, q 0 ) for a large range of frequency ω, by setting the two excitations and solving the two corresponding systems of equations through FEM simulations. Then, the complex function F (ω, q) = ρ(ω, q)χ(ω, q)ω 2 -q 2 and its derivative ∂F (ω, q)/∂q are calculated at these initial values q = q 0 , to start a Newton iteration q n+1 = q n -F (ω, q n )/[∂F (ω, q)/∂q] q=qn which converges quickly to yield the value of q at which F vanishes. As such, the phase velocity can be immediately determined. In this study, for the chosen frequencies we take for the initial values of q, the fundamental mode obtained by multiple scattering with a discrepancy of 20%. The fundamental modes obtained in such a way by nonlocal theory are found to be in remarkable agreement with those corresponding to multiple scattering method.

In Figures (3. For more concentrated media the discrepancies between local theory and multiple scattering predictions are larger and commence at lower frequencies. As a matter of fact, when the medium becomes more concentrated, the band gaps include larger frequency intervals and the resonance phenomena becomes more influential. These are considered as the signatures of spatial dispersion effects which can be precisely described by quasi-exact multiple scattering calculation.

Contrary to local theory, results issued from nonlocal approach show excellent agreement with those from multiple scattering, regarding the above Figures. These agreements appear to be insensitive to the frequency in which the phase velocity is computed. That was expected by the fact that in nonlocal approach, no length constraint, such as qL ≪ λ, has been considered. However, we note that as the frequency is increased and the medium becomes more concentrated, because of the significant effects due to of spatial dispersion in these cases, the behaviour of wave propagation become more complicated to describe and consequently more precision is required concerning the FEM computations to have an appropriate convergence stability, in the framework of nonlocal theory.

Conclusion

The nonlocal theory in rigid-framed porous media, recently proposed [1], has been here put into evidence through comparing the phase velocity of the least attenuated mode predicted by this theory, and a quasi-exact multiple scattering method, in the case where the microgeometry of the porous medium is in the form of a two-dimensional array of rigid cylinders and the propagation is along one of the perpendicular axis with which the square lattice can be constructed. The phase velocity according to local theory is computed as well in order to observe its domain of validity and its limits in terms of frequency, comparing to aforementioned approaches capable to take into account the spatial dispersion effects.

The quasi-exact multiple scattering calculation of the least attenuated wavenumber including viscothermal effects, which has been already validated [3], is described here in more details. Concerning the local and nonlocal approaches, the two different ways to compute the corresponding wavenumber, by applying the corresponding local and nonlocal procedures to determine effective density and bulk modulus, have been reviewed. We have seen that these procedures lead to solving four independent action-response microscopic problems, each of which associated with a local or nonlocal effective property of the porous medium. These microscopic cell-problems have been solved here through direct numerical simulations by Finite Element Method using FreeFem++, to give then the frequency dependent phase velocities according to local and nonlocal theories.

The results of computations show that with the geometry considered in this chapter, contrary to local theory, nonlocal theory successfully describes the whole dynamics, including high frequency one, where the band gaps are present. In the future we will try to show that by nonlocal theory we can obtain the other axisymmetric modes of the medium, as expected by the nonlocal theory through its resulting dispersion relation.

medium. Given these functions, a Newton-scheme solution of the dispersion equation ρ(ω, k)χ(ω, k)ω 2 = k 2 , allows finding the medium wavenumbers k(ω).

We focus here on the least attenuated mode. The Zwikker and Kosten approximation consists in assuming that in the various elementary duct portions, the propagation of the sound field may be described in terms of Zwikker and Kosten effective density and compressibility [3,4].

Within the same approximations a direct simplified modelling of the propagation of normal modes in the medium will be made in section 4. This modelling will be called here the Bloch-wave modelling. It predicts a typical resonance behaviour of the least attenuated solution. The classical local theory cannot describe this behaviour. On the contrary, the calculations based on nonlocal theory are found to accurately predict the resonance behaviour. This, once again, will provide an unambiguous validation of the proposed general theory.

The two different problems, the action-response one related to nonlocal theory, and the other, which is an eigenvalue problem related to Bloch-wave modelling, are solved in a single periodic square cell of the medium of length L. In Fig. 4.2 this cell is illustrated with the lengths corresponding to its different parts. The widths of the main tube, neck and cavity of the resonator are denoted by Σ, σ, and L -Σ -2l, respectively. The lengths of the tube, neck, and cavity are L, l, and L -l.

Determination of nonlocal effective density

Considering the periodic cell of Fig. 4.2 and the corresponding cell average operation , we look for the response ofthe fluid when a harmonic stirring force f (t, x) = f 0 e -iωt+ikx in the direction of e x is applied in the medium. The direction y goes downward in the neck. For later convenience, we note that our coordinate system (x, y) will be taken such that the point (0, 0) is located at the center of the neck portion. The indicated positions ( 5) and ( 6) in Fig. 4.3 will have the coordinates (0, -l/2) and (0, l/2) respectively.

If we can determine the microscopic response fields velocity v and pressure p then we will have the function ρ(ω, k) through the relation (see chapter 1, Eq.(1.53))

ρ(ω, k) = f 0 -ikP -iω v (4.1) with P(ω, k) = pv v (4.2)
where the v is the x-component of the microscopic velocity v.

In what follows, we make this calculation in analytical simplified manner. We proceed to determine the functions ρ(ω, k) and χ -1 (ω, k) sufficiently precise to give an appropriate modelling of the least attenuated mode. To this aim, we need not consider in full detail the microscopic fields v and p. In the waveguide and cavity, instead of the microscopic fields, we can use with the mean values V x = v • e x and P = p, where the overline denotes the average (at a given x) over the waveguide or the cavity width; and in the neck, we can use with the mean values V y = v • e y and P = p, where the overline denotes the average (at a given y) over the neck width. At the same time, we make some simplifications consistent with describing the propagation of these averaged quantities in terms of the Zwikker and Kosten densities ρ(ω) and bulk modulii χ -1 (ω), in the different slit portions. These depend only on the slit half-widths, which we shall denote by s t , s n , and s c , in the tube, neck, and cavity. The different slit-like tube portions are illustrated in Fig. 4.3. The main tube t is divided in two Zwikker and Kosten ducts, a left duct, and a right duct, oriented in the x direction. The same separation is made for the cavity c, whereas the neck n is not divided but seen as one Zwikker and Kosten duct oriented in y direction.

In Appendix A, the Zwikker and Kosten local theory is expressed for tubes of circular cross-section. For 2D slits, exactly the same general principles of modelling may be used; only some details of the calculations are changed. The cross-section averages become the aforementioned overlines, and the bessel functions J 0 and J 1 are replaced by cosh and sinh functions. Zwikker and Kosten's effective densities ρ α (ω) and bulk moduluii χ -1 α (ω) in the guide, neck and cavity, will be [4]

ρ α (ω) = ρ 0 1 - tanh -iωρ 0 s 2 α /η -iωρ 0 s 2 α /η -1 , α = t, n, c (4.3a) χ -1 α (ω) = γP 0 1 + (γ -1) tanh -iωρ 0 c p s 2 α /κ -iωρ 0 c p s 2 α /κ -1 , α = t, n, c (4.3b)
where the indexes t, n, and c are related to the tube, neck, and cavity respectively, c p is the heat capacity at constant pressure, γ = c p /c v the ratio of the heat capacities at constant pressure and constant volume, κ the coefficient of thermal conductivity, ρ 0 the fluid density at rest, and P 0 the fluid pressure at rest. The corresponding wavenumbers k α (ω) and characteristic admittances Y α (ω) are expressed as

k α = ω c α , α = t, n, c (4.4) Y α (ω) = 2s α ρ α c α , α = t, n, c (4.5)
where c α = 1/ √ ρ α χ α , is the corresponding Zwiker and Kosten's phase velocity.

Notice that we include the slit width 2s α in the definition of the characteristic admittance, because it simplifies the subsequent writing of continuity conditions; in what follows we replace directly the 2s α by their values Σ, σ, and L -Σ -2l. Now, we start writing the Zwikker and Kosten's equations in the different parts of the periodic cell. In the main tube, we have

-iω ρ t (ω) Σ V t = - ∂P t ∂x + f (4.6a) iωΣχ t (ω)P t = ∂V t ∂x (4.6b)
where, V t = V x Σ is the flow rate field in the tube, with V x the x-component of the velocity in the sense of Zwikker and Kosten (averaged over the section), and P t is the Zwikker and Kosten's pressure in the tube. In the neck, the external excitation having no y-component

iω ρ n (ω) σ V n = ∂P n ∂y (4.7a) iωσχ n (ω)P n = ∂V n ∂y (4.7b)
where, V n = V y σ is the flow rate, with V y the y-component of the velocity, and P n is the Zwikker and Kosten's pressure in the neck. In the cavity

-iω ρ c (ω) L -Σ -2l V c = - ∂P c ∂x + f (4.8a) iω(L -Σ -2l)χ c (ω)P c = ∂V c ∂x (4.8b)
where, V c = V y (L -Σ -2l) is the flow rate and P c the Zwikker and Kosten's pressure in the cavity.

The general solution of the non homogeneous equations (4.6) is written as the sum of the general solution (P t,h , V t,h ) of the homogeneous equations and a particular solution (P t,p , V t,p ) of the non homogeneous equations

P t V t = P t,h V t,h + P t,p V t,p (4.9)
A general solution of the homogeneous equations (4.6) is written as

P t,h V t,h = 1 Y t A + e iktx + 1 -Y t A -e -iktx (4.10)
where A + and A -are the amplitudes of the plane waves in direction of the positive x-axis and negative x-axis, respectively. The following particular solution can be considered

P t,p V t,p = B t C t f 0 e ikx (4.11)
where B t and C t are two constants (for each ω) to be determined. Substituting (4.11) in (4.6) gives the two constants

B t = ikχ -1 t k 2 χ -1 t -ω 2 ρ t (4.12a) C t = iωΣ k 2 χ -1 t -ω 2 ρ t (4.12b)
The particular solution is the same in the left and right portions. On the contrary and because of the presence of the neck, the general solution will have different amplitude constants in the left and right portions. Thus, the general solution of Eq.(4.6) in the main tube can be written as

P t V t = 1 Y t A 1 f 0 e iktx + 1 -Y t A 2 f 0 e -iktx + B t C t f 0 e ikx (4.13a) P t V t = 1 Y t A 3 f 0 e iktx + 1 -Y t A 4 f 0 e -iktx + B t C t f 0 e ikx (4.13b)
where (4.13a) corresponds to the left part of the tube, and (4.13b) to the right part. The constants A 1 , A 2 , A 3 , and A 4 are the amplitude-relating constants to be determined.

The general solution of Eqs.(4.7), (P n , V n ) has the form

P n V n = 1 Y n A 5 f 0 e ikny + 1 -Y n A 6 f 0 e -ikny (4.14)
where A 5 and A 6 are the neck amplitude-relating constants to be determined. Similar to the case of the tube, the general solution of the non homogeneous equations (4.8) is written as the sum of the general solution of the homogeneous equation (P c,h , V c,h ) and a particular solution of the non homogeneous equations (P c,p , V c,p ) in the right or left part of the cavity

P c V c = P c,h V c,h + P c,p V c,p (4.15)
The general solution (P c , V c ) of the homogeneous equation is of the same form as (4.10), both in the left and right part of the cavity but with different amplitudes.

We can find a particular solution as

P c,p V c,p = B c C c f 0 e ikx (4.16)
where B c and C c are two constants to be determined. Substituting (4.16) in (4.8) will give the two constants

B c = ikχ -1 c k 2 χ -1 c -ω 2 ρ c (4.17a) C c = iω(L -Σ -2l) k 2 χ -1 t -ω 2 ρ t (4.17b)
Thus, the general solution of Eq.(4.8) in the cavity can be written as

P c V c = 1 Y c A 7 f 0 e ikcx + 1 -Y c A 8 f 0 e -ikcx + B c C c f 0 e ikx (4.18a) P c V c = 1 Y c A 9 f 0 e ikcx + 1 -Y c A 10 f 0 e -ikcx + B c C c f 0 e ikx (4.18b)
where (4.18a) corresponds to the left part of the cavity, and (4.18b) to the right part. The constants A 7 , A 8 , A 9 and A 10 are the amplitude-relating constants to be determined.

Indeed, in the framework of our simple plane-wave modelling, there are 10 relations concerning the flow rate and pressure, which are assumed to be verified. These continuity relations involve the values of the fields at different locations indicated by numbers (n = 1, ..., 10) in Fig. 4.3. We now proceed to write them.

1-The Bloch condition results in P

4) t = e ikL P (
t , then

A 3 e i k t L 2 + A 4 e -i k t L 2 = e ikL A 1 e -i k t L 2 + A 2 e i k t L 2 (4.19)
t = e ikL V (4)
t , therefore

A 3 e i k t L 2 -A 4 e -i k t L 2 = e ikL A 1 e -i k t L 2 -A 2 e i k t L 2 (4.20)
3-We assume the continuity of the pressure at the junction (2)-(3), P

= P (3) t
t , then

A 3 + A 4 = A 1 + A 2 (4.21)
4-We assume the continuity of the pressure at the junction ( 5)-( 2), P

n = P (5)
t , then

A 5 e -i knl 2 + A 6 e i knl 2 = A 1 + A 2 + B t (4.22)
5-The flow rate at the junction ( 2)-( 3)-( 5) is assumed to verify V

t -V

3) t = V (
n , which yields

Y t A 1 -A 2 -A 3 + A 4 = Y n A 5 e -i knl 2 -A 6 e i knl 2 (4.23)
6-The continuity of the pressure at the junction ( 6)-( 7), P

n = P (7) c results in A 5 e i knl 2 + A 6 e -i knl 2 = A 7 + A 8 + B c (4.24) (6)
7-The flow rate at the junction ( 6)-( 7)-( 8) is assumed to verify V

n +V

n = V (8) c Y n A 5 e i knl 2 -A 6 e -i knl 2 + Y c (A 7 -A 8 ) = Y c (A 9 -A 10 ) (4.25) (7)
8-The pressure is continuous at (7)-( 8) P

= 0, we have

Y c A 7 e -i kc(L-l) 2 -A 8 e i kc(L-l) 2 = -C c e -i k(L-l) 2 (4.27)
10-The flow rate vanishes at the interface solid fluid, V

= 0, we have

Y c A 9 e i kc(L-l) 2 -A 10 e -i kc (L-l) 2 = -C c e -i k(L-l) 2 (4.28)
As such, we have 10 equations (4.19-4.28) for 10 unknowns amplitudes A 1 , .., A 10 . Once these are determined, we will have all the Zwikker and Kosten's fields

Determination of nonlocal effective bulk modulus

Considering the periodic cell (Fig. 4.3), when a harmonic stirring heating Q(t, x) = Q0 e -iωt+ikx = -iωβ 0 T 0 P 0 e -iωt+ikx is applied in the medium, we write the Zwikker and Kosten's equations, in each part of the resonator: tube, neck, and cavity. The aim is to obtain the function χ -1 (ω, k) (see chapter 1)

χ -1 (ω, k) = P(ω, k) + P 0 b ′ (ω, k, x) + φγχ 0 P 0 (4.32)
In the main tube, we write

-iω ρ t (ω) Σ V t = - ∂P t ∂x (4.33a) iω (Σχ t (ω) -γΣχ 0 ) P 0 + iωΣχ t (ω)P t = ∂V t ∂x (4.33b)
The first term of the second equation might not seem to be obvious but follows the very procedure (1.60) seen in nonlocal theory. In the neck, the equations are written as

iω ρ n (ω) σ V n = ∂P n ∂y (4.34a) iω (σχ n (ω) -γσχ 0 ) P 0 e ikx + iωσχ n (ω)P n = ∂V n ∂y (4.34b)
where the term P 0 e ikx comes from the averaging of Q over the neck section. Here also, the second equation might not appear obvious, but follows the procedure (1.60) seen in nonlocal theory. In the cavity

-iω ρ c (ω) L -Σ -2l V c = - ∂P c ∂x (4.35a) iω(L -Σ -2l) [(χ c (ω) -γχ 0 ) P 0 + χ c (ω)P c ] = ∂V c ∂x (4.35b)
The general solution of the non homogeneous equations (4.33) is written as the sum of the general solution of the homogeneous equation (P t,h , V t,h ) and a particular solution of the non homogeneous equations (P t,p , V t,p ) in the right or left part of the tube

P t V t = P t,h V t,h + P t,p V t,p (
B t = ikχ -1 t k 2 χ -1 t -ω 2 ρ t 1 - γχ 0 χ t (4.39a) C t = iωΣ k 2 χ -1 t -ω 2 ρ t 1 - γχ 0 χ t ik (4.39b)
Thus, the general solution Eq. (4.33) in the tube can be written as

P t V t = 1 Y t A 1 P 0 e iktx + 1 -Y t A 2 P 0 e -iktx + B t C t P 0 e ikx (4.40a) P t V t = 1 Y t A 3 P 0 e iktx + 1 -Y t A 4 P 0 e -iktx + B t C t P 0 e ikx (4.40b)
where (4.40a) corresponds to the left part of the tube, and (4.40b) to the right part. The constants A 1 , A 2 , A 3 and A 4 are the amplitude-relating constants to be determined.

As for the tube, the general solution of the non homogeneous equations (4.34) in the neck, is written as the sum of the general solution (P n,h , V n,h ) of the homogeneous equations and a particular solution (P n,p , V n,p ) of the non homogeneous equations

P n V n = P n,h V n,h + P n,p V n,p (4.41)
We can find a particular solution in the following form Thus, the general solution of Eq.(4.34) in the neck can be written as

P n,p V n,p = B n C n P 0 e ikx (
P n V n = 1 Y n A 5 P 0 e ikny + 1 -Y c A 6 P 0 e -ikny + B n 0 P 0 (4.45)
where A 5 and A 6 are amplitude-relating constants to be determined.

In a similar manner, the general solution of Eq.(4.35) is written as the sum of the general solution for homogeneous equation and a particular solution which can be found as

P c,p V c,p = B c C c P 0 e ikx (4.46)
where B c and C c are two constants to be determined. As it has been done before, substituting (4.46) in (4.35) gives the two constants

B t = ikχ -1 c k 2 χ -1 c -ω 2 ρ c 1 - γχ 0 χ c (4.47a) C t = iω(L -Σ -2L) k 2 χ -1 c -ω 2 ρ c 1 - γχ 0 χ c ik (4.47b)
The general solution, then, is expressed as

P c V c = 1 Y c A 7 P 0 e ikcx + 1 -Y c A 8 P 0 e -ikcx + B c C c P 0 e ikx (4.48a) P c V c = 1 Y c A 9 P 0 e ikcx + 1 -Y c A 10 P 0 e -ikcx + B c C c P 0 e ikx (
P (5) n = P (2) t (4.49d) V (2) t -V (3) t
= V (5) n (4.49e) P (6) n = P (7) c (4.49f)

V (6) n + V (7) n = V (8) c (4.49g)

P (7) c = P (8) c (4.49h) V (9) c = 0 (4.49i) V (10) c = 0 (4.49j)
which, respectively, result in the following relations

A 3 e i k t L 2 + A 4 e -i k t L 2 = e ikL A 1 e -i k t L 2 + A 2 e i k t L 2 (4.50a) A 3 e i k t L 2 -A 4 e -i k t L 2 = e ikL A 1 e -i k t L 2 -A 2 e i k t L 2 (4.50b) A 3 + A 4 = A 1 + A 2 (4.50c) A 5 e -i knl 2 + A 6 e i knl 2 + B n = A 1 + A 2 + B t (4.50d) Y t A 1 -A 2 -A 3 + A 4 = Y n A 5 e -i knl 2 -A 6 e i knl 2 (4.50e) A 5 e i knl 2 + A 6 e -i knl 2 + B n = A 7 + A 8 + B c (4.50f) Y n A 5 e i knl 2 -A 6 e -i knl 2 + Y c (A 7 -A 8 ) = Y c (A 9 -A 10 ) (4.50g) A 7 + A 8 = A 9 + A 10 (4.50h) Y c A 7 e -i kc(L-l) 2 -A 8 e i kc(L-l) 2 = -C c e -i k(L-l) 2 (4.50i) Y c A 9 e i kc(L-l) 2 -A 10 e -i kc(L-l) 2 = -C c e -i k(L-l) 2 ( 4

Results

For the geometry considered, the functions ρ(ω, k) and χ -1 (ω, k), with k involved in the excitation terms, can be determined within the approximations in the framework of our modelling. Given these expressions, we know that according to nonlocal theory the possible wavenumbers in the medium will be the solutions of the following dispersion relation as it has been mentioned in previous chapters ρ(ω, q)χ(ω, q)ω 2 = q 2 (4.72) Solving the equation (4.72) by a Newton scheme, we have checked that the obtained expressions of ρ(ω, k) and χ -1 (ω, q) are such that a complex solution q(ω) to (4.72) exists, very close to the value k B (ω) in (4.71).

To perform the computations, we have set L = 1cm, Σ = 0.2L and σ = 0.015L. For the initial values of k(ω), we have chosen the values of k B (ω) with a 10% discrepancy. Fluid properties for all computations are presented in Table 4.1. We see in Fig. 4.4 that the real and imaginary parts of q computed by nonlocal theory via Newton's method converges exactly to the real and imaginary parts of k B which has been computed by a simple Bloch-wave modelling without any use of nonlocal theory. The horizental axis is a reduced frequency equal to ωL/c 0 π, with c 0 the adiabatic sound speed. As such, through nonlocal theory, the 'metamaterial' resonance behaviour of the medium can be perfectly described. T 0 c 0 η κ χ 0 c p γ (kg/m 3 ) (K) (m/s) (kg ms -1 ) (W m -1 K -1 ) (P a -1 ) (J kg -1 K -1 ) 1.205 293.5 340.1391 1.8369 × 10 -5 2.57 × 10 -2 7.173 × 10 -6 997.5422 1.4

Conclusion

The porous matrix with a microgeometry in the form of Two-dimensional arrays of Helmholtz resonators has been considered in this chapter to investigate the validity of the proposed theory of sound propagation through porous media. We have used the homogenization method in nonlocal theory and taking advantage of a plane wave modelling to obtain the effective density and bulk modulus functions in Fourier space. In the framework of the homogenization method, we have employed Zwikker and Kosten's equations governing the pressure and velocity fields' dynamics averaged over the cross-sections of the different parts of Helmholtz resonators, in order to coarse-grain them to the scale of a periodic cell containing one resonator. Once these two effective properties have been determined, the corresponding least attenuated wavenumber of the medium could be obtained through a dispersion equation established via nonlocal theory. The frequency range has been chosen such that the structure-based resonance phenomena could appear.

Indeed, an analytical modelling, has been performed to obtain the least attenuated Bloch mode propagating in the medium. It has been shown that, 120 4 Case of two-dimensional arrays of rigid Helmholtz Resonators the values of Bloch modes obtained in such a way, match exactly those computed by the nonlocal approach. Consequently, we have observed that not only the Bloch wave modelling, but also, especially, the modelling based on the new theory could describe the resonance phenomena, which can be interpreted as a demonstration of the influential effects of the spatial dispersion in the medium. The Finite Element numerical simulations allowing to compute the wavenumbers, in the same manner as has been done in chapter 3, are in progress to confirm the approximations which has been applied in our modelling here.

Conclusion

In this work, inspired by the electromagnetic theory and a particular thermodynamic consideration relating to Umov-Poynting-Heaviside-Schoch, concept of acoustic part of energy current density, we have presented a new nonlocal theory of sound propagation in unbounded homogeneous rigid-frame porous media saturated with a viscothermal fluid. Contrary to the classical local theory, this new approach allows both for temporal and spatial dispersion, which appears in the fact that the acoustic susceptibilities depend both on the frequency and wavenumber. The theory has been formulated to be applied to either isotropic materials, or to periodic materials having a symmetry axis along which the propagation is considered. In the framework of this theory, we have proposed a homogenization procedure to upscale the dynamics of sound propagation from Navier-Stokes-Fourier scale to the volume-average scale, through solving two independent microscopic action-response problems.

An important aspect of the new homogenization method is that, contrary to classical method, there is no lengh-constraint to be considered alongside of its development, thus, in principle, there is no frequency limit for the medium effective properties to be valid. In absence of solid matrix, this procedure leads to Kirchhoff-Langevin's dispersion equation for sound propagation in viscothermal fluids. This theory can be extended, in the future, to anisotropic, bounded media.

The new theory and upscaling procedure has been validated in three cases corresponding to three different microgeometries of the porous structure. A successful test of this theory has been made by a semi-analytical method, in the simple case of cylindrical circular tubes filled with a viscothermal fluid. It has been found that the wavenumbers and impedances predicted by nonlocal theory match with those of the long-known Kirchhoff's exact solution, while, the results by local theory (Zwikker and Kosten's) yield only the wavenumber of the least attenuated mode, in addition, with a small discrepancy compared to Kirchhoff's. 123 124

Conclusion

Another validation case concerned the microgeometry of a two-dimensional array of rigid cylinders. The propagation has been considered along one of the perpendicular axis with which the square lattice can be constructed. The microscopic equations relating to different action-response problems, which lead to the effective density and bulk modulus of the medium according to nonlocal and also local theories, have been solved by direct Finite Element numerical simulations. These have been allowed to compute the corresponding phase velocities of the least attenuated mode according to local and nonlocal approaches. On the other hand, the phase velocity of the least attenuated Bloch wave has been computed through a completely different quasi-exact multiple scattering method taking into account the viscothermal effects. The results of computations based on these three approaches have shown a remarkable agreement between the nonlocal and multiple scattering phase velocity predictions in a wide frequency range. Furthermore, the local theory which takes into account only the temporal dispersion, has shown its limits relating to the frequency, to predict correctly the phase velocity. There, we have observed clearly that the new upscaling procedure, in fact, has imposed no length-constraint; what has been expressed in the correct predictions at high frequency regime. It is conceivable in the future, by improving the performance of the numerical method, to obtain the phase velocities of the higher modes by nonlocal theory and compare them with those obtained by multiple scattering method.

The last case which has been investigated in order to validate the nonocal theory has been related to the microgeometry in the form of a daisy chained Helmholtz resonators. Using the upscaling procedure in nonlocal theory and a plane wave modelling led to two effective density and bulk modulus functions in Fourier space (ω, k). In the framework of this upscaling procedure, we have employed Zwikker and Kosten's equations governing the pressure and velocity fields' dynamics averaged over the cross-sections of the different parts of Helmholtz resonators, in order to coarse-grain them to the scale of a periodic cell containing one resonator. Once these two effective properties have been determined, the corresponding least attenuated wavenumber of the medium could be obtained through a dispersion equation established via nonlocal theory. The frequency range has been chosen such that the structure-based resonance phenomena could appear. Indeed, an analytical modelling, then, has been performed to obtain the least attenuated Bloch mode propagating in the medium. It has been shown that, the values of Bloch modes obtained in such a way, match exactly those computed by the nonlocal approach. Consequently, it has been observed that not only the Bloch wave modelling, but also, especially, the modelling based on the new theory could describe the resonance phenomena, which can be interpreted as a demonstration of the influential effects of the spatial dispersion in the medium. As a matter of We note first that, in addition to Eqs.(A.2), it must also be assumed a relation between the cross-section average of excess temperature and pressure. This relation, between τ and p , would play the role of (A.2a) between V and P.

Using the similarity between (A.1a) and (A.1d) we write it as

ρ ′ ZW (ω) ∂τ ∂t = β 0 T 0 ∂P ∂t (A.5)
We note next that, combining (A.1b) and (A.1c) and averaging over a cross section a general relation between ∂P/∂t, ∂τ /∂t, and ∂V /∂x can be obtained dt ′ ρ ZW (t -t ′ )V (t ′ , x) (A.8c)

H(t, x) = t -∞ dt ′ χ -1 ZW (t -t ′ )B(t ′ , x), H ≡ P (A.8d)
We explain later on in concluding paragraph, why, in this context, the field H coincides with the mean pressure. For the moment, we note that Eqs.(A.2-A.8), local in space, are incompatible with Eqs.(A.1). Indeed, as we have seen, Eqs.(A.1) consistently lead to Maxwellian acoustic equations which are -contrary to the above -nonlocal in time and also in space.

Thus, what Zwikker and Kosten's theory is doing to arrive at equations having the local form (A.2), is not to solve the complete equations (A.1), expressed in the complete action-response problems (2.62-2.64) and (2.77-2.79), but, by introducing various simplifying approximations and idealizations, to solve only some truncated simplified versions of these equations.

The approximations are made in a way to capture the characteristics of the planewave component fields, in the limit where the wavelengths are large compared to the duct transverse dimensions. It is only in this limit that the functions ρ ZW and χ -1 ZW of Zwikker and Kosten's theory allow to describe with high precision the propagation of the least-attenuated, plane wave mode. The simplified versions of the equations, and action-response problems determining ρ ZW and χ -1 ZW , can be directly guessed through the fact that they have to neglect spatial dispersion.

To compute the density, we consider that, since the wavelengths are very large compared to the duct transverse dimensions, the spatial variation of the pressure gradient term in Eq.(A.1a) can be neglected for the purpose of determining the fluid velocity pattern across a section. We thus look at the response of the fluid subjected to the action of an external driving force-per-unit-volume f , which is a pure spatial constant, while harmonic in time (in any real physical wave propagation problem, the temporal variation of the pressure gradient would mean that this gradient is also, to some extent, spatially variable). In this circumstance we have not only to replace the driving f in (2.64) by a spatial constant f 0 e -iωt , but also to drop the response pressure gradient term in (2.62). Indeed, in the cylindrical duct geometry and with constant driving force, no response pressure is generated and no compression-dilatation of the fluid occurs. Thus, we also have to drop the two other fields b and τ . The resulting fictitious problem reads In conclusion, we note that, within the simplifying approximations made in Zwikker and Kosten's local theory, as the excess pressure gradient term in Eq.(A.1a) is represented by the constant term f , and as the excess pressure time derivative term in Eq.(A.1d) is represented by the constant term Q, the excess pressure is replaced by a constant over the cross-section. Thus, when writing the definition up = H u of the H field, the pressure can be extracted from the averaging operation and it turns out that H = p, which also yields H = P since the pressure is constant. We now see why, in the framework of Zwikker and Kosten's local theory, no distinction is to be made between the mean pressure P and the effective macroscopic pressure H. This remark extends more generally to the case of the usual local description [8, Appendix A], which generalizes Zwikker and Kosten's local solution to the case of arbitrary geometries.

Macroscopic theory of sound propagation in rigid-framed porous materials allowing for spatial dispersion: principle and validation

Navid Nemati

Doctoral thesis in Acoustics, Université du Maine, Le Mans, France

This work is dedicated to present and validate a new and generalized macroscopic nonlocal theory of sound propagation in rigid-framed porous media saturated with a viscothermal fluid. This theory allows to go beyond the limits of the classical local theory and within the limits of linear theory, to take not only temporal dispersion, but also spatial dispersion into account. In the framework of the new approach, a homogenization procedure is proposed to upscale the dynamics of sound propagation from Navier-Stokes-Fourier scale to the volume-average scale, through solving two independent microscopic action-response problems. Contrary to the classical method of homogenization, there is no length-constraint to be considered alongside of the development of the new method, thus, there is no frequency limit for the medium effective properties to be valid. In absence of solid matrix, this procedure leads to Kirchhoff-Langevin's dispersion equation for sound propagation in viscothermal fluids.

The new theory and upscaling procedure are validated in three cases corresponding to three different periodic microgeometries of the porous structure. Employing a semianalytical method in the simple case of cylindrical circular tubes filled with a viscothermal fluid, it is found that the wavenumbers and impedances predicted by nonlocal theory match with those of the long-known Kirchhoff's exact solution, while the results by local theory (Zwikker and Kosten's) yield only the wavenumber of the least attenuated mode, in addition, with a small discrepancy compared to Kirchhoff's.

In the case where the porous medium is made of a 2D square network of cylindrical solid inclusions, the frequency-dependent phase velocities of the least attenuated mode are computed based on the local and nonlocal approaches, by using direct Finite Element numerical simulations. The phase velocity of the least attenuated Bloch wave computed through a completely different quasi-exact multiple scattering method taking into account the viscothermal effects, shows a remarkable agreement with those obtained by the nonlocal theory in a wide frequency range.

When the microgeometry is in the form of daisy chained Helmholtz resonators, using the upscaling procedure in nonlocal theory and a plane wave modelling lead to two effective density and bulk modulus functions in Fourier space. In the framework of the new upscaling procedure, Zwikker and Kosten's equations governing the pressure and velocity fields' dynamics averaged over the cross-sections of the different parts of Helmholtz resonators, are employed in order to coarse-grain them to the scale of a periodic cell containing one resonator. The least attenuated wavenumber of the medium is obtained through a dispersion equation established via nonlocal theory, while an analytical modelling is performed, independently, to obtain the least attenuated Bloch mode propagating in the medium, in a frequency range where the resonance phenomena can be observed. The results corresponding to these two different methods show that not only the Bloch wave modelling, but also, especially, the modelling based on the new theory can describe the resonance phenomena originating from the spatial dispersion effects present in the macroscopic dynamics of the matarial.

Introduction 1

 1 Nonlocal theory of sound propagation in porous media 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Electromagnetic recasting of the acoustic equations . . . . . . . . . 1.3 Generalization to fluid-saturated rigid-framed porous media . . . . 1.3.1 Basic definitions . . . . . . . . . . . . . . . . . . . . . .

1. 3 . 4

 34 Identification of the field H It results directly from writing the following macroscopic version of the relation (1.7) used in the viscothermal fluid to identify the field h S = HV (1.43)

Chapter 2 Nonlocal

 2 theory of sound propagation in porous media: case of circular pores Elaborating on a Maxwellian representation of longitudinal wave propagation in a viscothermal fluid, a new general nonlocal macroscopic theory of sound propagation in homogeneous porous media saturated with viscothermal fluid has been recently proposed. The present paper validates this new nonlocal Maxwellian theory by showing that, in the case of the propagation in straight circular tubes, it is in complete agreement with the long known Kirchhoff-Langevin's full solutions.

Figure 3 . 1 :

 31 Figure 3.1: Two-dimensional array of rigid cylinders with identical radius R. The nearest neighbours in this lattice are distanced with the length L. The periodic cell considered is shown by the square of length L.

  3.12) and (3.13) represent the field equations which are completed by two constitutive equations (3.14) and(3.15). The operators ρ and χ-1 are called density and bulk modulus, describing the effective properties of the medium.

  3.1) and bounded in x ∈ [0, L], y ∈ [-L/2, L/2], there are unique amplitude fields v(ω, r) and p(ω, r) solutions to Eqs.(3.21a-3.22), which are periodic with the period L, such that they give the same values on the cell boundaries x = 0 and x = L, ∀y; and also on the cell boundaries y = -L/2 and y = L/2, ∀x.

  3.33) and (3.34) represent the field equations which are completed by two constitutive equations (3.35) and (3.36).

90 3

 3 Case of two-dimensional arrays of rigid cylinders Substituting (3.81) in (3.77), (3.80) for α = a, e, v and β = a, e, v, yields φ αβ n,inc (r, r l ) = ∞ m=0 γ=a,e,v D γβ m C γβ mn G α m (r, r l ) + δ αβ φ β n,in (r, r l ) (3.83)

( 1 -δ 4 L

 14 e, v By identifying the Eqs.(3.75) and the resulting (3.78), we conclude the explicit expressions for the reflection and transmission coefficients R αβ pn and T αβ pn for all βn,R (0, y) cos 2πpy L dy, α = a, e αβ δ pn (1 -δ p0 )e ik α nx L + (1 -δ p0 ) , reflection and transmission properties of one row are entirely determined.

  . The use of scattering matrix relation relation (3.73) and the Bloch condition (3.92) leads to the following eigenvalue problem

94 3 Figure 3 . 3 :

 333 Figure 3.3: Real and imaginary parts of the phase velocities according to local theory, normal mode nonlocal theory and normal mode multiple scattering method for φ = 0.99.

Figure 3 . 4 :

 34 Figure 3.4: Real and imaginary parts of the phase velocities according to local theory, normal mode nonlocal theory and normal mode multiple scattering method for φ = 0.9.

  3),(3.4), and (3.5) referring to the three values of porosity φ = 0.99, φ = 0.9, and φ = 0.7, respectively, real and imaginary parts of phase velocities predicted by local theory, nonlocal theory, and multiple scattering method, are depicted in function of the reduced frequency k 0 L/π = ωL/c 0 π. Regarding the multiple scattering related curves as those showing the most precise values, in all three Figures we observe that the phase velocity predicted by local theory is limited up to a frequency satisfying the condition qL ≪ λ, which has resulted in particular, in the microscale incompressibility of the fluid ∇.v = 0. For all three cases of porosity the rapid variations around reduced frequency k 0 L/π = 1 correspond to the location of the first band gap. This may be viewed as a cell resonance which occurs when the length of the cell is around λ/2.

Figure 3 . 5 :

 35 Figure 3.5: Real and imaginary parts of the phase velocities according to local theory, normal mode nonlocal theory and normal mode multiple scattering method for φ = 0.7.

Figure 4 . 2 :

 42 Figure 4.2: A periodic cell of the structure.

4. 2

 2 Determination of nonlocal effective density 103

Figure 4 . 3 :

 43 Figure 4.3: Illustration of slit portions. Different positions are indicated by (n), and different amplitudes by A n , n = 1, ..., 10.

108 4

 4 Case of two-dimensional arrays of rigid Helmholtz Resonators 2-Again, because of the Bloch condition, V

A 7 + A 8 = A 9 + A 10 ( 4 . 26 ) 9 -

 78104269 The flow rate vanishes at the interface solid-fluid, V

Figure 4 . 4 :

 44 Figure 4.4: Real and imaginary parts of the normal mode, computed by Blochwave calculations, and nonlocal theory.

128 Conclusion

 128 Knowing ρ ZW (ω) and χ ZW (ω) would entirely characterize the tube propagation characteristics. It would give the propagation constantk ZW = ω ρ ZW (ω)χ ZW (ω) (A.3)and the characteristic impedance of the progressive waveZ ZW = ρ ZW (ω)χ -1 ZW (ω) (A.4)

  in this equation the relation (A.5) and using the general thermodynamic identity(2.11), it is easy to verify that, once the functions χ ZW and ρ ′ ZW exist, they must be related byχ ZW (ω) = χ 0 γ -(γ -1) ρ 0 c p ρ ′ ZW (ω) (A.7)Now we observe that Eqs.(A.2) have exactly the same form as the Maxwellian equations (2.58-2.61) with, however, the crucial difference that they are written excluding spatial dispersion

  ξ Z (ω) is the previous relaxation function (A.15) and Pr = ηc p /κ is the Prandtl number. Finally, as ρ ′ ZW is related to the compressibility χ ZW by the relation (A.7), the end Zwikker and Kosten's expression for χ ZW (ω) is written as χ ZW (ω) = χ 0 [1 + (γ -1)ξ(ω Pr)] (A.21)

  Periodic media . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.7 Characteristic wavenumbers and impedances . . . . . . . . 1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Nonlocal theory of sound propagation in a circular pore 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Kirchhoff-Langevin's theory of sound propagation . . . . . . . . . . 2.2.1 Linearized equations, in the Navier-Stokes-Fourier model .

. . 1.3.2 Pore-level equations . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Macroscopic equations . . . . . . . . . . . . . . . . . . . . . 1.3.4 Identification of the field H . . . . . . . . . . . . . . . . . . 1.3.5 Identification of constitutive operators . . . . . . . . . . . . 1.3.6 2.2.2 Propagation in the circular tube . . . . . . . . . . . . . . . 2.3 Application of the general nonlocal theory . . . . . . . . . . . . . . 2.3.1 Linearized macroscopic equations . . . . . . . . . . . . . . . 2.3.2 Determination of the nonlocal density ρ(ω, k) . . . . . . . . 2.3.3 Determination of the nonlocal bulk modulus χ -1 (ω, k) . . . ix x

  26d) Equations (2.26a) and (2.26d) are the linearized Navier-Stokes and Fourier equations, respectively. Eqs.(2.26) are our starting linearized viscothermal equations. For simplicity in what follows the prime over the excess thermodynamic pressure p ′ is omitted.

  propagation of the averaged quantities V and B is described by the following field equations (see Eqs.(1.36), (1.37), (1.41) and (1.42))

	∂B ∂t	+	∂V ∂x	= 0	(2.58)
	∂D ∂t	= -	∂H ∂x	(2.59)
	and constitutive relations

  1 and A 2 , in terms of the arbitrary deriving pressure amplitude P 0 .

			(2.74)
	This assumption leads to Eqs.(1.50) and (1.53) in chapter 1, which give here
	ρ(ω, k) =	k(P + P 0 ) ω u	(2.75)

Knowing the response fields (2.71a) and (2.73) as functions of ω, k and r, is all we need to compute the density ρ(ω, k). According to the conjectured upscaling procedure we assume that in the action-response problem (2.62), the role of the macroscopic field H in Eqs.(2.58-2.61) is played by the field P + P where P is the macroscopic part of the response pressure field p, which is defined by (see chapter 1) pv = P v , i.e. pu = P u

  1 and A 2 , in terms of the arbitrary deriving amplitude P 0 .

Knowing the response fields (2.85a), (2.87) and (2.88) as functions of ω, k and r, is what we need to compute the bulk modulus χ -1 (ω, k). According to the conjectured upscaling procedure in chapter 1, we assume that in the actionresponse problem (2.77), the role of the macroscopic field H in Eqs.(2.58-2.61) is played by the field P + P where P is the macroscopic part of the response pressure field p, which is defined by pv = P v , i.e. pu = P u (2.89) and, at the same time, the role of the macroscopic field B is played by the averaged field b + γχ 0 P . This leads to Eqs.(1.59) and (1.61), which writes here

Table 2 . 1 :

 21 Fluid properties used in all computations.

  1.205 293.5 340.1391 1.8369 × 10 -5 0.6 η 2.57 × 10 -2 7.173 × 10 -6 997.5422 1.4

Table 3 . 1 :

 31 Fluid properties used in all computations.

  + and A -are the amplitudes of the plane waves in direction of the positive x-axis and negative x-axis, respectively. The following particular solution can be considered

	4.3 Determination of nonlocal effective bulk modulus	111
	A general solution of the homogeneous equations (4.33) is written as
	P t,h V t,h	=	1 Y t	A + e iktx +	1 -Y t	A -e -iktx	(4.37)
	where A P t,p V t,p	=	B t C t	P 0 e ikx			(4.38)
	where B t and C t are two constants to be determined. Substituting (4.38) in (4.33)
	gives the two constants
							4.36)

  4.42) where B n and C n are two constants to be determined. Substituting (4.42) in (4.34) gives the two constants To obtain the above expression for B n , the average e ikx has been easily calculated

	B n =	2 kσ	γχ 0 χ n	-1 sin	kσ 2			(4.43a)
	C n = 0							(4.43b)
	e ikx =	1 σ	σ/2 -σ/2	e ikx dx =	2 kσ	sin	kσ 2	(4.44)

  As in the previous section, in the framework of our modelling, there are 10 relations which are assumed to be verified, allowing to relate the flow rate and pressures at different indicated points in Fig.4.3. The assumptions for the flow rate and pressure at different junctions are the same as considered in the previous section:

	P t (4)	= e ikL P t (1)	(4.49a)
	V	(4) t	= e ikL V t (1)	(4.49b)
	P t (3)	= P t (2)	(4.49c)

4.48b)

where (4.48a) corresponds to the left part of the cavity, and (4.48b) to the right part. The constants A 7 , A 8 , A 9 , and A 10 are the amplitude-relating constants to be determined.

Table 4 . 1 :

 41 Fluid properties used in all computations.

	ρ 0

the chapters 2, 3, and 4, chapter 1 has been cited as a reference to this submitted paper.

This chapter in its current form has been recently submitted as a paper to Wave Motion.

Shear motions will be present at the pore scale in the material, but only as a result of longitudinal into shear wave conversion, occurring at the pore walls. This conversion is due to the fact that the longitudinal motions in the bulk fluid cannot satisfy the no-slip conditions at the pore walls, if they keep their initial type of motion.

Nonlocal theory of sound propagation in porous media

Acknowledgements

It gives an immense pleasure to express my gratitude here to people without their help, this work would not be finalized.

Acknowledgements 1 Nonlocal theory of sound propagation in porous media

Bibliography

Chapter 4

Nonlocal theory of sound propagation in porous media; case of two-dimensional arrays of rigid Helmholtz Resonators

Introduction

In this section we will show that by the nonlocal theory presented in chapter 1 [1], we will be able to predict the behaviour of sound propagation in the structures exhibiting the resonance phenomena. To illustrate this we consider -see Fig. 4.1 -a 2D medium with embedded structures of Helmholtz-resonators' type [2]. Macroscopic propagation is considered along the waveguide axis x. In sections 2 and 3, using a Zwikker and Kosten approximation [3], usual in duct acoustics, we will show how to make a simplified modelling of the frequency and wavenumber dependent density ρ(ω, k) and bulk modulus χ -1 (ω, k) of the through Eqs.(4.13), (4.14), and (4.18). At this point, we can easily obtain the cell averages v and pv . Let us start with v regarding the fact that the Zwikker and Kosten's flow rate has no component along the x-axis

Similarly, we can compute pv through the following relation

thereby, we have

Now, we can obtain explicitly the effective density function ρ(ω, k) through Eq.(4.1). In the next section, the effective bulk modulus is computed in a similar way but with a different excitation term, and with exactly the same conditions on the flow rate and pressure fields at different junctions.

can obtain the averages v and pv through the following expressions

Similarly, pv will be computed by

thereby, we have

Normal Bloch modes 115

We need also the expression for b ′ to obtain χ -1 (ω, k). We obtain it in the following way

where n is the normal unit vector outward from the border of integration. In the above, the first line comes from the microscopic mass balance equation. In the second line the divergence theorem has been applied. The integral over the normal component of the microscopic velocity is the difference of the outgoing and ingoing flow rates at the exit and entrance sections of the tube. Now, we can obtain explicitly the effective bulk modulus function χ -1 (ω, k) through Eq.(4.32).

Normal Bloch modes

In this section, we seek the macroscopic Bloch wavenumber k B of the least attenuated wave propagating in the direction of positive x-axis, such that

With a field constituted of 10 Zwikker and Kosten's slit waves, as illustrated in Fig. 4.3, will be associated 10 complex amplitudes A 1 , ..., A 10 . As before, on these 10 amplitudes there are the 2 relations (4.55) expressing the Bloch conditions, and the 8 relations expressing the continuity equations. All these relations are now homogeneous relations, so that nontrivial solutions will be obtained only if the determinant vanishes. This condition will give the wavenumber k B . The first step is to determine the entrance admittance of the resonator Y r = V

(5)

n . The general solution of the homogeneous form of Eqs.(4.8) without the forcing term, is written as

A 2 e -ikcx (4.56a)

where (4.56b) corresponds to the left part of the cavity, and (4.56b) to the right part. Regarding the above equations, the three conditions P (7) c = P (8) c (4.57a)

V (10) c = 0 (4.57c) result in the three following relations

1 + e ikc(L-l) (4.58b)

Using (4.56) combining with (4.58), gives

Then, we can obtain the exprssions for P (6) n and V

n , through the following already indicated continuity conditions P (6) n = P (7) c (4.60a)

which, subsequently, yields the impedance

n and V

Now, we look for the macroscopic wavenumber k B . The following relations are satisfied in the right and left part of the tube

Making use of Eq.(4.55), the above equations result in

On the other hand, as we have seen before, the three following conditions are assumed in the resonator

We have immediately

Writing the two equations resulting from (4.65), and eliminating P

(3) t and P

(2) t in these equations, gives

The determinant of the coefficient matrix must vanish, if the above equations have non-zeros solutions. This yields a second degree algebraic equation

fact, the need for a correct description of spatial dispersion effects, has been the motivation of the new macroscopic theory presented and validated here.

Appendix A

Zwikker and Kosten's simplified local theory

The complete equations of the wave propagation problem in the circular tube are

The aim of Zwikker and Kosten's theory is to find ρ ZW (ω) and χ ZW (ω), such that in harmonic regime

where V and P are the cross-section averages of velocity and pressure.

Conclusion

at r = R, and the deriving force is given by f = -∇P = e x f 0 e -iωt (A.11)

The corresponding density ρ ZW (ω), such that

or, to compare with (2.75)

is the wanted Zwikker and Kosten density. Straightforward calculations result in

where ξ(ω) is the following relaxation function

In a similar manner, to compute the compressibility, we now consider that in the long-wavelength limit, the pressure driving term in Eq.(A.1d) may be viewed as a pure spatial constant for the purpose of determining the excess temperature response pattern, across a section. Thus we consider that this driving term acts as a spatial constant β 0 T 0 ∂p/∂t ≡ Q0 e -iωt with Q0 a constant. The resulting fictitious heat conduction problem reads

for r < R, and τ = 0 (A.17 The corresponding function ρ ′ ZW (ω), such that