Issue 
Acta Acust.
Volume 7, 2023



Article Number  44  
Number of page(s)  29  
Section  General Linear Acoustics  
DOI  https://doi.org/10.1051/aacus/2023034  
Published online  18 September 2023 
Scientific Article
Brownian motion with radioactive decay to calculate the dynamic bulk modulus of gases saturating porous media according to Biot theory
^{1}
Laboratoire d’Acoustique de l’Université du Mans, UMR 6613, 72085 Le Mans, France
^{2}
MSME, Univ Gustave Eiffel, CNRS UMR 8208, Univ Paris Est Creteil, 77454 MarnelaVallée, France
^{*} Corresponding author: denis.lafarge@univlemans.fr
Received:
8
February
2023
Accepted:
5
July
2023
We present a new stochastic simulation method for determining the longwavelength effective dynamic bulk modulus of gases, such as ambient air, saturating porous media with relatively arbitrary microgeometries, i.e., simple enough to warrant Biot’s simplification that the fluid and solid motions are quasiincompressible motions at the pore scale. The simulation method is based on the mathematical isomorphism between two different physical problems. One of them is the actual Fourier heat exchange problem between gas and solid in the context of Biot theory. The other is a diffusiondisintegrationcontrolled problem that considers Brownian motion of diffusing particles undergoing radioactivetype decay in the pore volume and instant decay at the pore walls. By appropriately choosing the decay time and the diffusion coefficient, the stochastic algorithm we develop to determine the average lifetime of the diffusing particles, directly gives the effective apparent modulus of the saturating fluid. We show how it leads to purely geometric stochastic constructions to determine a number of geometrical parameters. After validating the algorithm for cylindrical circular pores, its power is illustrated for the case of fibrous materials of the type used in noise control. The results agree well with a model of the effective modulus with three purely geometric parameters of the pore space: static thermal permeability divided by porosity, static thermal tortuosity, and thermal characteristic length.
Key words: Biot theory / Porous media / Bulk modulus / Brownian motion / Thermal permeability
© The Author(s), Published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
One of the elements of Biot’s macroscopic linear theory of sound propagation [1], as extended for airsaturated porous materials by Attenborough [2] and Allard and Atalla [3], is the complex effective modulus of air, K_{f}(ω), which accounts for heat exchange between the air and the solid constituents taking place at the microscopic level [4] in time harmonic regime. In case where the material can have macroscopic gradients of properties due to inhomogeneous microstructure, a dependence on macroscopic position x will also be present. The index f is for fluid, here the air: recall that in Biot’s theory one has a solid and fluid phase, each forming its own infinite cluster, resp. noted with indexes s and f, and whose coupled macroscopic displacement motions are followed separately [5]. This modulus is of concern in many airborne noise control applications when a precise prediction of sound reflection and transmission versus frequency is needed [3, 4]. Let us define it as the factor relating in harmonic regime e^{−iωt} the excess macroscopic pressure P, to the mean condensation: P/K_{f}(ω, x) = 〈ρ′〉_{p}/ρ_{0}, where 〈ρ′〉_{p} is the porescaleaveraged excess density in the air, and ρ_{0} is the air equilibrium density (ρ_{0} + ρ′ the actual porescale variable air density ρ). See Appendix A and discussion in Section 3 for the definition of totalvolume and porespacevolume averages 〈·〉 and 〈·〉_{p} and the notion of macroscopic pressure. In the time domain, this implies an operator or integral relation of the form, $P(t,\mathit{x})={\widehat{K}}_{f}[\langle \rho \mathrm{\prime}{\rangle}_{p}/{\rho}_{0}](t,\mathit{x})={\int}_{\mathrm{\infty}}^{t}\mathrm{}{K}_{f}(tt\mathrm{\prime},\mathit{x})[\langle \rho \mathrm{\prime}{\rangle}_{p}/{\rho}_{0}](t\mathrm{\prime},\mathit{x})\mathrm{d}t\mathrm{\prime}$, where, using by analogy a language familiar in electromagnetism [6–9], we see the presence of temporal dispersion effects but no spatial dispersion (nonlocalities in time but not in space). Actually, this absence indicates that some drastic simplifications of principle have been made in the Biot model. It is important to realize that the effective bulk modulus kernel, or the operator itself, owes its existence to these simplifications, as do other elements of the theory. Here, we will clearly establish what these simplifications are, and exploit them to devise a very simple method of calculation of the bulk modulus kernel in its Fourier K_{f}(ω) or Laplace ${K}_{f}^{\mathcal{L}}\left(s\right)$ representation. A general theory that would contradict conventional Biot’s theory would have to deal with underlying nonlocal aspects and could lead to very new quantities and operators^{1}. As started in [8, 9] we know that the construction and transition from the spatially local Biot theory to a fully nonlocal generalized theory will not be simple and straightforward. It will have to be the subject of further theoretical and experimental work.
Recall that Biot theory assumes that the acoustic effective wavelengths in the material are large compared to the dimension of averaging volumes. Less wellknown in the literature, it also makes restrictions on the geometries that explain the absence of spatial dispersion. Only those geometries that lead to the characteristic incompressibility of fluid and solid motions, see Appendix A and Sections 2 and 3, are considered. At long wavelengths, as long as the geometries are not too complicate to lead to resonances or localized imbalances originating from other causes (e.g. pressure diffusion process mentioned in Appendix A), the framework of Biot theory will tend to apply and an effective dynamic bulk modulus of the saturating air will tend to exist in the sense of the definition sketched above and developed in Section 3. But as soon as the microgeometries are sufficiently complex to lead to resonances or imbalances incompatible with incompressibility, the macroscopic model under consideration no longer works, and we enter into the realm of the much more complicated nonlocal general description involving both temporal dispersion and spatial dispersion. In this paper, we restrict ourselves to the microgeometric limits of Biot theory when the condition of incompressibility of fluid and solid is satisfied. The reader should keep in mind, however, that there are scenarios (e.g., with Helmholtz resonators [10], but not limited to), which would require a complete redefinition of the treatment as opposed to the simple analysis made here.
In a lowfrequency limit heat exchange between the air and the solid has time to fully establish itself (quasistationary regime). Then, provided that the porous solid structure is thermally inert enough to act as a thermostat (i.e. essentially remain at ambient temperature), the air will have to undergo quasiisothermal expansion and compression [11]. Therefore in this limit K_{f}(ω) will tend to the air isothermal bulk modulus K_{0} (the atmospheric pressure).
Conversely, in a highfrequency limit, heat exchange does not has time to occur, except in the immediate vicinity of the pore walls (inside a socalled thermal boundary layer of characteristic thickness (2ν′/ω)^{1/2}, where ν′ = κ/ρ_{0}c_{P} is the thermal diffusivity [4], κ the thermal conductivity and c_{P} the specific heat coefficient at constant presssure, of air). The air will undergo adiabatic expansion and compression except in this evershrinking vicinity. In this limit K_{f}(ω) tends to the adiabatic bulk modulus K_{a} = γK_{0} (with γ = c_{P}/c_{V} ≈ 1.4, the ratio of specific heat coefficients of air at constant pressure or constant volume). In between, at intermediate frequencies, a smooth relaxation transition is observed, which connects the “relaxed” state (lowfrequency or “static” isothermal limit), to the “frozen” state (highfrequency adiabatic limit), and is determined by the geometry and dimensions of the pore space and the thermal diffusivity ν′ of air^{2}.
More generally, within Biot theory and when the saturating fluid is not air but any trivariate fluid [12] with γ significantly departing from 1, meaning a gas as opposed to a liquid, the effective fluid bulk modulus is in principle entirely determined by the solution of a “Fourier”^{3} diffusiondisintegrationcontrolled problem^{4}, which can be stated and solved in the pore space – see Appendix in [9] and the comprehensive presentation in Section 3, leading to equations (33)–(37), equations (29)–(30) and equations (42)–(43). In general, because of the complications in the geometries, no closedform analytic solution can be given, and numerical simulations such as finite differences or finite elements have been used to determine the precise solution.
In practice the geometries present random variations and a given type of material can be best represented by an ensemble of realizations. Thus, the numerical process of obtaining an estimate of K_{f}(ω) can become significantly cumbersome: first, the local governing equations for the Fourier diffusiondisintegration problem with appropriate boundary conditions have to be solved using numerical techniques such as mentioned; this will be repeated in an ensemble of configurations; and then the results will be configurationally averaged as the effective bulk modulus will rely averaged fields^{5}. As noted by Torquato and Kim [13] in a larger context, this can be a wasteful calculation process as there is a significant amount of information that is lost in passing from the local to the averaged fields. To quote them: “such calculations can become exorbitant even when performed on a supercomputer”.
In this paper, we present a meshless stochastic numerical algorithm to perform the above operations and compute – with a possible shunting of the local field determination step – the effective modulus K_{f}(ω) that appears in Biot theory. The idea of the algorithm directly inspired by Torquato and Kim [13] is not new, see [14], but was given in the latter Ref. with one faulty reasoning (and a resulting faulty Eq. (9)). After due correction of the latter, it was, through personal communications, applied by Perrot et al. [15] to the description of thermal dissipation properties of threedimensional reconstructed unit cells of aluminum foams saturated by the ambient air. Here we describe precisely the Biot context of the effective fluid bulkmodulus model and present in detail the (corrected) calculation procedure^{6}. We give a clear validation on cylindrical circular tubes and redo properly the simple illustration cases of Ref. [14].
The paper is organized as follows. First, we analyze in Section 2 and Appendices A and B the physical position of Biot theory. We highlight the role played by the fluid/solid incompressibility characteristics, their relation to hypotheses on microgeometry, and the nature of the consequences. In Section 3, we analyze how the notions of dynamicviscous and dynamicthermal tortuosities/permeabilities response functions emerge within Biot’s simplifications, and show how the thermal ones determine the effective fluid bulkmodulus. In the remaining sections, we then establish how efficiently these thermal response functions can be calculated using the mean survival times of particles undergoing Brownian motion with radioactive decay in the volume and instantaneous absorption at the pore walls. To this end, in Section 4 we generalize a series of considerations of Lifson and Jackson [16], who develop a method originally due to Pontrjagin et al. [17], leading to a rigorous and simple derivation of the differential equation obeyed by the above survival times. In Section 5 we use the results from Section 4 and derive the relevant statistical properties of Brownian motion trajectories with radioactive decay, which appear in the trajectory ensemble construction procedure proposed by Torquato and Kim [13]. It is finally revealed in the proposed algorithm, which provides a stochastic purely geometric way to obtain the Fourier or Laplace dynamic thermal tortuosity/thermal permeability and the dynamic bulk modulus. We obtain simple pure geometric procedures to calculate some pure geometric parameters embedded in these functions and appearing in the lowfrequency limit, such as static thermal permeability and static thermal tortuosity. Higher order parameters could be obtained sequentially. Also, we show how the parameters appearing in the opposite highfrequency limit can be computed in another recursive process taking a limit ω → ∞. To validate our simulation method we calculate in Section 6 the effective bulk modulus for the case of parallel cylindrical circular pores for which there are wellknown closedform analytical results (Zwikker and Kosten’s). The power of the method is illustrated by deriving the effective dynamic bulk modulus of ambient air saturating a, regular or random, idealized simple fibrous material of the type used in noise control applications. We show that in both regular and random cases the results support a longproposed frequencydependent model [18], recently reexposed ([9], Appendix), that includes geometric parameters, in addition to porosity ϕ, called (static) thermal permeability ${k}_{0}^{\text{'}}$ (the inverse trapping constant [19], Lafarge [18]) (static) thermal tortuosity ${\alpha}_{0}^{\text{'}}$ (Lafarge [18]), and thermal characteristic length Λ′ (Allard 1991 [20]). Finally, we give conclusions in Section 7.
2 Special physical position of Biot theory and limited scope of the present work
Biot’s (1956) classical linear theory [1], complemented to account for thermal effects [2, 3] and losses in the solid [3], is a longwavelength semiempirical model allowing to describe lowamplitude sound propagation in many different porous elastic solids saturated with compressible viscothermal fluids. For a long time, it was not clearly realized that Biot theory relies on tacit limitations made on microstructure, leading to the incompressibility of fluid and solid motions on a small scale, i.e. in a representative averagingvolume, or in other words, on a simplified view of dispersion effects that completely ignores spatial dispersion and simplifies temporal dispersion.
The model was originally developed in the context of geophysical applications where one is concerned with porous media (rocks) saturated by heavy fluids (oil). This model describes a situation where the fluid and the skeleton, each forming its own infinite cluster connected through the sample, move simultaneously at different macroscopic velocities, while microscopically their motions – significantly distributed at the pore scale for the fluid because of the tortuous microgeometry, but not for the solid as it cannot slide like a fluid – are nearly divergencefree due to long wavelengths. To our knowledge, this hypothesis has only recently been clearly identified as a constraint on the allowable microgeometries in conventional Biot theory [8, 9].
In the late 1970s and in the 1980s, important works contributed, if not a fully clarified understanding, at least a better knowledge of Biot theory. It became established as a relatively general theory of poroelastic fluidsaturated media, which describes the coupled solid and fluid motions, and expresses in coherent manner the inertial, viscous, and elastic interactions between the two phases^{7}. It was extended and used by Attenborough [2], Allard and collaborators [20–26], and by others, e.g., Bolton and collaborators [27], as the basis for a general description of sound propagation in many common porous acoustic materials saturated by ambient air (polyurethane foams, glass wool, rock wool, etc.), used either singly or in multiple layers with other materials to effectively absorb sound, in noise control applications [3].
Here, we wish to highlight the cornerstone importance, in Biot theory, of the simplification of “quasiincompressibility” (in a representative volume) mentioned in Abstract and Section 1 (more shortly referred to as the incompressibility hypothesis, or incompressibility, below). Once again, it means quasi divergencefree motion for the fluid and quasi uniform motion of the solid, at the pore scale, i.e. in a representative averagingvolume.
As described by Johnson et al. [5, 28] in the geophysical heavy fluid context, all the detailed physics of the system microgeometry is hidden in a few parameters; Allard and collaborators complete the list to describe the frequencydependent bulk modulus that arises from thermal exchanges when the fluid is a gas. In the lowfrequency limit where Biot theory was first formulated, the important parameters are porosity ϕ, Darcy permeability k_{0}^{8}, the static tortuosity α_{0} [9, 18, 29], and the abstract BiotWillis elastic constants P, Q, R, N [5, 30]. The latter elastic constants are identified through the three quasistatic BiotWillis Gedanken Experiments [30]. From the first Gedanken Experiment where the material is subjected to pure shear, N is identified as the shear modulus of the frame itself unaffected by the presence of the viscous saturating fluid (no restoring shear force arises between the frame and the fluid when the latter is macroscopically sheared). From the other two Gedanken Experiments – one in which a representative volume, jacketed and drained with small capillaries, is subjected to hydrostatic overpressure, and one in which, without a jacket, it is directly subjected to the same hydrostatic overpressure – the remaining macroscopic coefficients P, Q, R are expressed by the socalled BiotWillis relations in terms of: ϕ, N, and K_{b}, K_{s}, K_{f} (see e.g. Eqs. (2.11 a, b, c) in [5] or Eqs. (6.18, 19, 20) in [3]), where K_{b} is the bulk modulus of the frame in vacuum, K_{s} is the bulk modulus of the elastic solid from which the frame is made, and K_{f} is the bulk modulus of the saturating fluid itself.
Starting from the quasistatic limit, the theory directly generalizes in arbitrary harmonic regime by using the cornerstone incompressibility hypothesis: To our knowledge, this has never been noted before in the literature, however, it is through this hypothesis that one can consider extending to arbitrary frequencies the quasistatic thought experiments. The BiotWillis relations, giving the elastic Biot coefficients, remain valid in the harmonic regime by using the complex values for the elastic constants arising from the viscoelastic nature of the solid phase, and, when the fluid is a gas, the heat exchanges. The effect of internal friction in solids is accounted for by substituting simple frequencydependent values (in practice complex amplitude values with socalled small loss angles that are approximately constant), for the original real coefficients N, K_{b} and K_{s}. The effect of thermal exchanges between solid and gas is accounted for by substituting the complex function K_{f}(ω) for the original real K_{f}. The fact that the incompressibility hypothesis is the reason for the transfer of the static BiotWillis experiments into the harmonic domain with the abovementioned substitutions makes the importance of this hypothesis clear.
At arbitrary frequencies it has been described how the parameters k_{0} and α_{0} combine in a single notion, the frequencydependent dynamic permeability/tortuosity function, k(ω) and α(ω) (see [18] and Appendix in [9]), wellstudied in Johnson et al.’s landmark works, e.g. [28] (α(ω) ≡ νϕ/(−iωk(ω)), ν = η/ρ_{0}, with η dynamic viscosity, ν kinematic viscosity). In particular, Johnson et al., 1986, clarified the highfrequency limit properties of these functions which gave rise to the classical [31] notion of ideal fluid tortuosity α_{∞}, and, when the pore wall surface is smooth, to the important notion of characteristic viscous or electrical length Λ [28, 32]. As explained in [28], the dynamic tortuosity α(ω) directly determines the various effective densities (ρ_{11}, ρ_{22}, ρ_{12}), that appear in Biot theory. It is a complex frequencydependent function to account in detail for the inertialviscous interaction between the solid and the fluid. Explicitly, the densities verify, ρ_{11} + ρ_{12} = (1 − ϕ)ρ_{s}, ρ_{22} + ρ_{12} = ϕρ_{0}, and, ρ_{12} = −(α(ω) − 1)ϕρ_{0}, where ρ_{s} is the density of the elastic solid on which the frame is built; ρ_{12} describes the inertialviscous coupling coefficient for the drag that the fluid exerts on the solid, and vice versa by the principle of action and reaction. It arises as both fluid and solid constituents are in different macroscopic motions, one moving and accelerating with respect to the other, which generates viscous and inertial reaction forces at the pore walls. Johnson has explained how these functions α(ω) and k(ω) are strongly constrained by the condition that the fluid motion at the pore scale is divergencefree to a first approximation in the long wavelength limit inherent to Biot theory, and as a consequence, are very well represented in terms of a few lowfrequency and highfrequency geometric parameters (see Appendix A in [28] and also, Appendix in [9]). Actually, the above relations on the Biot effective densities ρ_{11}, ρ_{12}, ρ_{22,} are obtained on the condition that, not only the fluid motion is divergencefree, but also the solid moves entirely as a whole on the pore scale, in first approximation. Thus, we have to say that all Biot densities are defined and strongly constrained by the simultaneous conditions of fluid and solid incompressibility. It is this simplification that solid and fluid motions are locally divergencefree whether or not the frame oscillations characteristic of Biot theory manifest themselves strongly, which explains that there is a unique set of frequencydependent effective densities that must be considered in the theory. They derive from the unique notion of the dynamic permeability/tortuosity response function. The latter can be studied for a rigidframe motionless structure of same microgeometry as the actual, undeformed, material.
Like the densities and for the same reasons the function K_{f}(ω) is unique, whether the frame oscillations characteristic of Biot theory manifest themselves or not. It can be likewise studied for a rigidframe motionless structure. It needs to be considered only when the saturating fluid is not a liquid, as in the original Biot considerations, but a fluid with nonnegligible thermal expansion. This is the case of gases in general, such as air here. Gases have an expansioncompression behavior well described by the perfect gas law, and that means they have a coefficient of thermal expansion, β_{0} ≡ −[∂ρ/(ρ∂P)]P_{0}, very close to 1/T_{0} (T_{0}, absolute ambient temperature, P_{0}, atmospheric ambient pressure, ρ, the density). By the general thermodynamic relation, $\gamma 1={T}_{0}{\beta}_{0}^{2}{c}_{0}^{2}/{c}_{P}$, (${c}_{0}=\sqrt{{K}_{a}/{\rho}_{0}}$, adiabatic speed of sound), that can be found e.g. in Landau and Lifshitz [33], it means that the factor γ − 1 deviates significantly from zero (general trivariate fluid [12]). For liquids, this is not the case because β_{0} is very small and, as the thermodynamic identity shows, the deviation is a quadratic effect on β_{0}. Simple kinetic theory considerations give: γ – 1 = 2/n_{d}, which is not small, where ${n}_{d}$ is the number of excited degrees of freedom of a typical molecule (for air, mainly composed of diatomic molecules, n_{d} ≈ 5 and γ − 1 ≈ 0.4). In this case, the isothermal K_{0} and adiabatic K_{a} = γK_{0} bulk moduli differ significantly and there is a frequency dependence of the effective fluid bulk modulus that must be taken into account as already explained in the Introduction.
The frequency dependence of the effective fluid bulk modulus is directly obtained through the introduction of thermal dynamic permeability/tortuosity functions k′(ω) and α′(ω), kind of thermal counterparts of the notions of viscous dynamic permeability/tortuosity [4, 9, 18]. Johnson et al.’s entire analysis of the kind of frequency dependence that appears in these functions [28] translates directly for them, e.g. [4] (Appendix C) and [9] (Appendix). It is the model of K_{f}(ω) [9, 18] constructed using parameters k′_{0}/ϕ, α′_{0}, Λ′, and based on these general considerations, that will be checked in our numerical simulations in Section 6.
Historically, the lowfrequency limit was first clarified for the case of heavy fluids (part I of Biot’s 1956 series of two papers [1]), and the theory generalized in harmonic regime with a description of inertialviscous effects thanks to the notion of dynamic viscous permeability/tortuosity (Johnson et al. paper [28]). With the work of Allard and collaborators and the consideration of additional loss phenomena associated with internal friction in solids and heat exchange between fluids and solids, the theory finally extended its already tremendous predictive power, reported in [5]. The theory was found to be particularly useful in noise control applications, see [3].
This procedure allowed one to express consistently the temporal dispersion that appears when the basic fluid/solid incompressibility assumption of Biot theory is respected. In the following, we will base the description of the function K_{f}(ω) directly on incompressibility, but we should not forget that this does not hold in the general case. We refer to the Appendices A and B for additional comments on the physical position of Biot theory and the present study.
3 Definition of the dynamic bulkmodulus and related functions
For the class of materials well described by Biot theory, the acoustic vibrations of the poroelastic structures have no effect on the apparent bulk modulus K_{f}(ω) of the saturation fluid. In what follows, it is assumed without loss of generality that the porous skeleton is sufficiently inert/rigid not to be set in motion. When modeling the geometrically constrained effective bulk modulus of air in the context of Biot’s “local theory”, this makes it possible to simplify the considerations and calculations in a meaningful way.
We exclude microgeometries that would include structures such as Helmholtz resonators: Their resonances would trigger the spatial dispersion phenomena that lie outside the conventional Biot description (see [9], previous Section 2, and Appendices A and B). If we were to perform a general nonlocal homogenization that accounts for all types of microgeometries and is thus able to account for the most general spatial and temporal dispersion, we would have to start at the pore scale with the full set of governing coupled equations of sound propagation. In the present case of an inert solid frame, these would be equations (7.16) from Ref. [9]. In these equations, the zero velocity boundary condition is set at the pore walls because the porous skeleton is assumed to be immobile; the zero excess temperature boundary condition is set because we also assume it to be thermally inert [4] (Appendix B.3).
However, for the Biot description we are interested in here, since we assume incompressibility of the fluid at the pore scale, we must replace this set of coupled equations for inertia, viscosity, elasticity and heat, with a simplified set. This is a set of simpler and disconnected parts, one describing the inertia and viscosity effects, leading to the notion of dynamic viscous permeability and tortuosity, one describing heat exchange, leading to the notion of dynamic thermal permeability and tortuosity, and one expressing only the equation of state from which the relationship between the latter thermal terms and the effective bulk modulus of the fluid is derived. The principle of this simplification has been studied, e.g. in [4] and [9] (Appendix). We revisit it in detail here and take the opportunity to further elaborate on what was said earlier about the basic assumption of incompressibility of the theory and the relation that this assumption has with the rejection of the description of spatial dispersion and a related part of temporal dispersion. We think that this will be useful in understanding the position of the Biot theory that we stated in the Section 2 and the Appendices A and B.
First, we will examine the simplified equations, detached from the other equations, that we can use to describe the inertia/viscosity effects. Let us direct our attention to a finite, coarsegrained spherical region centered at a certain position x within the material and whose radius is chosen large enough to obtain smoothed spatial averages, but not so large as to cancel the actual pointtopoint changes that can be meaningfully observed – see the general remarks of Lorentz [34] pp. 133–134, §§ 113 and 114. That is, in addition to a longwavelength limit implying a scale separation between the (small) radius of the averaging sphere and the (large) characteristic distance over which the acoustic macroscopic fields vary, we also assume a scale separation between this homogenization radius and the characteristic distance over which the macroscopic parameters of the material can vary^{9}. Under these conditions, the internal structure inside the sphere or averaging ball (b) (where “(b)” stands for “ball” or “bounded”, since its important property is mainly its bounded spatial extent^{10}) can be considered as that of a macroscopically homogeneous medium. We will approximate it as a finite part of a periodic or stationary random medium representing the structure of the material around the considered macroscopic position.
Let this solidfluid averaging ball (b) be a bounded spherical region with a connected pore space ${V}_{p}^{\left(b\right)}$ (filled by the fluid) and a solidfluid pore wall surface ${S}_{p}^{\left(b\right)}$. Underlying Biot’s local theoretical approach is the notion that the fluid velocity pattern v(t, r) in this region of central position x and thus its macroscopic mean 〈v〉(t, x), which is the quantity of interest to us, can be viewed as responding to the history of macroscopic pressure gradients applied (in the present and in the past) at this very position^{11}.
Before proceeding, let us distinguish two possible manifestations of the spatial derivative operator, $\nabla =\widehat{\mathbf{x}}\partial /\partial x+\widehat{\mathbf{y}}\partial /\partial y+\widehat{\mathbf{z}}\partial /\partial z$, which is also written in a classical shorthand notation, ∇ = ∂/∂x (e.g., Landau and Lifshitz, Mechanics [35], § 5, The Lagrangian for a system of particles), where $\mathbf{x}=x\widehat{\mathbf{x}}+y\widehat{\mathbf{y}}+z\widehat{\mathbf{z}}$. For the sake of clarity we adopt two different notations for a position vector x: r for a microscopic position vector indicating a position in the pore space, and x for a macroscopic position vector referring to a macroscopic quantity. With this in mind, in what follows, ∇ will mean ∂/∂r when acting on a microscopic variable, and ∇ = ∂/∂x when acting on a coarsegrained variable such as P or 〈v〉, where x denotes the “central position” of the averaging sphere that is used to define the macroscopic quantity in question (with, for the specific case of macroscopic pressure in the fluid, either of the definitions that we mentioned in Appendix A, ϕP = 〈p〉 or 〈pv〉 = P〈v〉).
Since the actual microscopic air motion in the spherical spatial averaging region satisfies, at any time, the equation,
$${\rho}_{0}\frac{\mathrm{\partial}\mathit{v}}{\mathrm{\partial}t}=\nabla p+\eta \Delta \mathit{v}+\left(\frac{\eta}{3}+\zeta \right)\nabla \left(\nabla \xb7\mathit{v}\right),\hspace{1em}\mathrm{in}{V}_{p}^{\left(b\right)},$$(1)
(which is one among the complete set of governing coupled equations mentioned above (7.16) in [9], η and ζ, first and second viscosities, and p a thermodynamic variable as explained in [8]), with boundary conditions:
$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on}{S}_{p}^{\left(b\right)},$$(2)
and since in Biot theory we apply a divergencefree simplification of the representation of air motion at the pore scale, we must conclude that in this framework we are considering only a situation in which the pattern of air velocity occurring in the spherical averaging region must very nearly satisfy equations having the following form, which is drastically simplified:
$${\rho}_{0}\frac{\mathrm{\partial}\mathit{v}}{\mathrm{\partial}t}=\nabla p+\eta \Delta \mathit{v},\hspace{1em}\mathrm{in}{V}_{p}^{\left(b\right)},$$(3)
$$\nabla \xb7\mathit{v}=0,\hspace{1em}\mathrm{in}{V}_{p}^{\left(b\right)},$$(4)
with boundary conditions:
$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on}{S}_{p}^{\left(b\right)}.$$(5)
Because of the lack of definition of variable p in (3) (which is not the same variable as p in (1), and in particular is no longer a thermodynamic variable), and boundary conditions on the remaining surfaces bounding the averaging sphere in the fluid, these equations do not correspond to a welldefined problem. However, a welldefined interpretation/solution of these equations that makes sense as a representation of the velocity field pattern inside the sphere and is consistent with the above view of a velocity pattern that is a response to the time evolution of the applied macroscopic pressure gradients can be achieved by completing the following.
We will imagine that we extend the regions ${V}_{p}^{\left(b\right)}$ and ${S}_{p}^{\left(b\right)}$ into the space outside the sphere in an appropriately defined “uniform manner”, as if the material were macroscopically homogeneous. Using the periodic or stationaryrandom approximation seen above, we will think of the periodic extension or stationaryrandom extension of the internal structure over the entire space, which generates the regions denoted below ${V}_{p}^{\left(\mathrm{ub}\right)}$ and ${S}_{p}^{\left(\mathrm{ub}\right)}$, where “(ub)” stands for “unbounded”. When the material is macroscopically inhomogeneous, these regions should be called ${V}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and ${S}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ because their construction refers specifically to the macroscopic position x at which a certain type of local structure of the material is assumed. By construction, the regions ${V}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and ${S}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ are macroscopically homogeneous in nature^{12}.
Let −∇P(t′, x) for t′ ∈ [−∞, t] define the values that the physical macroscopic pressure gradients have taken at the x position, and recall that we can define these gradients using either of the definitions of macroscopic pressure that we mentioned in Appendix A. Stated and solved in the construct unrestricted domains ${V}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and ${S}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, the equations (3)–(5) have a unique abstract mathematical solution for the velocity v(t, r), that will be a valid representation of the physical velocities in a representative volume around x, if we constrain the abstract pressure, p(t′, r) = P + π, to be split in a uniformly variable abstract macroscopic part P defined from the values −∇P(t′, x) by: P(t′, r) = (r − x) · ∇P(t′, x) + P(t′, x) (note that the knowledge of −∇P(t′, x) does not give that of P(t′, x) but the latter value will play no role here), and a deviatoric part, π(t′, r), that will automatically be stationaryrandom or periodic as a function of r, to be coherent with the spatial constancy of the gradient of the abovedefined abstract P(t′, r), (−∂P(t′, r)/∂r = −∇P(t′, x)), which constancy is itself the result of the incompressibility of the considered motion. Obviously this condition that π must be stationaryrandom or periodic must hold in same manner whether we take, on the abstract field p or abstract fields p, v, the definition ϕP = 〈p〉 or the definition P〈v〉 = 〈pv〉, to obtain our representative macroscopic pressure P. Finally, with one or the other definition, we solve the same problem where the spatially uniform, arbitrarily time variable forcing −∇P, is given by specified values. With one or the other definition, the equations of the problem are at all times t′ ∈ [−∞, t]:
$${\rho}_{0}\frac{\mathrm{\partial}\mathit{v}}{\mathrm{\partial}t}=\nabla \pi +\eta \Delta \mathit{v}\nabla P,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(6)
$$\nabla \xb7\mathit{v}=0,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(7)
$$\pi :\mathrm{stationary}\mathrm{random}/\mathrm{periodic},\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(8)
with boundary conditions:
$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on}{S}_{p}^{\left(\mathrm{ub}\right)},$$(9)
with only the difference that the stationaryrandom π(t′, r) solution, that is unique up to an additive integration constant C(t′), can be fixed in one or the other case either with 〈π〉 = 0 or with 〈πv〉 = 0. However, as these two ways to remove the additive constant are obviously equivalent we must conclude that the same solution, v(t′, r), π(t′, r), is obtained, and thus the same macroscopic pressure P is obtained, in both the direct and indirect averagingconceptions of Appendix A^{13}. We note also that, the spatial constancy of the considered macroscopic pressure gradient is a signature that there can be no account of spatial dispersion in the description to be developed. Since this constancy is directly related to the incompressibility condition in the posed mathematical problem, it shows the already mentioned connection between the assumption of incompressibility at the pore level and the rejection of spatial dispersion at macroscopic level.
As for the time variations of the macroscopic pressure gradient, they can be assumed to be arbitrary (as long as this is compatible with the long wavelength considerations), i.e., $\nabla P=f\left(t\right)\widehat{\mathit{e}}\left(t\right)$, with arbitrary amplitude function f(t) and arbitrary timedependent orientation of the unit vector $\widehat{\mathit{e}}$, ${\widehat{\mathit{e}}}^{2}=1$. It should be clear that, if we allow for arbitrary time variation here, this in no way means that the effects of temporal dispersion fully appear in the description. Rather, as mentioned in Section 2, the allowed temporal dispersion effect is very strictly limited. The use of the incompressibility assumption not only precludes the occurrence of spatial dispersion effects, but also severely restricts the temporal dispersion effects: as previously mentioned in Section 2 it causes the zeros and singularities of the response functions in the complex frequency plane to lie on the imaginary halfaxis [28] (Appendix A), as opposed to the full halfplane that is allowed in pure causality considerations [6, 36]. This becomes clear in Avellaneda and Torquato’s, 1991, relaxationtimes representation of the viscous dynamical permeability [19]: their relaxation times Θ_{n} correspond to purely imaginary singular angular frequencies, ω_{n} = −i/Θ_{n}. It follows that these functions are constructed with basis elements which have strictly monotone behavior on the real frequency axis (see [18, 37] and Appendix in [9]).
Note that, exactly the same response fields, v and π, will appear in another abstract problem extended as above:
$${\rho}_{0}\frac{\mathrm{\partial}\mathit{v}}{\mathrm{\partial}t}=\nabla \pi +\eta \Delta \mathit{v}+\mathit{f},\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(10)
$$\nabla \xb7\mathit{v}=0,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(11)
$$\pi :\mathrm{stationary}\mathrm{random}/\mathrm{periodic},\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(12)
with boundary conditions:
$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on}{S}_{p}^{\left(\mathrm{ub}\right)},$$(13)
and in which the imposed spatial constant vector f, instead of representing a macroscopic pressure gradient, now represents an imposed external, bulk, body force, directly acting in the fluid. The same motion as before will result when this spatial constant force will be taken, at the different instants t′ ∈ [−∞, t], numerically equal to the previously considered pressure gradient: $\mathit{f}\left(t\mathrm{\prime}\right)=f\left(t\mathrm{\prime}\right)\widehat{\mathit{e}}\left(t\mathrm{\prime}\right)=\nabla P(t\mathrm{\prime},\mathit{x})$. In this second problem, the quantity π, in a way, represents the pressure itself, which appears in the fluid in the given nontrivial microgeometry, in reaction to the imposed, past and present, forces f. In the preceding problem, π, in a way, was the local deviation to the macroscopic excess pressure. We say “in a way”, because, in these mathematical problems, since the fluid is viewed as completely incompressible which is a purely mental view, the quantities that seem to act as actual physical overpressures (p = P + π in the first problem, p = π in the second) have only symbolic meaning, not physical meaning (remember that the thermodynamic meaning has been lost). We cannot expect these abstract mathematical quantities to yield a truly quantitatively meaningful representation of the physical excess pressure in the actual physical problem of sound propagation. Recall that the unphysical spatial uniformity of the source terms, related to the divergencefree nature of the motion, is an expression of the fact that the spatial dispersion effects (intrinsically connected with the spatial variations of source fields [6]), are, by force here, entirely left outside the description framework. This is why the π, small, abstract pressure variations, will necessarily be truncated representations of the physical ones; they will be sufficient for our purpose of accurately representing the velocity pattern and computing its spatial average 〈v〉 in a representative volume but we believe that they will not very accurately capture the pattern of the pressure^{14}.
With arbitrary time variations of the source terms (the macroscopic pressure gradient or the external bulk force), the solution is a general convolution of the elementary solution of the impulse Stokes problem considered in [19] (response of the viscous incompressible fluid to a temporal Diracdelta stirring force, set along a given direction $\widehat{\mathit{e}}$ in the stationary random or periodic medium). Of course, this solution exhibits the limitations indicated in [28] (Appendix A), in that it is constructed with purely damped normal modes^{15}. By explicitly writing the convolution we will learn that there are local “tortuosity” macroscopic operators ${\widehat{\alpha}}_{\mathrm{ij}}\left(\mathit{x}\right)$, or response factors α_{ij}(t, x), such that,
$${\rho}_{0}\frac{\mathrm{\partial}\langle {\mathit{v}}_{i}{\rangle}_{p}}{\mathrm{\partial}t}\left(t\right)={\widehat{\alpha}}_{\mathrm{ij}}^{1}\left(\mathit{x}\right){\nabla}_{j}P\left(t\right)={\int}_{\mathrm{\infty}}^{t}\mathrm{}\mathrm{d}{t}^{\mathrm{\prime}}{\alpha}_{\mathrm{ij}}^{1}\left(t{t}^{\mathrm{\prime}},\mathit{x}\right){\nabla}_{j}P\left({t}^{\mathrm{\prime}}\right),$$(14)
or equivalently, local “permeability” macroscopic operators ${\widehat{k}}_{\mathrm{ij}}\left(\mathit{x}\right)$, or response factors k_{ij}(t, x), such that,
$$\varphi \left(\mathit{x}\right)\langle {\mathit{v}}_{i}{\rangle}_{p}\left(t\right)=\frac{{\widehat{k}}_{\mathrm{ij}}\left(\mathit{x}\right)}{\eta}{\nabla}_{j}P\left(t\right)=\frac{1}{\eta}{\int}_{\mathrm{\infty}}^{t}\mathrm{}\mathrm{d}{t}^{\mathrm{\prime}}{k}_{\mathrm{ij}}(t{t}^{\mathrm{\prime}},\mathit{x}){\nabla}_{j}P\left({t}^{\mathrm{\prime}}\right),$$(15)
(with symmetry on the exchange of i = 1, 2, 3 and j = 1, 2, 3, see below; the general case of an anisotropic medium is considered, the indices i and j refer to a particular arbitrary choice of three orthogonal coordinate axes; the summation convention on repeated indices is used). The above relations are written with 〈 〉_{p} performed in the averagingvolume of central position x where the values $f\left(t\right)\widehat{\mathit{e}}\left(t\right)$ of the macroscopic pressure gradient were recorded. Recall that this volume served to construct a representative periodic or stationary random homogeneous medium, so that the values of the material parameters, porosity ϕ and kernels α_{ij} or k_{ij}, depend on x if the existing starting material has gradients of properties.
In harmonic regime the operator relations (14) and (15) become multiplications with the hat quantities becoming the usual frequencydependent functions, dynamic inertial/viscous tortuosity:
$$\mathrm{i\omega}{\rho}_{0}{\alpha}_{\mathrm{ij}}(\omega ,\mathit{x})\langle {\mathit{v}}_{j}{\rangle}_{p}={\nabla}_{i}P,$$(16)
and dynamic viscous/inertial permeability:
$$\varphi \left(\mathit{x}\right)\langle {\mathit{v}}_{i}{\rangle}_{p}=\frac{{k}_{\mathrm{ij}}\left(\omega ,\mathit{x}\right)}{\eta}{\nabla}_{j}P.$$(17)
At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem:
$$\mathrm{i\omega}{\rho}_{0}\mathit{v}=\nabla \pi +\eta \Delta \mathit{v}\nabla P,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(18)
$$\nabla \xb7\mathit{v}=0,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(19)
$$\pi :\mathrm{stationary}\mathrm{random}/\mathrm{periodic},\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(20)
with boundary conditions:
$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on}{S}_{p}^{\left(\mathrm{ub}\right)},$$(21)
and −∇P, a given spatial vector constant. From the three different unique solutions v^{(l)} associated with the source term $\nabla P=\frac{\delta P}{\delta \mathcal{l}}{\widehat{\mathit{e}}}^{\left(l\right)}$ oriented along the three different mutually orthogonal directions ${\widehat{\mathit{e}}}^{\left(l\right)}$ of the coordinate axes (${\widehat{\mathit{e}}}_{i}^{\left(l\right)}={\delta}_{\mathrm{il}}$, Kronecker symbol; δP variation of pressure along distance $\delta \mathcal{l}$ taken along ${\widehat{\mathit{e}}}^{\left(l\right)}$), we get, after averaging in the fluid and using the relations/definitions (16) and (17) the values of the desired functions k_{ij}(ω, x) or ${\alpha}_{\mathrm{ij}}^{1}(\omega ,\mathit{x})$:
$${k}_{\mathrm{ij}}(\omega ,\mathit{x})=\frac{\eta \varphi \left(\mathit{x}\right)}{\delta P/\delta \mathcal{l}}\langle {\mathit{v}}^{\left(j\right)}{\rangle}_{p}\xb7{\widehat{\mathit{e}}}^{\left(i\right)},\hspace{1em}{\alpha}_{\mathrm{ij}}^{1}(\omega ,\mathit{x})=\frac{\left(\mathrm{i\omega}{\rho}_{0}\right)}{\delta P/\delta \mathcal{l}}\langle {\mathit{v}}^{\left(j\right)}{\rangle}_{p}\xb7{\widehat{\mathit{e}}}^{\left(i\right)}.$$(22)
The obtained functions α_{ij}(ω, x) or k_{ij}(ω, x) are automatically symmetric [38, 39], expressing an OnsagerCasimir symmetry property, and they are related by:
$${k}_{\mathrm{ij}}\left(\omega ,\mathit{x}\right)=\frac{\nu \varphi \left(\mathit{x}\right)}{\mathrm{i\omega}}{\alpha}_{\mathrm{ij}}^{1}\left(\omega ,\mathit{x}\right).$$(23)
The frequencydependent dynamic permeability relation (17) between the macroscopic flow velocity 〈v〉 = ϕ〈v〉_{p} and the macroscopic pressure gradient, is referred to as a dynamic regime Darcy’s law. The other relation (16) is a dynamic regime Newton’s law.
To introduce the analogous notions of “dynamic thermal tortuosity and permeability”^{16}, we now proceed in a somewhat parallel manner. For heat exchange, we are interested in the description of the excesstemperaturepattern τ in air, appearing in the above finite coarsegrainingball of volume ${V}_{p}^{\left(b\right)}$ and pore surface ${S}_{p}^{\left(b\right)}$ centered at x. It obeys the following Fourier equation in ${V}_{p}^{\left(b\right)}$ (which is one among the complete set of governing coupled equations mentioned above (7.16) [9]),
$${\rho}_{0}{c}_{P}\frac{\mathrm{\partial}\tau}{\mathrm{\partial}t}={\beta}_{0}{T}_{0}\frac{\mathrm{\partial}p}{\mathrm{\partial}t}+\kappa \Delta \tau ,\hspace{1em}\mathrm{in}{V}_{p}^{\left(b\right)},$$(24)
with boundary conditions:
$$\tau =0,\hspace{1em}\mathrm{on}{S}_{p}^{\left(b\right)}.$$(25)
In the full set of equations of which (24) and (25) is an excerpt, p is a thermodynamic variable; it will now lose this status in the simplified way in which we will consider (24), isolated from the rest. Since in Biot theory we apply a divergencefree simplification of the representation of air motion at the pore scale, we must note at this time that we are considering only a situation where the pressure p within the averagingvolume can be viewed as the sum of a macroscopic part P (defined with ϕP = 〈p〉 or equivalently 〈pv〉 = P〈v〉), which is uniformly variable (P(t, r) = (r − x) · ∇P(t, x) + P(t, x), where ∇P(t, x) is a vector constant – it has null derivative ∂/∂r), and a deviatoric part p − P, with zero mean value (in either sense, 〈(p − P)〉 = 0 or 〈(p − P)v〉 = 0).
Parallel to what was said earlier for the velocity field, we now claim that, underlying Biot’s local theoretical approach is the notion that the excess temperature pattern τ(t, r) in the region of central position x and thus its macroscopic mean 〈τ〉(t, x), which is the quantity of interest to us, can be viewed as responding to the history of macroscopic terms β_{0}T_{0}∂P/∂t applied (in the present and in the past) at this very position. To define this response, we first expand, as before, in a suitable “uniform manner” (periodic or stationary random), the domains ${V}_{p}^{\left(b\right)}$ and ${S}_{p}^{\left(b\right)}$ so that they encompass the entire space outside the sphere, as if the material were macroscopically homogeneous. Next, we assume that the excess temperature pattern we seek, is approximately the same as that (periodic or stationary random) determined in the following wellposed problem:
$${\rho}_{0}{c}_{P}\frac{\mathrm{\partial}\tau}{\mathrm{\partial}t}={\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}+\kappa \Delta \tau ,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(26)
$$P=f\left(t\right),\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(27)
where f(t) is a spatial constant which represents the history of the macroscopic pressure at the central position x of the original averaging ball that served to construct the abstract homogeneous unbounded medium representation (${V}_{p}^{\left(\mathrm{ub}\right)}$ and ${S}_{p}^{\left(\mathrm{ub}\right)}$, periodic or stationary random), with boundary conditions:
$$\tau =0,\hspace{1em}\mathrm{on}{S}_{p}^{\left(\mathrm{ub}\right)}.$$(28)
We can justify the way this problem is written as follows. Since P represents now the time variations of the macroscopic pressure at the central position x of the original averaging sphere, the equation (26) is reminiscent of the Stokes equation (6) or (10) used to define the velocity pattern after incompressibility was introduced. It should be recalled that the π field did not exactly mimic the deviatoric overpressure, p − P, i.e. the porescale distribution of the physical overpressure p around its macroscopic average, but this was not important. The purpose of the plot was to provide an estimate, which could be considered precise in sufficiently simple geometries, of the relationship between the velocities and their mean, and the course of the macroscopic pressure gradients. Similarly here, by replacing the microscopic variable p in (26) with the macroscopic average at the original central sphere position, P, in (27), we do not attempt to mimic the physical overpressure distribution. The difference between the real overpressure p and the averaged P is however very small within Biot theory framework (see Appendix A) and moreover it will be of no consequence at all in estimating the relationship between the excess temperature τ (and thus its macroscopic mean) in the original sphere, and the course of the macroscopic term β_{0}T_{0}∂P/∂t with the spatial constant P as the averaged central variable^{17}. Finally, the problem (24) and (25) becomes wellposed if it is simplified and extended over the whole space in the way indicated to give (26)–(28).
As for the time variations of the uniform excitation P, they can be assumed to be arbitrary (as long as they are compatible with long wavelengths), but again, this does not mean that the intervening temporal dispersion will be general. For arbitrary time variations, the solution is a convolution of the elementary solution of the impulse Fourier problem considered in [18] (response to a Dirac delta heat source in time). This solution has the constraints given in Appendix C in [4], since it is constructed with purely damped normal modes (Avellaneda and Torquato’s viscous relaxationtimes analysis is easily transposed in a thermal relaxationtimes analysis, see [18]).
By writing the convolution we will learn that there are local macroscopic operators of “thermal tortuosity” $\widehat{\alpha}\mathrm{\prime}\left(\mathit{x}\right)$, associated with response, kernel factors α′(t, x), such that,
$${\rho}_{0}{c}_{P}\frac{\mathrm{\partial}\langle \tau {\rangle}_{p}}{\mathrm{\partial}t}\left(t\right)={\widehat{\alpha}}^{{\mathrm{\prime}}^{1}}\left(\mathit{x}\right){\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}\left(t\right)={\int}_{\mathrm{\infty}}^{t}\mathrm{}\mathrm{d}{t}^{\mathrm{\prime}}{\alpha}^{{\mathrm{\prime}}^{1}}\left(t{t}^{\mathrm{\prime}},\mathit{x}\right){\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}\left({t}^{\mathrm{\prime}}\right),$$(29)
or equivalently, local macroscopic operators of “thermal permeability” $\widehat{k}\mathrm{\prime}\left(\mathit{x}\right)$, associated with response, kernel factors k′(t, x), such that,
$$\varphi \left(\mathit{x}\right)\langle \tau {\rangle}_{p}\left(t\right)=\frac{{\widehat{k}}^{\mathrm{\prime}}\left(\mathit{x}\right)}{\kappa}{\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}\left(t\right)=\frac{1}{\kappa}{\int}_{\mathrm{\infty}}^{t}\mathrm{}\mathrm{d}t\mathrm{\prime}k\mathrm{\prime}(tt\mathrm{\prime},\mathit{x}){\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}\left(t\mathrm{\prime}\right).$$(30)
The above relations are written with 〈 〉_{p} performed in the averagingvolume of the central position x where the f(t) values of the macroscopic source term β_{0}T_{0}∂P/∂t have been recorded. Recall that the values of the material parameters, porosity ϕ and kernels α′ and k′, depend on x if there are gradients in macroscopic properties; note that for shortness this argument x is omitted in the fields in the equations ((14)–(15), (29)–(30)).
In harmonic regime the operator relations (29) and (30) become multiplications with the hat quantities becoming the usual frequencydependent functions, dynamic thermal tortuosity α′(ω, x):
$$\mathrm{i\omega}{\rho}_{0}{c}_{P}\alpha \mathrm{\prime}(\omega ,\mathit{x})\langle \tau {\rangle}_{p}={\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t},$$(31)
and dynamic thermal permeability k′(ω, x):
$$\varphi \left(\mathit{x}\right)\langle \tau {\rangle}_{p}=\frac{{k}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)}{\kappa}{\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}.$$(32)
At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem:
$$\mathrm{i\omega}{\rho}_{0}{c}_{P}\tau =\mathrm{i\omega}{\beta}_{0}{T}_{0}P+\kappa \Delta \tau ,\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(33)
$$P=\mathrm{a}\mathrm{given}\mathrm{spatial}\mathrm{constant},\hspace{1em}\mathrm{in}{V}_{p}^{\left(\mathrm{ub}\right)},$$(34)
with boundary conditions:
$$\tau =0,\hspace{1em}\mathrm{on}{S}_{p}^{\left(\mathrm{ub}\right)}.$$(35)
From the unique solution τ of this diffusioncontrolled Fourier problem, that will be automatically stationary random or periodic, will derive, by averaging in the fluid and by using the relations/definitions (31) and (32) the values of the wanted functions α′(ω, x) or k′(ω, x):
$$k\mathrm{\prime}(\omega ,\mathit{x})=\frac{\kappa \varphi \left(\mathit{x}\right)}{{\beta}_{0}{T}_{0}\mathrm{\partial}P/\mathrm{\partial}t}\langle \tau {\rangle}_{p},{\alpha}^{\mathrm{\prime}1}(\omega ,\mathit{x})=\frac{\left(\mathrm{i\omega}{\rho}_{0}{c}_{P}\right)}{{\beta}_{0}{T}_{0}\mathrm{\partial}P/\mathrm{\partial}t}\langle \tau {\rangle}_{p},$$(36)
which are related by:
$${k}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)=\frac{{\nu}^{\mathrm{\prime}}\varphi \left(\mathit{x}\right)}{\mathrm{i\omega}}{\alpha}^{\mathrm{\prime}1}\left(\omega ,\mathit{x}\right).$$(37)
The frequencydependent dynamic thermal permeability relation (32) between the macroscopic excess temperature 〈τ〉 = ϕ〈τ〉_{p} and the macroscopic source term, is referred to as a dynamic regime “thermal Darcy’s law”. The other relation (31) may be referred to as a dynamic regime “thermal Newton’s law”. It remains here to explain how these notions of thermal dynamic permeability and tortuosity determine the effective bulk modulus of the air. The link between thermal dynamic permeability/tortuosity and effective bulk modulus has been presented e.g. in [4, 9, 18]. For completeness we briefly recall here the relation.
As derived in [8, 9], the thermodynamic equation of state of the saturating fluid leads, after linearization, to the following relation between excess pressure, condensation b ≡ ρ′/ρ_{0} where ρ′ is the excess density, and excess temperature:
$$\gamma {\chi}_{0}p=b+{\beta}_{0}\tau ,\hspace{1em}\mathrm{in}{V}_{p}^{\left(b\right)},$$(38)
where the notation χ_{0} is employed for the adiabatic compressibility 1/K_{a} (setting, τ = 0, in (38), the resulting relation, γχ_{0}p = b, must describe the isothermal excess pressurecondensation relation, so that, obviously, γχ_{0} = 1/K_{0}). Taking the average and the time derivative it gives the relation, since 〈p〉_{p} = P,
$$\gamma {\chi}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t}=\frac{\mathrm{\partial}\langle b{\rangle}_{p}}{\mathrm{\partial}t}+{\beta}_{0}\frac{\mathrm{\partial}\langle \tau {\rangle}_{p}}{\mathrm{\partial}t}.$$(39)
In harmonic regime e^{−iωt} and working with the complex amplitudes we can use (31) to express the righthand side in (39). Using the general thermodynamic identity seen in Section 2 and suppressing the common factor (−iω) this gives after straightforward calculation the relation,
$${\chi}_{0}\left[\gamma \frac{\gamma 1}{{\alpha}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)}\right]P=\langle b{\rangle}_{p}.$$(40)
Therefore in harmonic regime we find a relation having the form,
$$P={K}_{f}(\omega ,\mathit{x})\langle b{\rangle}_{p}={\chi}_{f}^{1}(\omega ,\mathit{x})\langle b{\rangle}_{p},$$(41)
with
$${K}_{f}^{1}\left(\omega ,\mathit{x}\right)={\chi}_{f}\left(\omega ,\mathit{x}\right)={\chi}_{0}\left[\gamma \frac{\gamma 1}{{\alpha}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)}\right],$$(42)
$$={\chi}_{0}\left[\gamma \left(\gamma 1\right)\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}\varphi}{k}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)\right].$$(43)
It gives the Fourier coefficients of the kernel functions ${\chi}_{f}^{1}(t,\mathit{x})$ or K_{f}(t, x) which appear in the arbitrary timevariableregime, operator or integral equation,
$$P(t,\mathit{x})={\widehat{\chi}}_{f}^{1}\langle b{\rangle}_{p}(t,\mathit{x})={\int}_{\mathrm{\infty}}^{t}\mathrm{}{\chi}_{f}^{1}(tt\mathrm{\prime},\mathit{x})\langle b{\rangle}_{p}(t\mathrm{\prime},\mathit{x})\mathrm{d}t\mathrm{\prime},$$(44)
$$={\widehat{K}}_{f}\langle b{\rangle}_{p}(t,\mathit{x})={\int}_{\mathrm{\infty}}^{t}\mathrm{}{K}_{f}(tt\mathrm{\prime},\mathit{x})\langle b{\rangle}_{p}(t\mathrm{\prime},\mathit{x})\mathrm{d}t\mathrm{\prime}.$$(45)
We recognize here, the relation announced in Section 1. The function χ_{f}(ω, x), resp. its inverse K_{f}(ω, x), represents a local dynamic compressibility, resp. dynamic bulk modulus, of the saturating fluid, that depends on frequency because of the thermal exchanges between fluid and solid.
Finally, we note that it is customary to define a normalized dynamic compressibility, β(ω, x), such that, in harmonic regime,
$${\chi}_{0}\beta \left(\omega ,\mathit{x}\right)\varphi \left(\mathit{x}\right)\frac{\mathrm{\partial}P}{\mathrm{\partial}t}={\chi}_{0}\beta \left(\omega ,\mathit{x}\right)\frac{\mathrm{\partial}\langle p\rangle}{\mathrm{\partial}t}=\frac{\mathrm{\partial}\langle b\rangle}{\mathrm{\partial}t}=\nabla \xb7\langle \mathit{v}\rangle .$$(46)
It is given by,
$$\beta \left(\omega ,\mathit{x}\right)=\frac{{\chi}_{f}\left(\omega ,\mathit{x}\right)}{{\chi}_{0}}=\gamma \frac{\gamma 1}{{\alpha}^{\mathrm{\prime}}\left(\omega ,\mathit{x}\right)}.$$(47)
Note that in (46) we have directly written that, ∂〈b〉/∂t = −∇ · 〈v〉, using the identity 〈∇ · v〉 = ∇ · 〈v〉, which is valid (with no indice p) even if the porosity varies spatially. We refer the reader to [9] where the explanations are given on how this identity is obtained (where in the left, ∇ = ∂/∂r, and in the right, ∇ = ∂/∂x).
4 Generalization of the method of Pontrjagin, Andronow and Witt, in the presence of radioactive decay of diffusing species
As will be seen in Section 5, there is a direct proportionality between the distribution of excess temperature field τ(r) in thermal problem defined in (33)–(35) and the mean survival time $\stackrel{\u0305}{t}\left(\mathit{r}\right)$ of a tiny particle released at position r, moving and decaying under the combined action of thermal agitation (Brownian motion) and radioactive disintegration, and instantly absorbed by reaching a specified boundary (that of the pore walls). (The overbar designates the average over realizations of the random walk in the limit of large numbers). In this section, we will first study how to compute the mean survival time $\stackrel{\u0305}{t}\left(\mathit{r}\right)$ of the particle with instant absorption (instant decay), at the pore walls. In a paper on the statistical treatment of dynamic systems, Pontrjagin et al. [17], have discussed a general method, allowing in particular to establish a differential equation for the average time it takes a Brownian motion particle moving under the combined effect of thermal agitation and a stationary field of force, to reach a defined boundary. Focusing on this particular problem, an helpful detailed account of Pontrjagin, Andronow and Witt’s general method has been presented by Lifson and Jackson in [16]. Here, closely following their exposition, we will establish a differential equation for the average time it takes for a Brownian motion particle to reach the pore walls, moving under the combined effects of thermal agitation and superimposed radioactive decay of the diffusing species. The method of calculation of $\stackrel{\u0305}{t}\left(\mathit{r}\right)$ in a complex geometry will be based on the known solutions of the differential equation in the most simple open (disk and ball in 2D and 3D, see Eqs. (73) and (74)).
In their presentation, Lifson and Jakson [16], as in the general considerations in [17], started with a closed region V bounded by a surface S. Let us substitute for volume V our pore space V_{p}, occupied by air, and for enclosing surface S, the pore wall surface S_{p} between air and porous solid. Our pore space V_{p} is not a closed region but this is not a significant difference. What is important is that this space is only bounded by the surface S_{p}. If it happens in our following reasonings that V_{p} and S_{p} are not macroscopically homogeneous regions, we will understand to identify them with the homogeneous domains we referred previously by the notations ${V}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, ${S}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, with the dependence (x) not indicated, for simplicity. Then, with the substitution of ${V}_{p}^{\left(\mathrm{ub}\right)}$, ${S}_{p}^{\left(\mathrm{ub}\right)}$, for V, S, in [16], only minimal change in the considerations will occur. This is why we take the liberty, in the following, to paraphrase almost completely parts of Lifson and Jakson’s exposition, mutatis mutandis only what is necessary in the presence of the superimposed radioactive decay of the diffusing particles.
As we will see, varying the radioactive decay time τ_{b}, from ∞ (no disintegration, low frequency limit) to 0 (instant disintegration, high frequency limit), will allow us to explore the effective bulkmodulus function in its LaplaceTransform variable representation, ${K}_{f}^{\mathcal{L}}\left(s\right)$, from s = 0 to s = ∞, where s is related to ω according to s = −iω using the convention e^{−iωt}. In [16], the conservative force field, F(r), acts on the diffusing species and is required to be finite everywhere in the V_{p} region. For the sake of generality in this section we retain it, though we need not consider it in the sequel. We will later (Sect. 5) set, F(r) = 0.
We consider Brownian motion particles, moving, in the region V_{p}, under the combined influence of thermal agitation and the force F; the new element compared to [16] is that each diffusing particle undergoes “radioactive decay” with time characteristic τ_{b} in the volume V_{p}, i.e., if present (in bulk) at time t the particle has probability κ_{b} dt = dt/τ_{b} of being destroyed in the infinitesimal time interval (t, t +dt).
“Disintegration” may refer to the particle itself that can be removed e.g. by chemical reactions, or to some physical excitation carried by the particle that can be deexcited by interaction in the bulk. In all cases we speak of “the particle” and of its disappearance. We are not concerned here with the physical implementation but only with the posed mathematics of the diffusiondisintegration process which is assumed to mimick Brownian motion and radioactive decay laws. We define the cumulative probability function $\mathcal{D}(\mathit{r},t)$ (replacing the Lifson and Jackson cumulative probability function $\mathcal{W}(\mathit{r},t)$), as the cumulative probability function that a particle located at point r at time t = 0 will reach the bounding surface S_{p} during the time interval (0, t), or else, will disappear in the bulk V_{p} without reaching S_{p} during the same interval.
The “$\mathcal{D}$” – for “$\mathcal{D}$isappearance” (or “$\mathcal{D}$ecay”, “$\mathcal{D}$isintegration”, “$\mathcal{D}$elete” or $\mathcal{D}$eath”) – qualification for this probability is appropriate if any particle touching the surface S_{p}, or the excitation it carries, is instantly removed; $\mathcal{D}(\mathit{r},t)$ only differs from $\mathcal{W}(\mathit{r},t)$ in [16] in that the superimposed radioactive decay (of diffusing particles or their excitation), provides an additional spontaneous means of removal (in bulk). A number of simple properties previously listed for $\mathcal{W}(\mathit{r},t)$ in [16], still hold for $\mathcal{D}(\mathit{r},t)$. Evidently, for a point s on the surface S_{p}, $\mathcal{D}(\mathit{s},t)=1$ for all t. For an interior point we have $\mathcal{D}(\mathit{r},\mathrm{\infty})=1$, while $\mathcal{D}(\mathit{r},0)=0$. A difference with [16] is that the condition,
$$\underset{\tau \to 0}{\mathrm{lim}}\left[\mathcal{W}\left(\mathit{r},\tau \right)/\tau \right]=0,$$(48)
now becomes, after accounting for the finite rate of radioactive decay,
$$\underset{\tau \to 0}{\mathrm{lim}}\left[\mathcal{D}\left(\mathit{r},\tau \right)/\tau \right]=\frac{1}{{\tau}_{b}}={\kappa}_{b}.$$(49)
To see this, recall that $\mathcal{D}(\mathit{r},\tau )$ is the probability that a particle (or its excitation), starting from r at t = 0, dies before t = τ, either by touching S_{p}, or by spontaneously decaying in the bulk. So to calculate $\mathcal{D}$, we launch N → ∞ particles at r at t = 0, and, after a time interval $\tau $, we see N″ paths of particles that have not touched S_{p} or died in the bulk. The remaining N′ = N − N″ paths, are dead paths, either stopped in the bulk when particles died spontaneously, or stopped at the boundary S_{p} where instant decay occurs. Then, $\mathcal{D}(\mathit{r},\tau )={\mathrm{lim}}_{N\to \mathrm{\infty}}(N\mathrm{\prime}/N)$. When τ → 0 there is not enough time to reach S_{p}, so the nondead path amounts to N″ predicted on the only basis of mass “radioactive” absorption: N″ = N(1 − κ_{b}τ); consequently N′ = Nκ_{b}τ and $\mathcal{D}={\kappa}_{b}\tau $, as stated in (49).
We can now relate (as Lifson and Jackson did for probability $\mathcal{W}$) the probability $\mathcal{D}$ at time τ, to its value at the next time t + τ. For this purpose we define the function ι(r, t, ρ), such that, ι(r, t, ρ)dρ, is the transition probability for a particle, to go from point r, in V_{p}, at time t = 0, to a volume element dρ about another point, ρ, in V_{p}, at time t, without having touched the surface S_{p} during the time interval (0, t). This is not the same as υ(r, t, ρ)dρ in [16], because, as our particle now decays in a manner typical of radioactive decay, ι is therefore reduced compared to its counterpart υ without radioactive decay:
$$\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)={e}^{\tau /{\tau}_{b}}\upsilon \left(\mathit{r},\tau ,\mathit{\rho}\right).$$(50)
We may call $\mathcal{L}(\mathit{r},t)\equiv {\int}_{{V}_{p}}\mathrm{}\iota (\mathit{r},t,\mathit{\rho})\mathrm{d}\mathit{\rho}$ the “$\mathcal{L}$ife expectancy” from going to r, at t = 0, to anywhere else in the pore space, at time t. As a particle must be counted as either deceased or living, the function ι(r, t, ρ) is normalized in the sense that,
$$\mathcal{L}\left(\mathit{r},t\right)\equiv {\int}_{{V}_{p}}\mathrm{}\iota \left(\mathit{r},t,\mathit{\rho}\right)d\mathit{\rho}=1\mathcal{D}\left(\mathit{r},t\right).$$(51)
From (51) and (49) we also note that,
$$\underset{\tau \to 0}{\mathrm{lim}}{\int}_{{V}_{p}}\mathrm{}\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)d\mathit{\rho}=1{\kappa}_{b}\tau +\dots $$(52)
It is now clear that the following integral equation holds:
$$\mathcal{D}\left(\mathit{r},t+\tau \right)=\mathcal{D}\left(\mathit{r},\tau \right)+{\int}_{{V}_{p}}\mathrm{}\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)\mathcal{D}\left(\mathit{\rho},t\right)\mathrm{d}\mathit{\rho}.$$(53)
The meaning of the above equation is as follows: If a particle has reached the surface S_{p} during the time interval (0, t + τ) or has died in the volume V_{p} during this time interval without touching it, then, either it has reached S_{p} or died in V_{p} during the time interval (0, τ), or else, it gets to the point ρ at time τ without touching the surface S_{p} or decaying in V_{p}, and then, reaches from this point the surface S_{p} or decays in V_{p} without touching it, in the time interval (τ, τ + t)^{18}.
We now recall known results for the successive moments of the ordinary transition probability υ_{0}(r, t, ρ) in an infinite medium, which will immediately translate into identical results for the successive moments of ι(r, t, ρ). For Brownian motion without disintegration in an infinite medium in the presence of a field of force, the particle has an average velocity at a point, proportional to the force, i.e., the following limit exists:
$$\underset{\tau \to 0}{\mathrm{lim}}\frac{\langle {\rho}_{i}{r}_{i}\rangle}{\tau}=\frac{{F}_{i}}{f}=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}\left({\rho}_{i}{r}_{i}\right){\upsilon}_{0}\left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho},$$(54)
where f is the hydrodynamic friction coefficient and υ_{0}(r, τ, ρ) the transition probability in the infinite medium. Furthermore, the diffusion coefficient, D = kT/f, is related to the second moment of the displacement, according to,
$$\underset{\tau \to 0}{\mathrm{lim}}\frac{\langle \left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right)\rangle}{\tau}=2D{\delta}_{\mathrm{ij}},$$(55)
$$=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}\left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right){\upsilon}_{0}\left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho}.$$(56)
To see this we note that the effect of the, forcedriven, constant drift velocity F_{i}/f, which brings in (ρ_{i} − r_{i}) a term τF_{i}/f which superposes to a random, nondriven one, will bring a cross term, τ^{−1}(F_{i}F_{j})τ^{2}/f^{2}, which does not contribute in the limit τ → 0, and two noncrossterms, that disappear by the absence of direction of the random nondriven motion. Thus, no forcedriven effect appears in the second moment, and the transition probability υ_{0}(r, t, ρ) can be replaced by the one without force, Gaussian in each direction, $={e}^{\mathbf{\rho}\mathbf{r}{}^{2}/2\mathrm{dDt}}/\left(2\pi \mathrm{dDt}\right)$, where d is the dimension. The result then expresses Einstein’s translational Brownian motion formula, $\stackrel{\u0305}{\left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right)}=2\mathrm{D\tau}{\delta}_{\mathrm{ij}}$. For all higher moments it can be checked that,
$$\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}\langle {\left({\rho}_{i}{r}_{i}\right)}^{l}{\left({\rho}_{j}{r}_{j}\right)}^{m}{\left({\rho}_{k}{r}_{k}\right)}^{n}\rangle =0,$$(57)
$$=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}{\left({\rho}_{i}{r}_{i}\right)}^{l}{\left({\rho}_{j}{r}_{j}\right)}^{m}{\left({\rho}_{k}{r}_{k}\right)}^{n}{\upsilon}_{0}\left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho},$$(58)
Now, in view of the fact that the function ι(r, t, ρ) must approach the ordinary (i.e. infinite medium) transition probability as τ approaches to zero, the moments of ι(r, t, ρ) must approach the moments of the ordinary transition probability. We then have, from equations (54), (55), and (57),
$$\underset{\tau \to 0}{\mathrm{lim}}\frac{\langle {\rho}_{i}{r}_{i}\rangle}{\tau}=\frac{{F}_{i}}{f}=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}\left({\rho}_{i}{r}_{i}\right)\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho},$$(60)
$$2D{\delta}_{\mathrm{ij}}=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}\left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right)\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho},$$(61)
$$0=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}{\int}_{{V}_{p}}\mathrm{}{\left({\rho}_{i}{r}_{i}\right)}^{l}{\left({\rho}_{j}{r}_{j}\right)}^{m}{\left({\rho}_{k}{r}_{k}\right)}^{n}\iota \left(\mathit{r},\tau ,\mathit{\rho}\right)\mathrm{d}\mathit{\rho},\left(l+m+n\right)>2.$$(62)
We now use the integral equation for $\mathcal{D}(\mathit{r},t)$ (Eq. (53)), to derive a partial differential equation. Expanding $\mathcal{D}(\mathit{\rho},t)$ about the point r we have
$$\frac{\mathrm{\partial}\mathcal{D}\left(\mathit{r},t\right)}{\mathrm{\partial}t}=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}\left[\mathcal{D}\left(\mathit{r},t+\tau \right)\mathcal{D}\left(\mathit{r},t\right)\right],$$(63)
$$=\underset{\tau \to 0}{\mathrm{lim}}{\tau}^{1}\{\mathcal{D}(\mathit{r},t)+\mathcal{D}(\mathit{r},\tau )+{\int}_{{V}_{p}}\mathrm{}d\mathit{\rho}\iota (\mathit{r},\tau ,\mathit{\rho})[\mathcal{D}(\mathit{r},t)+\sum _{i}\mathrm{}({\rho}_{i}{r}_{i})\frac{\mathrm{\partial}\mathcal{D}}{\mathrm{\partial}{r}_{i}}+\frac{1}{2!}\sum _{i,j}\mathrm{}\left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right)\frac{{\mathrm{\partial}}^{2}\mathcal{D}}{\mathrm{\partial}{r}_{i}\mathrm{\partial}{r}_{j}}+\frac{1}{3!}\sum _{i,j,k}\mathrm{}\left({\rho}_{i}{r}_{i}\right)\left({\rho}_{j}{r}_{j}\right)\times ({\rho}_{k}{r}_{k})\frac{{\mathrm{\partial}}^{3}\mathcal{D}}{\mathrm{\partial}{r}_{i}\mathrm{\partial}{r}_{j}\mathrm{\partial}{r}_{k}}+\dots \left]\right\}.$$(64)
Utilizing the various limiting properties of ι and $\mathcal{D}$ presented in (49), (52), (60), (61), (62), we get a differential equation for the cumulative probability $\mathcal{D}(\mathit{r},t)$ that a particle launched at r in V_{p} at t = 0, vanishes before t, either by touching S_{p} or decaying in the bulk:
$$\frac{\mathrm{\partial}\mathcal{D}}{\mathrm{\partial}t}={f}^{1}\mathit{F}\xb7\nabla \mathcal{D}+D\Delta \mathcal{D}{\kappa}_{b}\mathcal{D}+{\kappa}_{b}.$$(65)
This interprets as the result obtained for $\mathcal{W}$ in Lifson and Jackson [16] equation (9), modified by the presence of a term ${\kappa}_{b}\mathcal{L}$ added in the lefthand side to account for an additional rate of change coming from the radioactive disintegration.
Now observe that the equation, $\mathcal{D}(\mathit{r},t+\mathrm{d}t)=\mathcal{D}(\mathit{r},t)+\mathrm{d}t\mathrm{\partial}\mathcal{D}(\mathit{r},t)/\mathrm{\partial}t,$ should be interpreted as follows: If a particle starting from r at t = 0 vanishes during the time interval (0, t + dt) (either by reaching S_{p} or decaying in the bulk V_{p}), then, either it vanishes during the time interval (0, t) (by reaching S_{p} or decaying in the bulk V_{p}), or it is still alive at t, and then, during the time interval (t, t + dt), reaches for the first time S_{p} or vanishes in the bulk. This shows that $\mathrm{d}t\mathrm{\partial}\mathcal{D}(\mathit{r},t)/\mathrm{\partial}t$ is the probability of a particle to vanish (reach for the first time S_{p} or vanish in the bulk) during (t, t + dt), when starting from r at time t = 0. Therefore, the average survival time for a particle is:
$$\stackrel{\u0305}{t}\left(\mathit{r}\right)={\int}_{0}^{\mathrm{\infty}}\mathrm{}t\frac{\mathrm{\partial}\mathcal{D}}{\mathrm{\partial}t}\left(\mathit{r},t\right)\mathrm{d}t.$$(66)
From this we finally obtain the following desired differential equation satisfied by $\stackrel{\u0305}{t}\left(\mathit{r}\right)$:
$$\left(1/f\right)\mathit{F}\xb7\nabla \stackrel{\u0305}{t}+D\Delta \stackrel{\u0305}{t}{\kappa}_{b}\stackrel{\u0305}{t}=1,\hspace{1em}\mathrm{in}{V}_{p},$$(67)
with boundary condition:
$$\stackrel{\u0305}{t}=0,\hspace{1em}\mathrm{on}{S}_{p}.$$(68)
When the field F is zero one may directly write the solutions for a line, cylinder, and sphere. One has, with notation ${\mu}_{a}\equiv a\sqrt{{\kappa}_{b}/D}$,
One dimensional case: length of line 2L, x ∈ [−L, L],
$$\stackrel{\u0305}{t}\left(x\right)={\tau}_{b}\left[1\frac{\mathrm{cosh}\left({\mu}_{x}\right)}{\mathrm{cosh}\left({\mu}_{L}\right)}\right];$$(69)
Twodimensional case: radius of cylinder R, r ∈ [0, R],
$$\stackrel{\u0305}{t}\left(r\right)={\tau}_{b}\left[1\frac{{I}_{0}\left({\mu}_{r}\right)}{{I}_{0}\left({\mu}_{R}\right)}\right];$$(70)
Threedimensional case: radius of sphere R, r ∈ [0, R],
$$\stackrel{\u0305}{t}\left(r\right)={\tau}_{b}\left[1\frac{{\mu}_{R}\mathrm{sinh}({\mu}_{r})}{{\mu}_{r}\mathrm{sinh}\left({\mu}_{R}\right)}\right].$$(71)
5 Efficient algorithm to compute ${K}_{f}\left(\omega \right)$ and other functions and intervening parameters
With F = 0 we can directly compare the average survival time problem (67)–(68) with the excess temperature pattern problem (33)–(34). Examination of both shows that the following quantities, with the same dimensions, can be put into correspondence:
$$\stackrel{\u0305}{t}\leftrightarrow \frac{1}{\mathrm{i\omega}\frac{{\beta}_{0}{T}_{0}}{{\rho}_{0}{c}_{P}}P}\tau ;$$(72)
$$D\leftrightarrow \frac{\kappa}{{\rho}_{0}{c}_{P}}\equiv {\nu}^{\text{'}};$$(73)
$${\kappa}_{b}\equiv \frac{1}{{\tau}_{b}}\leftrightarrow \mathrm{i\omega}.$$(74)
Therefore, when we have solved the mean survival time problem, we have solved the thermal problem. For example we can compute the thermal dynamic permeability by, first, writing the relation (e.g. combine (32) with the correspondence (72)–(74)):
$${k}^{\mathrm{\prime}}\left(\omega \right)=\varphi D\langle \stackrel{\u0305}{t}\rangle ,$$(75)
and next, substituting in the calculations that are made to evaluate the product $D\langle \stackrel{\u0305}{t}\rangle $, the values D = ν′ and κ_{b} = −iω of the coefficients D and κ_{b}.
Likewise, we can compute the bulk modulus K_{f}(ω) by, first, writing the relation (e.g. combine (43) with (75) and the correspondence (72)–(74)):
$${K}_{f}^{1}\left(\omega \right)={K}_{a}^{1}\left[\gamma \left(\gamma 1\right){\kappa}_{b}\langle \stackrel{\u0305}{t}\rangle \right],$$(76)
and next, substituting the above values of D and κ_{b} in the expressions obtained for the product ${\kappa}_{b}\langle \stackrel{\u0305}{t}\rangle $. In fact, it means that in the Laplacevariable representation of the kernel K_{f}(t) of the bulk modulus operator ${\widehat{K}}_{f}$, (45)^{19}, we simply have, ${\left({K}_{f}^{\mathcal{L}}\right)}^{1}\left(s\right)={K}_{a}^{1}\left[\gamma (\gamma 1){\kappa}_{b}\langle \stackrel{\u0305}{t}\rangle \right]$, with this time, D = ν′ and κ_{b} = s: The problem of mean survival time, most naturally translates in the description of thermal response factors expressed in Laplace representation.
Having established the connection between the mean survival time problem and the thermal problem, let us now focus on the former and show, how a simple geometric construction of the Brownian random motion path considered in Torquato and Kim [13], see Figure 1, leads us to a simple geometric algorithm to calculate the mean survival time $\stackrel{\u0305}{t}$ and its average $\langle \stackrel{\u0305}{t}\rangle $ over position or over configuration (see in Appendix A, the respective viewpoints of Lorentz and Gibbs).
Figure 1 Schematic representation of the construction of Brownian motion paths r_{1}, r_{2}, …, r_{n} with succession of characteristic radii [R_{1}, R_{2}, …, R_{n}]. 
For the sake of simplicity of Figure 1 and the discussion here we take our material to be an airsaturated array of parallel fibers with identical radii. Thus only the simulation of diffusion reactions in a medium with static identical circular cylindrical traps in 2D needs to be performed; the generalization to arbitrary microgeometries and 3D is obvious and poses no difficulty in principle. The algorithm to calculate the mean survival time $\stackrel{\u0305}{t}\left(\mathit{r}\right)$ of random walkers released in one configuration at a given place, or to calculate its average $\langle \stackrel{\u0305}{t}\rangle $ either viewed as a volume average (in a given configuration) when walkers are released at any place with equal probability (Lorentz viewpoint), or as an ensemble average when they are released at the same position but the configuration is varied appropriately (Gibbs viewpoint), is now described.
Let r_{1} be one arbitrary initial position of the random walker. As observed by Zheng and Chiew [40] and Torquato and Kim [13] the zigzag motion of the random walker need not be simulated step by step. Instead, one constructs the largest concentric circle/sphere of radius R_{1} which does not overlap any solid trap. As the random motion establishes no preference on the directions, the next position r_{1} of the walker is chosen at random on this circle/sphere. In this way, Brownian motion paths, r_{1}, r_{2}, …, r_{n}, are generated (see Fig. 1), which can be stopped when r_{n} is within ϵ of a solid trap (R_{n} < ϵ), with ϵ taken small enough to ensure that in the problem at hand, the truncation results in insignificant errors.
We now make the necessary modifications to the scheme utilized by the above authors to account for the “radioactive decay” of our random walkers. (These authors utilized the Brownian motion path construction to compute the trapping constant, which is the inverse of the static thermal permeability, whereas we will use it to compute the dynamic thermal permeability or even more directly, the associated relaxation function, Eq. (90)).
Consider diffusion in one of the circular 2D regions of radius R_{i} defined in the above construction. With R the radius of one such region, we define a hitting probability p(R) that a random walker initially at the center, survives until reaching the boundary. We will start by evaluating this probability, which is unity in the absence of “radioactive decay” (the diffusing walker goes to infinity in an unbounded medium), and smaller than 1 if it exists (not all walkers survive to reach the boundary).
For a random walker starting at the origin O of the region, there are two types of outcomes to the random walk with decay, see Figure 2: (i) the decay does not occur before reaching the boundary of radius R, the associated probability is that we want to determine, p(R), and we know that the corresponding mean time (${\stackrel{\u0305}{t}}_{R}^{\left(s\right)}$, s for surface), to reach R, will be given by the inverse of Einstein’s relation, i.e.^{20},
$${\stackrel{\u0305}{t}}_{R}^{\left(s\right)}=\frac{{R}^{2}}{4\mathrm{D}},\left(2\mathrm{D}\right);{\stackrel{\u0305}{t}}_{R}^{\left(s\right)}=\frac{{R}^{2}}{6\mathrm{D}},\left(3\mathrm{D}\right),$$(77)
(ii) the decay occurs before radius R is reached, the associated probability is 1 − p(R) and we assume that these aborted random walks are performed in an average time ${\stackrel{\u0305}{t}}_{R}^{\left(b\right)}$ (b for bulk). Note that ${\stackrel{\u0305}{t}}_{\mathrm{\infty}}^{\left(b\right)}$ is the average survival time τ_{b} of the walker in the infinite medium.
Figure 2 Possible outcomes of the random walk with disintegration, associated probabilities and survival times. 
Let ${\stackrel{\u0305}{t}}_{R}$ be the mean survival time of the random walker starting at the origin, when the latter is instantly absorbed at the radius R, and subject to the radioactive decay in the bulk. We had it determined in Section 4 by equation (70) (at r = 0) in 2D and by equation (71) (at r = 0) in 3D:
$${\stackrel{\u0305}{t}}_{R}={\tau}_{b}\left[1\frac{1}{{I}_{0}\left({\mu}_{R}\right)}\right],\left(2\mathrm{D}\right);\hspace{1em}{\stackrel{\u0305}{t}}_{R}={\tau}_{b}\left[1\frac{{\mu}_{R}}{\mathrm{sinh}\left({\mu}_{R}\right)}\right],\left(3\mathrm{D}\right).$$(78)
The above two types of outcomes allow us to write it, also, in the following form:
$${\stackrel{\u0305}{t}}_{R}=\left[1p\left(R\right)\right]{\stackrel{\u0305}{t}}_{R}^{\left(b\right)}+p\left(R\right){\stackrel{\u0305}{t}}_{R}^{\left(s\right)}.$$(79)
To get another relationship between the same two, but yet unknown quantities, p(R) and ${\stackrel{\u0305}{t}}_{R}^{\left(b\right)}$, let us now imagine that the boundary at radius R is fictitious and drawn in an infinite medium. Let us again follow the particles released in the center O. We know that on average they live the time τ_{b}. A proportion 1 − p(R) of them do not reach radius R and live the average time ${\stackrel{\u0305}{t}}_{R}^{\left(b\right)}$. The remaining proportion p(R) reach R and then will later disappear at some instant. The latter particles therefore live first the average time ${\stackrel{\u0305}{t}}_{R}^{\left(s\right)}$ to reach R, and then take the additional averagetime τ_{b} to disappear. Thus, we have the obvious relation,
$${\tau}_{b}=\left[1p\left(R\right)\right]{\stackrel{\u0305}{t}}_{R}^{\left(b\right)}+p\left(R\right)\left[{\stackrel{\u0305}{t}}_{R}^{\left(s\right)}+{\tau}_{b}\right].$$(80)
Combined with the preceding, it reads, ${\tau}_{b}={\stackrel{\u0305}{t}}_{R}+p\left(R\right){\tau}_{b}$ (irrespective of the dimension 2 or 3), that is:
$${\stackrel{\u0305}{t}}_{R}={\tau}_{b}\left[1p\left(R\right)\right].$$(81)
Inserting in (81) the mean survival time expressions (78) we therefore find for the hitting probability p(R) the following expressions in the 2D and 3D cases:
$$p\left(R\right)=\frac{1}{{I}_{0}\left({\mu}_{R}\right)},\left(2\mathrm{D}\right);\hspace{1em}p\left(R\right)=\frac{{\mu}_{R}}{\mathrm{sinh}\left({\mu}_{R}\right)},\left(3\mathrm{D}\right);\hspace{1em}{\mu}_{R}\equiv R\sqrt{\frac{{\kappa}_{b}}{D}}.$$(82)
Having determined the probability p(R) (for having survived till crossing for the first time the radius r = R, after having started at the center r = 0) and on the other hand the mean time ${\stackrel{\u0305}{t}}_{R}$ (for crossing for the first time the radius r = R or decaying before touching it, after having started at r = 0) we now can write the following simple expression:
$$\begin{array}{c}\stackrel{\u0305}{t}\left({\mathit{r}}_{1},{\mathit{r}}_{2},\dots ,{\mathit{r}}_{n}\right)={\stackrel{\u0305}{t}}_{{R}_{1}}+p\left({R}_{1}\right){\stackrel{\u0305}{t}}_{{R}_{2}}+p\left({R}_{1}\right)p\left({R}_{2}\right){\stackrel{\u0305}{t}}_{{R}_{3}}+\dots \\ \dots +p\left({R}_{1}\right)\dots p\left({R}_{n1}\right){\stackrel{\u0305}{t}}_{{R}_{n}},\end{array}$$(83)
for the mean time for trapping, associated with the path construction, r_{1}, r_{2}, …, r_{n}, with R_{i} = r_{i+1} − r_{i}. Trapping means here, spontaneously decaying in bulk, or else, instantly decaying by reaching the solid or representative end point of the construction, after having started at r_{1} at t = 0. Finally, by inserting in (83) the expression (81) of times ${\stackrel{\u0305}{t}}_{{R}_{i}}$, a systematic two by two cancellation of all intermediate terms occurs, leaving only the two ends, and giving a result in simpler form:
$$\stackrel{\u0305}{t}\left({\mathit{r}}_{1},{\mathit{r}}_{2},\dots ,{\mathit{r}}_{n}\right)={\tau}_{b}\left[1p\left({R}_{1}\right)p\left({R}_{2}\right)\dots p\left({R}_{n}\right)\right].$$(84)
The exact result without truncation will be with ${\prod}_{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)$, for the probability product.
In a given realization of the medium, the final mean survival time ${\stackrel{\u0305}{t}}_{{\mathit{r}}_{1}}$ at position r_{1}, will be obtained by repeating the path constructions and taking the average of the above mean trapping times, over all these constructions:
$${\stackrel{\u0305}{t}}_{{\mathbf{r}}_{1}}=\underset{N\to \mathrm{\infty}}{\mathrm{lim}}\frac{\sum \mathrm{}\stackrel{\u0305}{t}\left({\mathit{r}}_{1},{\mathit{r}}_{2},\dots ,{\mathit{r}}_{n}\right)}{N}.$$(85)
The sum corresponds to the systematic repetition of the constructions and N is the number of times the constructions are repeated. The r_{2}, …, r_{n} including the final index value n vary from one construction to another. Likewise in the vectors [R_{1}, R_{2}, …, R_{n}] specifying one construction, only the first value R_{1} repeats itself, all the others, R_{2}, …, R_{n}, vary, including the vector dimension n if ϵ is fixed. The exact result without any truncation would be, where the overline on the right is the average over an infinite set of constructions:
$${\stackrel{\u0305}{t}}_{{\mathit{r}}_{1}}={\tau}_{b}\left[1\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}\right].$$(86)
The microscopic excess temperature field at position r_{1},
$$\tau \left({\mathit{r}}_{1}\right)=\mathrm{i\omega}\left({\beta}_{0}{T}_{0}/{\rho}_{0}{c}_{P}\right)P{\stackrel{\u0305}{t}}_{{\mathit{r}}_{1}},$$(87)
(see Eq. (72)), will be directly obtained geometrically from the construction of N → ∞ vectors [R_{1}, R_{2}, …, R_{n}] by applying (84) and (85), and using the expression (82) of the hitting probability p(R). This determination will correspond to launching N random walkers at r_{1} and making the N different constructions for the given same sample of material. In the limit N → ∞ the mean survival time or equivalently the microscopic field τ at this position is determined.
To determine the dynamic thermal permeability we only need to have a macroscopic average 〈τ〉(x). To have it we can for example proceed as follows.
First consider Lorentz’s conception of the average. We work with one sample of the material and must consider the result of an average in a representative volume. To calculate this volume average 〈τ〉(x), where x is the cental position of the representative volume, we can randomly select an initial position r_{1}, accept it and increment the number of trials N if it lies in the fluid part of this volume, then perform the construction, [R_{1}, R_{2}, …, R_{n}], apply our formula (84), and take the average on repetition of this process. Convergence towards the sought macroscopic average will automatically occur as the number N of initial positions r_{1} increases, without predetermination of the microscopic field. Of course, as the number of walkers increases and the result converges (in $1/\sqrt{N}$), the microscopic field is also automatically extracted at any point with increasing accuracy; but this is an output, not an input, in the averaging process.
With a Gibbs’ conception of the average, the initial position x is fixed and we can randomly select, each time, a new configuration of the traps. We accept it and increment the number of trials N if x lies in the fluid, each time giving rise to one single new construction [R_{1}, R_{2}, …, R_{n}]. The sought average is given by the same formula (84), averaged by repeating the process. Convergence will automatically occur as N increases indefinitely, without any predetermination in any configuration, of the microscopic field. In both these Lorentz and Gibbs cases, as we are dealing with stochastic processes, we know that the convergence will be slow, in $1/\sqrt{N}$, with N the used number of constructions. Before moving on to validation and method examples, it is useful now to restate, further detail, and summarize our findings.
Let us return to the macroscopic average mean life time of diffusing radioactive species that are directly absorbed in the pore wall. We can consider it in the Lorentz conception when it is calculated in a single sample by means of constructions derived from random points of a representative fluid volume, or in the Gibbs conception when the point of origin is invariable and the construction each time is carried out in a new configuration including certain random variations. (Also, mixed conceptions that make partial use of both ways of averaging are possible). By indicating symbolically in all cases, with an overline (which must be interpreted appropriately), the exact way in which the average is performed, one can always designate the mean survival time in question $\stackrel{\u0305}{t}$ as used in equation (86)^{21}:
$$\stackrel{\u0305}{t}={\tau}_{b}\left[1\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}\right].$$(88)
The important thing is that whatever conception is considered, for the overline it will correspond to a certain vision of the macroscopic medium and a definite procedure to generate the construction [R_{1}, R_{2}, …, R_{n},…] and calculate the average $\stackrel{\u0305}{{\prod}_{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}$.
Now recall the correspondence (72)–(74). There follows immediately from it and the equations (88) and (75) that the dynamic thermal permeability is given by:
$${k}^{\mathrm{\prime}}\left(\omega \right)=\frac{{\nu}^{\mathrm{\prime}}\varphi}{\mathrm{i\omega}}\left[1\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}\right].$$(89)
In general, whatever the geometry tortuous or not, the highfrequency limit ${\alpha}_{\mathrm{\infty}}^{\text{'}}$ of α′(ω) is automatically equal to 1. It is customary to express the dynamic thermal tortuosity α′(ω) in terms of its highfrequency limit, 1, and a viscous relaxationfunction, χ′(ω), taking a value of 1 when the thermal effects are in relaxed state (isothermal motions originating from the pore walls having time to fully penetrate the fluid during the cycle time), and a value of 0 when they are in frozen state (isothermal motions having no time to penetrate the fluid during the cycle time):
$$\frac{1}{{\alpha}^{\mathrm{\prime}}\left(\omega \right)}=1{\chi}^{\mathrm{\prime}}\left(\omega \right).$$(90)
A comparison of this definition (90) with (89) and (37), immediately shows us that the averaged infinite product, $\stackrel{\u0305}{{\prod}_{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}$, is nothing but this thermal relaxation function:
$${\chi}^{\mathrm{\prime}}\left(\omega \right)=\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}.$$(91)
Once the relaxation function is conveniently directly determined through the average (91), the other thermal functions derive immediately. The effective bulk modulus, in particular, is given by
$${K}_{f}^{1}\left(\omega \right)={K}_{a}^{1}\left[1+\left(\gamma 1\right){\chi}^{\mathrm{\prime}}\left(\omega \right)\right].$$(92)
In the low frequency limit, it is a known result [9, 41], that α′(ω) can be expanded in the following Laurent series:
$${\alpha}^{\mathrm{\prime}}\left(\omega \right)=\frac{{\nu}^{\mathrm{\prime}}\varphi}{\mathrm{i\omega}{k}_{0}^{\mathrm{\prime}}}+{\alpha}_{0}^{\mathrm{\prime}}+{{\alpha}^{\mathrm{\prime}}}_{1}^{2}\frac{\mathrm{i\omega}{k}_{0}^{\mathrm{\prime}}}{{\nu}^{\mathrm{\prime}}\varphi}+{{\alpha}^{\mathrm{\prime}}}_{2}^{3}{\left(\frac{\mathrm{i\omega}{k}_{0}^{\mathrm{\prime}}}{{\nu}^{\mathrm{\prime}}\varphi}\right)}^{2}+\dots ,$$(93)
where ${k}_{0}^{\mathrm{\prime}}$ is the static thermal permeability or inverse trapping constant, ${\alpha}_{0}^{\mathrm{\prime}}$ is the static thermal tortuosity, and ${\alpha}_{1}^{\mathrm{\prime}}$, ${\alpha}_{2}^{\mathrm{\prime}}$, etc., are some higher order purely geometrical dimensionless parameters. Comparing this expansion with the expansions that follow from using equations (91) and (82) with ${\mu}_{R}=R\sqrt{\mathrm{i\omega}/\nu \mathrm{\prime}}$, as given by (73) and (74), and the known small argument expansion of the function I_{0} or sinh, it is easy to derive explicit geometric expressions for the above geometric parameters. For example in 2D, we find for ${k}_{0}^{\mathrm{\prime}}$ and ${\alpha}_{0}^{\mathrm{\prime}}$ the following purely geometric expressions:
$${k}_{0}^{\mathrm{\prime}}/\varphi =\frac{1}{4}\stackrel{\u0305}{\sum _{i=1}^{\mathrm{\infty}}\mathrm{}{R}_{i}^{2}},$$(94)
and,
$$\frac{{{k}_{0}^{\mathrm{\prime}}}^{2}{\alpha}_{0}^{\mathrm{\prime}}}{{\varphi}^{2}}=\stackrel{\u0305}{\frac{3}{64}\sum _{i=1}^{\mathrm{\infty}}\mathrm{}{R}_{i}^{4}+\frac{1}{16}\sum _{i=1}^{\mathrm{\infty}}\mathrm{}\sum _{j=i}^{\mathrm{\infty}}\mathrm{}{R}_{i}^{2}{R}_{j}^{2}}.$$(95)
We let the reader derive higherorder expressions for higherorder geometric parameters.
With respect to the highfrequency limit, we know that we have an expansion in terms of successive powers of the small thermal skin depth (as long as the pore walls are smooth and are therefore considered as locally flat surfaces at the boundary layer level) [42],
$${\chi}^{\mathrm{\prime}}\left(\omega \right)={\left(\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}+\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\frac{1}{{\mathrm{\Sigma}}^{\mathrm{\prime}}}+{\left(\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\right)}^{3/2}\frac{1}{4{V}^{\mathrm{\prime}}}+\dots ,$$(96)
with thermal characteristic length Λ′, given by, 2/Λ′ = S_{p}/V_{p}, where S_{p} and V_{p} are pore wall surface and pore volume. (The numerical factors have been chosen so that for cylindrical circular pores of radius R, one finds, Λ′ = R, Σ′ = R^{2} and V′ = R^{3}). To our knowledge there are no known general expressions for the higher order parameters, Σ′, V′, etc., which have surface, volume, etc., dimensions. In this limit the expression (91) cannot be expanded because for the first large radii in the product we are in the large argument limit, but for the last very small ones near the pore wall we are in the small argument limit. Nevertheless, given the HF data of the function χ′(ω) calculated using (91), we will have the means to estimate, sequentially, the values of the high frequency parameters, considering the ω → ∞ limits, e.g. for Λ′, Σ′ and V′:
$$\frac{2}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}=\underset{\omega \to \mathrm{\infty}}{\mathrm{lim}}{\left(\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}\right)}^{1/2}\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)},$$(97)
$$\frac{1}{{\mathrm{\Sigma}}^{\mathrm{\prime}}}=\underset{\omega \to \mathrm{\infty}}{\mathrm{lim}}\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}\left[\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}{\left(\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}\right],$$(98)
and,
$$\frac{1}{4{V}^{\mathrm{\prime}}}=\underset{\omega \to \mathrm{\infty}}{\mathrm{lim}}{\left(\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}\right)}^{3/2}\left[\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}{\left(\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\frac{1}{{\mathrm{\Sigma}}^{\mathrm{\prime}}}\right],$$(99)
etc., where the frequencydependent complete expressions of the p(R_{i}) are left. Provided the statistics is sufficient to accurately estimate the averaged infinite product $\stackrel{\u0305}{{\prod}_{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}$ the above highfrequency parameters, Λ′, Σ′, V′, etc., will be directly computable by the simulations from estimation of the above recursive limits.
6 Validation and examples of computation
We can now first validate our algorithm (Sect. 5) by considering a simple case where we have a known analytical solution within the Biot framework of the simplifications made. In particular, and for the special case where the material has airfilled pores aligned in identical cylinders, we can choose a situation where the crosssection has a very simple shape such as a slit, a disk, or equilateral triangle, so that we can figure out the analytical formula. To fix our ideas, consider an aligned cylindrical circular pore of radius R. For long wavelength sound propagation along the zaxis of the tube, and coherent with the incompressible simplification in Biot theory, the complete set of viscothermal differential equations governing fluid motion in a single tube, that normally is expressed in Kirchhoff’s theory [43] and leads to an exact^{22} nonlocal macroscopic description [44], can be simplified in a local theory in the manner done by Zwikker and Kosten [45], which is reexpressed below. Based on the simplification of the typical fluid incompressibility inherent in Biot theory, everything happens as if the following simplified governing equations can be combined:
(i) equations determining the velocity pattern v(x, y) = v(r):
$${\rho}_{0}\frac{\partial v}{\partial t}=\frac{\partial P}{\partial z}+\eta \left(\frac{{\partial}^{2}}{\partial {x}^{2}}+\frac{{\partial}^{2}}{\partial {y}^{2}}\right)v,\hspace{1em}\left(\mathrm{cross}\u2010\mathrm{section}\right),$$(100)
$$v=\mathbf{0},\hspace{1em}\left(\mathrm{pore}\u2010\mathrm{wall}\right),$$(101)
where ∂P/∂z is a given constant over the crosssection,
(ii) equations determining the excess temperature pattern τ(x, y) = τ(r):
$${\rho}_{0}{c}_{P}\frac{\partial \tau}{\partial t}={\beta}_{0}{T}_{0}\frac{\partial P}{\partial t}+\kappa (\frac{{\partial}^{2}}{\partial {x}^{2}}+\frac{{\partial}^{2}}{\partial {y}^{2}})\tau ,\hspace{1em}\left(\mathrm{cross}\u2010\mathrm{section}\right),$$(102)
$$\tau =0,\hspace{1em}\left(\mathrm{pore}\u2010\mathrm{wall}\right),$$(103)
where β_{0}T_{0}∂P/∂t is again a spatial constant over the crosssection, which will be combined with the following additional relations:
(iii) equation of continuity:
$$\frac{\partial b}{\partial t}+\frac{\partial v}{\partial z}=\frac{1}{{\rho}_{0}}\frac{\partial \rho \text{'}}{\partial t}+\frac{\partial v}{\partial z}=0,\hspace{1em}\left(\mathrm{cross}\u2010\mathrm{section}\right),$$(104)
and,(iv) thermodynamic equation of state:
$$\gamma {\chi}_{0}P=b+{\beta}_{0}\tau ,\hspace{1em}\left(\mathrm{cross}\u2010\mathrm{section}\right).$$(105)
In harmonic regime e^{−iωt} the solution fields v(r) and τ(r) are of the following form:
$$v=\frac{\nu}{\mathrm{i\omega}}\left[\frac{\mathrm{\partial}P}{\eta \mathrm{\partial}z}\right]\left[1\frac{{I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{\nu}}r\right)}{{I}_{0}(\sqrt{\frac{\mathrm{i\omega}}{\nu}}R})\right];\hspace{1em}\tau =\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\left[{\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\kappa \mathrm{\partial}t}\right]\left[1\frac{{I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}r\right)}{{I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right)}\right].$$(106)
For the case of other crosssectional shapes of cylindrical tubes, the same differential equations and boundary conditions will apply, but must be solved taking into account the particular shape of the geometry. With 〈 〉 denoting the average at the crosssection of the tube, we can define macroscopic variables such as average axial velocity 〈v〉, average excess temperature 〈τ〉, average overdensity 〈ρ′〉, or average condensation B = 〈b〉 = 〈ρ′〉/ρ_{0}.
Therefore, comparing the above simplified governing equations (possibly written in the general cylindrical case) with those discussed in Section 3, we see that in the harmonic regime we have, at the macroscopic, averaged level, a viscous, dynamic Darcy law, stating that:
$$\langle v\rangle =\frac{k\left(\omega \right)}{\eta}\frac{\mathrm{\partial}P}{\mathrm{\partial}z},$$(107)
where k(ω) can be directly computed by using this definition and the average at the crosssection, of the harmonic regime solution v(x, y). For the cylindrical circular tube, v(x, y) = v(r) which is given by (106), and a straightforward integration shows that:
$$k\left(\omega \right)=\frac{\nu}{\mathrm{i\omega}}\left[12\frac{{I}_{1}\left(\sqrt{\frac{\mathrm{i\omega}}{\nu}}R\right)}{\left(\sqrt{\frac{\mathrm{i\omega}}{\nu}}R\right){I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{\nu}}R\right)}\right].$$(108)
The bracket is just the inverse dynamic viscous tortuosity function 1/α(ω). In general, when the geometry is not trivial like here in cylindrical tubes, but tortuous, it is convenient to express the dynamic viscous tortuosity in terms of its highfrequency nontrivial limit, α_{∞}, and a viscous relaxationfunction, χ(ω), taking a value of 1 when the viscous effects are in relaxed state (vortical motions originating from the pore walls having time to fully penetrate the fluid during the cycle time), and a value of 0 when they are in frozen state (vortical motions having no time to penetrate the fluid during the cycle time):
$$\frac{1}{\alpha \left(\omega \right)}=\frac{1}{{\alpha}_{\mathrm{\infty}}}\left[1\chi \left(\omega \right)\right].$$(109)
For all cylindrical tubes (whatever the form of the perimeter including fractal) α_{∞} is trivially equal to 1. Here, for cylindrical circular tubes of fixed radius, we see, comparing (108) and (109), that the complex expression in the brackets (108) is just the above relaxationfunction χ(ω) for this cylindrical circular case.
In the same way, the comparison of the above simplified governing equations with those of Section 3, also leads us to the analogous thermal and dynamic “Darcy” law, which states that:
$$\langle \tau \rangle =\frac{{k}^{\mathrm{\prime}}\left(\omega \right)}{\kappa}{\beta}_{0}{T}_{0}\frac{\mathrm{\partial}P}{\mathrm{\partial}t},$$(110)
where k′(ω) can be directly computed by using this definition and the average at the crosssection, of the harmonic regime solution τ(x, y). For cylindrical circular tubes of fixed radius R, τ(x, y) = τ(r) is given by (106), and the integration shows that:
$${k}^{\mathrm{\prime}}\left(\omega \right)=\frac{{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}}\left[12\frac{{I}_{1}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right)}{\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right){I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right)}\right].$$(111)
The bracket is just the inverse of the socalled dynamic thermal tortuosity function 1/α′(ω), for this cylindrical circular case. Comparing (111) and (90) we see that the complex expression in the brackets (111) is just the thermal relaxation function for the case of cylindrical circular tubes of radius R:
$${\chi}_{\mathrm{circ}}^{\mathrm{\prime}}\left(\omega \right)=2\frac{{I}_{1}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right)}{\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right){I}_{0}\left(\sqrt{\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}}R\right)}.$$(112)
Note that for the special geometry of aligned cylindrical tubes having any type of crosssection, the viscous and thermal problems are essentially the same and the frequencydependences make appear essentially one and the same function. In this case, viscous and thermal functions are directly transformed in one another by the scaling relations illustrated in the above example, namely, χ′(ω) = χ(ωPr), or, k(ω) = k′(ω/Pr), or α(ω) = α′(ω/Pr), where Pr ≡ ν/ν′ ≡ ηc_{P}/κ is the Prandtl number. In arbitrary geometries, on the contrary, there is no general relationship between the viscous and thermal response functions.
Recalling the general identification (91) of the thermal relaxation function, we see that, to validate our simulation technique in the case of a cylindrical circular tube, we only need to check the following. When we randomly launch the initial position of a random walker inside a disk of radius R, then construct and average the probability product over a large number of trials, it should be observed that:
$$\stackrel{\u0305}{\prod _{i=1}^{\mathrm{\infty}}\mathrm{}p\left({R}_{i}\right)}={\chi}_{\mathrm{circ}}^{\mathrm{\prime}}\left(\omega \right).$$(113)
In Figure 3 we present an Argand plot of this relaxation function, as given by the theoretical closed form solution (112), and as estimated with N = 10^{5} walkers by our random motion simulation technique. For our programming in Matlab here, Lorentz volume averaging was used, where the position of the circular section is considered fixed, while the initial position of the walker is uniformly random within it. The consecutive absolute positions r_{i} are dummy except for their differences, that define the construction steps R_{i} = r_{i+1} − r_{i}. To ensure uniform filling of the section required for volume averaging, it is sufficient to define the first radius R_{1} directly by line code, R(:, 1) = R0*(1 − sqrt(rand(N, 1))), where N is the number of constructs and R0 is the pore radius. Constructions reach the limit (porewall at radius R0) relatively quickly, most often no more than thirty steps. By systematically taking 100 steps in one construction (without testing the final radius) we are almost sure of the success of the operation, rendering the matlab programming trivial. Obviously the conception of Gibbs averaging can be used equivalently, systematically starting the path construction at the coordinate origin and randomly distributing the distance L of the disk region center to the origin, with suitable probability ditribution, namely, L(:, 1) = R0*sqrt(rand(N, 1)), which reproduces the previous uniform filling. It corresponds to making all purely random translations of the section that let the origin lie inside the section.
Figure 3 Relaxation Function of cylindrical circular tubes: Theory, Random motion simulation with N = 10^{5} walkers, and Modeling. 
All calculations were performed for 51 reduced frequencies^{23}, spanning the low, intermediate, and high frequencies. As mentioned, a convergence in $1/\sqrt{N}$ is slow. With N = 10^{5} walkers as in Figure 2, the eye can see without difficulty, that, in the transition region, some of the red circle dots indicating the Brownian motion simulation are still not fully centered on the black crosses indicating the theoretical values. With N = 10^{7} walkers not shown, the eye no longer perceives the differences. Indeed, with N = 10^{7} walkers, the equations (93) and (95) were found to give respectively the value, 0.1249996R^{2}, for the parameter ${k}_{0}^{\mathrm{\prime}}={k}_{0}$, whose theoretical value is, R^{2}/8 = 0.125R^{2}, and the value, 1.3336, for the parameter ${\alpha}_{0}^{\mathrm{\prime}}={\alpha}_{0}$, whose theoretical value is, 4/3 = 1.33333….
Also plotted in the Figure 3, are the values (for the present case) of the following two general models applicable in arbitrary geometries.
Model 1, was derived by Lafarge [4, 18], in clarification of previous work by Allard and Champoux [20]. It is the thermal pure counterpart of Johnson et al. general model [28] of the viscous dynamic tortuositypermeability. For the effective bulk modulus it reads:
$$\frac{{K}_{a}}{{K}_{f}\left(\omega \right)}=\gamma \left(\gamma 1\right){\left[1+\frac{{\nu}^{\mathrm{\prime}}\varphi}{\mathrm{i\omega}{k}_{0}^{\mathrm{\prime}}}{\left(1\frac{4\mathrm{i\omega}{{k}^{\mathrm{\prime}}}_{0}^{2}}{{\nu}^{\mathrm{\prime}}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}{\varphi}^{2}}\right)}^{1/2}\right]}^{1}.$$(114)
The thermal permeability ${k}_{0}^{\mathrm{\prime}}$ and thermal characteristic length Λ′ are independent parameters. However, the following dimensionless ratio
$${M}^{\mathrm{\prime}}=\frac{8{k}_{0}^{\mathrm{\prime}}}{\varphi {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}},$$(115)
is generally found to depart from unity by less than one order of magnitude. For the circular tubes here, M is exactly equal to 1, as it is known results that, ${k}_{0}^{\mathrm{\prime}}={R}^{2}/8\varphi $, and, Λ′ = R. The expression (114) simplifies and happens to coincide with the socalled AllardChampoux simpler model [20]:
$$\frac{{K}_{a}}{{K}_{f}\left(\omega \right)}=\gamma \left(\gamma 1\right){\left[1+\frac{8{\nu}^{\mathrm{\prime}}}{\mathrm{i\omega}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}}{\left(1\frac{\mathrm{i\omega}{{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}}{16{\nu}^{\mathrm{\prime}}}\right)}^{1/2}\right]}^{1}.$$(116)
Thus here, Model 1 is not different from the usual widely used AllardChampoux model that considers the approximation, ${k}_{0}^{\mathrm{\prime}}=\varphi {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8$. Making apparent the ratio M′ in (114) it can be seen that the magnitude of ${k}_{0}^{\mathrm{\prime}}$, or ${{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}$, determines a rollover frequency where a transition occurs between low and highfrequencies. The shape of this transition is affected by the value of the “form factor” M′.
Model 2, is a refined model, also derived by Lafarge [9, 18] in clarification of previous work by Pride et al. [46]. It essentially gives an improved definition of the shape of the transition frequency region, by adding the information provided by the thermal static tortuosity parameter ${\alpha}_{0}^{\mathrm{\prime}}$, which determines an additional form factor,
$${m}^{\mathrm{\prime}}=\frac{{M}^{\mathrm{\prime}}}{4\left({\alpha}_{0}^{\mathrm{\prime}}1\right)}.$$(117)
The refined model replaces the square root in (114) by the expression:
$$1{m}^{\mathrm{\prime}}+{m}^{\mathrm{\prime}}{\left(1+\frac{{M}^{\mathrm{\prime}}}{2{{m}^{\mathrm{\prime}}}^{2}}\frac{\mathrm{i\omega}}{{\nu}^{\mathrm{\prime}}}\frac{{k}_{0}^{\mathrm{\prime}}}{\varphi}\right)}^{1/2}.$$(118)
In doing so it is capable to retrieve the exact two first terms in (93). For an understanding of these models we refer to [4, 18, 28, 42], and more recently [9] (Appendix).
As we can see, Model 2 performs better than Model 1, however, small discrepancies subsist. The discrepancies are magnified in the Argand Plot representation, which particularly emphasizes the transition region. A direct look at the quantity of interest, the effective dynamic apparent modulus as a function of frequency, is given in Figures 4 and 5 where we plot the real and imaginary parts of the normalized bulk modulus versus reduced frequency $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$. (This choice of the reduced frequency is because of the next illustration, see below). On these figures the number of walkers is N = 10^{7} and the Brownian motion simulation cannot be distinguished from the theoretical values.
Figure 4 Real part of K_{f}(ω) versus Reduced frequency Ω = ωR^{2}/8ν′: Theory, Random motion simulation, and Modeling. 
Figure 5 Imaginary part of K_{f}(ω) versus Reduced frequency Ω = ωR^{2}/8ν′: Theory, Random motion simulation, and Modeling. 
What we can conclude from these figures and the numerical findings on ${k}_{0}^{\mathrm{\prime}}$ and ${\alpha}_{0}^{\mathrm{\prime}}$ is that the Brownian motion method is validated. When the number of random walkers is not too high, it is possible to have very fast but not very precise estimates. The results presented above are prepared for the cylindrical circular pores, which is a very simple shape. One could check similarly that for other or more complex shapes for which we have the analytical solution, such as slitlike pores or triangular equilateral pores, a similar perfect matching is also found between the simulations and the analytical formulae.
We now further illustrate the proposed algorithm by calculating the dynamic bulk modulus of air in two different 2D arrays of fibers surrounded by air (see Fig. 6): a regular square array and a “random” array, both at porosity ϕ = 0.97. We put random into quotes because the fibers positions were not completely random as the fibers were impenetrable. As illustrated on the Figure 6, we based our random geometry on a 21 L × 21 L square region in which 441 fibers of radius ${R}_{0}=L\sqrt{(1\varphi )/\pi}$ were accomodated.
Figure 6 The two different, regular and “random”, arrangements of cylinders. There are 21 × 21 = 441 fibers in our representative cell. Each walker is released at the origin, where no fiber is present, and sees a new configuration – see the text for detail. N = 10^{7} trials have been taken. 
The 441 fibers were successively randomly introduced in the region (right part, “Random” array). When a fiber overlapped a preceding one, or when it overlapped the origin of the walkers random motion (central position of the region, and origin of the construction of the ray suite, [R_{1}, R_{2}, …, R_{n}, …] as sketched in Fig. 6), the position was rejected. A new random position was selected, and so on. In the regular case (square array), to save time, we adapt the same general programming scheme, using $441$ regularly spaced fibers but performing Gibbs random positioning of the whole square configuration (not shown), only avoiding overlap of the same central origin with the displaced central fiber (which remains in the central L × L square). In this manner the programming for the regular configuration was adapted in trivial manner from that done for the random configuration.
In both geometries the paths [0, r_{2}, …, r_{n}, …] were found to never go outside the 21 L × 21 L region, with typical drift, before reaching the pore walls, observed to be on the order of L. Very rarely, a drift as large as 8 L was observed in the random geometry. Note that the drift values give an idea of the representative cell dimensions of the material’s microstructure. Note also that, by eliminating the overlap between fibers, we automatically ensure that the porosity ϕ and the thermal characteristic length Λ′ in both geometries, random and regular, are the same (for instance Λ′ is here automatically Rϕ/(1 − ϕ)). Thus, using the reduced frequency, $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$, we know that we should get the exact same frequency dependence for both geometries in the HF asymptotic limit. All calculations are done for the same 51 reduced frequencies as before.
It will be shown in this illustration that the irregular position of the fibers in real glass wool has a direct influence on the thermal relaxation of the effective bulk modulus, i.e. the way of moving from a relaxed isothermal state to a frozen adiabatic state. Note that due to the large density of glass compared to air, a porosity as high as ϕ = 0.97 is still a safe value to ensure that the fibers essentially remain at ambient temperature, justifying the assumed boundary conditions (28). The real part of the dynamic bulk modulus for the regular and “random” settings at porosity ϕ = 0.97 is plotted in Figure 7 versus the reduced frequency. The corresponding imaginary part is plotted in Figure 8. The predictions of Model 1 and Model 2 are also indicated, using as input parameters for ${k}_{0}^{\mathrm{\prime}}$ and ${\alpha}_{0}^{\mathrm{\prime}}$ the values inferred from the random walk algorithm, that were calculated beforehand with the same number of trials, N = 10^{7}, using (94) and (95).
Figure 7 Real part of the normalized bulk modulus K_{f}(ω)/K_{a} versus reduced frequency $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$ for the two different regular and “random” arrangements of cylinders. Brownian motion method is applied here with N = 10^{7} walkers, both for the computation of $\mathfrak{R}{K}_{f}\left(\omega \right)$, and the estimation of parameters ${k}_{0}^{\mathrm{\prime}}/\varphi $ and ${\alpha}_{0}^{\mathrm{\prime}}$ to plot the Model 1 and Model 2. 
Figure 8 Imaginary part of the normalized bulk modulus K_{f}(ω)/K_{a} versus reduced frequency $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$ for the two different regular and “random” arrangements of cylinders. Brownian motion method is applied here with N = 10^{7} walkers, both for the computation of $\mathfrak{I}{K}_{f}\left(\omega \right)$, and the estimation of parameters ${k}_{0}^{\mathrm{\prime}}/\varphi $ and ${\alpha}_{0}^{\mathrm{\prime}}$ to plot the Model 1 and Model 2. 
We note that we have at our disposal prior determinations of the parameters ${k}_{0}^{\mathrm{\prime}}$ and α_{0}, in Cortis’ PhD thesis [47], for the regular square setting at porosity ϕ = 0.97. They were made by using FreeFM finite elements. Here, with N = 10^{7} random walkers and by systematically taking 150 steps for one construction (rare events were observed where the previous value of 100 steps was insufficient to reach the boundary), our values are, ${k}_{0}^{\mathrm{\prime}}/\varphi =0.17141{L}^{2}$ and ${\alpha}_{0}^{\mathrm{\prime}}=1.0660$. Those of Cortis were, ${k}_{0}^{\mathrm{\prime}}/\varphi =0.17144{L}^{2}$ and ${\alpha}_{0}^{\mathrm{\prime}}=1.0654$. The differences concern digits that are not significant within the given FreeFM finite elements calculations.
In the random case the value of the static thermal permability is found to be approximately two times (1.87) that of the regular arrangement (namely 0.3197 L^{2} versus 0.1714 L^{2}). Thus the transition frequency is shifted to lower frequencies. The transition is steeper in the regular case, corresponding to a geometry where microscopic scales are less distributed than in the random case. This effect is clearly better described by an enhanced model where a second form factor (117) allows accounting for the broadening of the transition. Indeed, it can be noted that the value of the thermal tortuosity ${\alpha}_{0}^{\mathrm{\prime}}$ which is a measure of disorder (we have the known formula, ${\alpha}_{0}^{\mathrm{\prime}}=\langle {\tau}_{0}^{2}{\rangle}_{p}/\langle {\tau}_{0}{\rangle}_{p}^{2},$ see [18] or Appendix in [9], where τ_{0} is the relaxed excesstemperature pattern), increases significantly in the random case (1.3833 versus 1.0660), leading to a rather small value of the shape factor m′ (0.1671 versus 0.5205). This is why model 1, which assumes m′ = 1, is significantly less precise in the random case. Model 2 is also less precise, but to a much lesser extent. The introduction of the lowfrequency information ${\alpha}_{0}^{\mathrm{\prime}}$ leads to an improvement of the description of the relaxation transition. Clearly, the proposed random walk method provides an efficient way to numerically calculate quasiexactly the dynamic bulk modulus in a given geometry. Model 2 was found to be very close to the exact result.
7 Conclusion
A very simple simulation technique utilizing the known properties of Brownian motion with radioactive decay has been proposed to estimate, within the framework of Biot theory, the dynamic bulk modulus of a gas filling a porous medium. This technique has been validated on the simplest case of a cylindrical circular tube, and shown to provide a convenient way to calculate thermal relaxation in arbitrary geometries. The case of regular and random parallel fiber arrangement is considered because it corresponds to two extreme cases of glass wool microgeometry. The results of the quasi exact Brownian motion determination of the relaxation pattern of the function K_{f}(ω) highlight a simple relaxation function model ([9], Appendix) in terms of the following three pore space geometric parameters: static thermal permeability ${k}_{0}^{\mathrm{\prime}}$ divided by porosity ϕ, thermal tortuosity ${\alpha}_{0}^{\mathrm{\prime}}$, and thermal characteristic length Λ′. They show how the fiber arrangement is reflected on the shape of the relaxation transition.
The proposed simulation technique can provide useful clues on questions relative to the notion of macroscopic averaging. It is easy to implement and can be directly applied to all problems related to diffusioncontrolled reactions between sinks, where the diffusing species are added in volume by using a spatially uniform, arbitrarily timevariable, source.
In upcoming work, we will demonstrate the performance of this meshfree method by quickly determining the thermal functions in pores with fractal boundaries. In this case, the highfrequency limit (96) is no longer valid because the surface of the pore walls is not smooth, and a revision of the relaxation function model will be necessary. Whereas the present model corresponds to a special CurievonSchweidler law with exponent 1/2 (see Appendix in [9]), the revised model will correspond to a general CurievonSchweidler law with exponent that will depend on the fractal dimension of the pore walls.
In conclusion, we note that given the very general considerations in the original seminal work of Pontrjagin et al. [17], on which this paper is based, and the analogy we have presented between the thermal and viscous/inertial dynamical problems, there might be a question of generalizing the method to viscous/inertial dynamic problems, which, by combining an ideal nonviscous divergencefree flow motion with a Brownian motion with disintegration, would be able to describe the viscous/inertial relaxation. This could be an advance in solving the incompressible dynamic linearized Stokes equations.
Conflict of interest
The authors declare that they have no conflicts of interest in relation to this article.
Acknowledgments
The authors would like to thank Professors Evgeny A. Pchelintsev and Claude Depollier for sending them a copy of the article by Pontrjagin et al. [17] on which the method in this paper is based. N. Nemati acknowledges financial support from a French government Grant managed by ANR within the framework of the National Program Investments for the Future ANR11LABX002201.
Data availability statement
The research data associated with this article are included within the article.
Appendix A
Additional comments on incompressibility, pressure, mean pressure, and averages
As we are concerned here only with rigidframed materials and motion in the air, we introduce a macroscopic averaging operation 〈 〉(x) in the air, as a total–volume coarse graining average, $\langle \xb7\rangle \left(\mathit{x}\right)=\frac{1}{V}{\int}_{{V}_{p}}\mathrm{}\xb7\mathrm{d}{V}_{p}$, given by the local air–phase integral of the considered quantities, divided by the corresponding total (air–phase + solid–phase) volume, where V_{p}(x) is the porespace volume (the fluid volume) and V the total volume, in the representative volume (averaging sphere) of central position x. An averaging with normalization taken with respect to airvolume V_{p} instead of totalvolume V, say $\langle \xb7{\rangle}_{p}\left(\mathit{x}\right)=\frac{1}{{V}_{p}\left(\mathit{x}\right)}{\int}_{{V}_{p}}\mathrm{}\xb7\mathrm{d}{V}_{p}$, where index p over 〈 〉 is for porespace, may be preferable when willing to have averaged values representative of the values in the fluid phase only. (Note that by abuse of language we use the same notation for the fluid domain V_{p}(x) and its corresponding volume; See also the two paragraphs before equation (1) for more details, and our two different notations of microscopic (r) and macroscopic (x) positions: The fluid volume element dV_{p} is d^{3}r whereas x is reserved for the central position of an averaging sphere).
As we assume that Biot theory is a relevant description of the actual physics, the fluid and solid motions are, as insisted in the text, longwavelength quasiincompressible motions at the pore scale (i.e. in an averaging sphere). In particular it means here that, the microscopic excesspressure acoustic field in the saturating fluid (the ambient air, here), p, is very nearly timevariable porescale constant. No important difference is, quantitatively, to be made between the microscopic, p, and spatial averaged macroscopic, P, excesspressure. The quantitative doublescale asymptotic homogenizationtheory justification for this [4, 48–50] is pertinent but not fully exact as explained in Sections 2 and 3. The latter sections provide qualitative justification of the incompressibility and approximate correspondence between p and P (not ∇p and ∇P!), as mere consequences of the absence of very different pore sizes (which excludes local resonance and pressure diffusion process [51]); see also the discussion around equations (7.160)–(7.161) in [9] (Appendix). The very small variations p − P that occur genuinely at the pore scale are removed in the macroscopic average P. Not explicitly mentioned in literature, within Biot’s theory the latter can be synonymously defined as: (i) the local direct coarsegrained average, P = 〈p〉_{p} (here local refers to making the average around a chosen central position x and is not to be confused with the notion of local and nonlocal theory), or, (ii) the local indirect coarsegrained average given by, 〈pv_{i}〉 = P〈v_{i}〉, where v_{i} is the ith component of the microscopic velocity in the fluid, i.e. this time, the HeavisidePoynting or Umov averaged pressure (see Refs. [8, 9]) that can also be termed as a lumped pressure (see [52]). In this indirect definition of P the average is made indifferently with totalvolume normalization, 〈 〉, or with airphase normalization, 〈 〉_{p} (both give same result for P, as we have the general identity, 〈·〉(x) = ϕ(x)〈·〉_{p}(x), with ϕ(x) = V_{p}(x)/V the porosity). Here and in the text, for shortness, we sometimes omit to mention “excess” or “over” before “pressure” for p or P, simply called pressures. We assert that as long as Biot theory remains relevant, i.e. the fluid/solid incompressibility applies, both definitions of macroscopic pressure P(x), direct and indirect, are completely equivalent and indistinguishable. We note that, in what has been said above we have interpreted our averages 〈 〉(x) or 〈 〉_{p}(x) as volume averages performed at some central position x in a given sample (the socalled Lorentz point of view expressed in Refs. [8, 9]). The discussions in Section 3 are conducted using this point of view. In another poweful point of view that can be used in place of the preceding if ergodicity applies (the socalled Gibbs point of view in Refs. [8, 9]), the averages 〈 〉(x) or 〈 〉_{p}(x) are more simply viewed as the expectation value of the considered quantity at the position x, in an ensemble of realizations of the macroscopic medium; the first is when the quantities are extended to zero in the solid and all realizations are taken into consideration, obtaining a total volume normalization, the second is when we discard realizations with x lying in the solid, obtaining a fluidvolume normalization. Our reasonings on quantities such as the effective bulk modulus or the thermal permeability/tortuosity, will be the same whether our averages 〈 〉, 〈 〉_{p}, are conceived in terms of volume average, ensemble average, or a mixture of the two. We note also that the averagingvolume utilized to perform the macroscopic averages can be replaced by an averaging surface (see e.g. Sect. 2 in [28]). Finally, let us mention that, whereas the above direct and indirect definitions of the macroscopic pressure are both possible and equivalent in the context of Biot theory which is a local theory neglecting spatial dispersion, see Sections 1–3, the indirect definition is actually deeper. We anticipate that it is the latter, suitably generalized as: −〈pv_{i}〉 = H_{ij}〈v_{j}〉, that will provide a relevant general definition of a macroscopic symmetric stress field H_{ij} in the fluid, in the case of nonlocal theory capable to describe arbitrary microgeometries. This last comment refers to the mentioned general theory nonlocal in time and space that still has to be developed in its full details as recalled in Sections 1 and 2.
Appendix B
Additional comments on the constraints of Biot’s theory
Loosing sight of the semiempirical origin of Biot theory it is sometimes written (e.g. [48]) that this theory is rigorously derived from the microscopic equations by the method of twoscale asymptotic homogenization developed and studied by SanchezPalencia [53], Keller [54], Bensoussan et al. [55], and others. However, this view is a misinterpretation of the character of this homogenization method, which automatically introduces the assumption of fluid/solid incompressibility as a result of its own procedure that hypothesizes the existence of one characteristic pore size, in microgeometry. The twoscale homogenization method does not derive Biot theory from microscopic equations only, but from the microgeometric constraints on which it is based.
Biot theory, as shown by the previously mentioned absence of spatial dispersion, and the related crucial role played by the incompressibility hypothesis, is not general in some potentially important aspects, such as the definitions with which it can be expressed (like the systematic assimilation of macroscopic values with direct volume averages, see comments in Appendix A), and the quantities involved. It constitutes a special (though widely applicable) simplified limit (nonlocal in time but local in space) of a much more general theory (applicable to much wider classes of microgeometries) nonlocal in both time and space. This more general theory will possibly lead, in other simplified limits, to more general local (local in space) descriptions. The fully general nonlocal theory or its more general local simplifications would not satisfy the characteristic assumption of fluid/solid incompressibility^{24} and would lead to much richer behavior than conventional Biot theory allows.
The characteristic reduction in the possibility of temporal dispersion in Biot model is seen in the fact highlighted by Johnson [28] (Appendix A) and Lafarge [18] and [4] (Appendix C), that the singularities of the response functions α(ω), k(ω), α′(ω), k′(ω), in the complex ωplane, lie only on the imaginary halfaxis, as opposed to the full halfplane allowed by the general causality principle; this, also, expresses in the monotonic behaviors mentioned in Section 3. Using a divergencefree representation of all fluid displacements combined with a uniform representation of all solid displacements, at the pore scale, Biot’s model cannot account for spatial dispersion, at all. The latter remains entirely outside its framework. With spatial dispersion present as soon as significant divergent motions appear at pore scale [9], the usual simple basis for the definition of local response factors (that is detailed in Sect. 3), simply vanishes. In these circumstances the operators and field quantities must be thoroughly generalized, possibly enhancing even their tensor rank (see Appendix A and ^{Footnotes 1} and ^{16}). In a forthcoming work (the first mentioned in ^{Footnote 1}), we will illustrate how an imbalanced divergent fluid motion at the small scale can generate a much more general temporal dispersion in connection with underlying presence of spatial dispersion.
In general it is the relatively good satisfaction of the incompressibility of solid and fluid motion at long wavelengths in a variety of microgeometries that explains the large scope and experimental success of Biot theory. These simplifications are expected to hold true in the real propagation problem and to justify Biot theory, not only when there is a welldefined pore size, but more generally as long as the complications in the geometries do not lead to locally resonant structures [8, 9] or other effects related to fluid flow divergence, such as pressure diffusion [51]. In all these cases, both the longwavelength spatial dispersion (related to the nonuniform or divergent parts of the solid or the fluid motions on the pore scale) and the mentioned additional temporal dispersion (also triggered by these motions) are practically absent. This is an expression of the fact that, compared to the nearly uniform (solid) or nondivergent (fluid) motions on the pore scale in Biot theory, the actual nonuniform/divergent distributed motions of the solid/fluid on the pore scale usually have a very weak amplitude. Their complete absence in Biot theory is usually not a problem, but a reasonable simplification.
Therefore, instead of assuming a very welldefined pore size, as is the case with twoscale asymptotic homogenization, it would be more appropriate to characterize the physical position of Biot theory by recognizing that it corresponds to situations in which one can neglect the longwavelength spatial dispersion phenomena and a corresponding part of the temporal dispersion, which allows considerable simplifications. We note that the fact that the asymptotic twoscale homogenization method automatically preserves, to a first approximation, the incompressibility of the solid and fluid motion at the pore scale is an indication that this method is not suitable to describe the spatial dispersion and part of the temporal dispersion. For this reason, the higher order terms in this method do not accurately describe the onset of the spatial dispersion corrections that occur at decreasing wavelengths, regardless of the microgeometries, simple or complex, and the part of the temporal dispersion effects associated with this onset. This indicates the limited physical consistency of the expansion principle used in this method. Without the unnecessary detour via this method, as made in [4], we can base our description of the function K_{f}(ω, x) directly on the incompressibility of fluid and solid motion, which we introduce here as an explicit hypothesis/simplification justified by the type of geometry considered.
References
 M.A. Biot: The theory of propagation of elastic waves in a fluidsaturated porous solid. I. Low frequency range, II. Higher frequency range. Journal of the Acoustical Society of America 28 (1956) 168–191. [CrossRef] [Google Scholar]
 K. Attenborough: Acoustical characteristics of porous materials. Physics Reports 82 (1982) 179–227. [CrossRef] [Google Scholar]
 J.F. Allard, N. Atalla: Propagation of sound in porous media: Modelling sound absorbing materials. 2nd edn., WileyBlackwell, 2009. [CrossRef] [Google Scholar]
 D. Lafarge, P. Lemarinier, J.F. Allard, V. Tarnow: Dynamic compressibility of air in porous structures at audible frequencies. Journal of the Acoustical Society of America 102 (1997) 1995–2006. [CrossRef] [Google Scholar]
 D.L. Johnson, T.J. Plona: Recent developments in the acoustic properties of porous media. In: D. Sette, Ed. Proceedings of the International School of Physics Enrico Fermi, Course XCIII, vol. 17, North Holland, Amsterdam, 1996, pp. 255–290. [Google Scholar]
 L.D. Landau, E. Lifshitz: Electrodynamics of continuous media. Pergamon Press, Oxford, 1960. [Google Scholar]
 V.M. Agranovich, V.I. Ginzburg: Spatial dispersion in crystal optics and the theory of excitons. Interscience Publishers, London, 1966. [Google Scholar]
 D. Lafarge: Acoustic wave propagation in viscothermal fluids – an electromagnetic analogy. In: N. Jiménez, O. Umnova, J.P. Groby, Eds. Acoustic waves in periodic structures, metamaterials, and porous media, vol. 143, Springer, 2021, pp. 205–272. [CrossRef] [Google Scholar]
 D. Lafarge: Nonlocal dynamic homogenization of fluidsaturated metamaterials, in: N. Jiménez, O. Umnova, J.P. Groby, Eds. Acoustic waves in periodic structures, metamaterials and porous media, vol. 143, Springer, 2021, pp. 273–331; Appendix: Local dynamic homogenization of rigidframed fluidsaturated porous materials, ibid. pp. 316–331. [CrossRef] [Google Scholar]
 N. Nemati, A. Kumar, D. Lafarge, N.X. Fang: Nonlocal description of sound propagation through an array of Helmholtz resonators. Comptes Rendus Mécanique 433 (2015) 656–669. [CrossRef] [Google Scholar]
 L.L. Beranek: Acoustic impedance of porous materials, Journal of the Acoustical Society of America 13 (1942) 248–260. [CrossRef] [Google Scholar]
 C. Truesdell: Precise theory of the absorption and dispersion of forced plane infinitesimal waves according to the NavierStokes equations. Journal of Rational Mechanics and Analysis 2 (1953) 643–741. [Google Scholar]
 S. Torquato, I.C. Kim: Efficient simulation technique to compute effective properties of heterogeneous media. Applied Physics Letters 55 (1989) 1847–1849. [CrossRef] [Google Scholar]
 D. Lafarge: Determination of the dynamic bulk modulus of gases saturating porous media by brownian motion simulation. In: J.L. Auriault, Ed., Poromechanics II, Proceedings of the Second Biot Conference on Poromechanics, Swets and Zeitlinger, Grenoble, 2002, pp. 703–708. [Google Scholar]
 C. Perrot, R. Panneton, X. Olny: Computation of the dynamic thermal dissipation properties of porous media by brownian motion simulation: Application to an opencell aluminum foam. Journal of Applied Physics 102 (2007) 074917. [CrossRef] [Google Scholar]
 S. Lifson, J. Jackson: On the selfdiffusion of ions in a polyelectrolyte solution. Journal of Chemical Physics 36 (1962) 2410–2414. [CrossRef] [Google Scholar]
 L. Pontrjagin, A. Andronow, A. Witt: On the statistical consideration of dynamical systems. Journal of Experimental and Theoretical Physics (in Russian) 3 (1933) 172. [Google Scholar]
 D. Lafarge: Propagation du son dans les matériaux poreux à structure rigide saturés par un fluide viscothermique (Sound propagation in porous materials with rigid structure saturated by a viscothermal fluid). PhD thesis, Université du Maine, 1993. https://cyberdoc.univlemans.fr/theses/1993/1993LEMA1009.pdf. [Google Scholar]
 M. Avellaneda, S. Torquato: Rigorous link between fluid permeability, electrical conductivity, and relaxation times for transport in porous media. Physics of Fluids A: Fluid Dynamics 3 (1991) 2529–2540. [CrossRef] [Google Scholar]
 Y. Champoux, J.F. Allard: Dynamic tortuosity and bulk modulus in airsaturated porous media. Journal of Applied Physics 70 (1991) 1975–1979. [CrossRef] [Google Scholar]
 J.F. Allard, C. Depollier, A. L’Esperance: Observation of the Biot slow wave in a plastic foam of high flow resistance at acoustical frequencies. Journal of Applied Physics 59 (1986) 3367–3370. [CrossRef] [Google Scholar]
 J.F. Allard, Y. Champoux, C. Depollier: Modelization of layered sound absorbing materials with transfer matrices. Journal of the Acoustical Society of America 82 (1987) 1792–1796. [CrossRef] [Google Scholar]
 C. Depollier, J.F. Allard, W. Lauriks: Biot theory and stressstrain equations in porous sound absorbing materials. Journal of the Acoustical Society of America 84 (1988) 2277–2279. [CrossRef] [Google Scholar]
 W. Lauriks, A. Cops, J.F. Allard, C. Depollier, P. Rebillard: Modelization at oblique incidence of layered porous materials with impervious screens. Journal of the Acoustical Society of America 87 (1990) 1200–1206. [CrossRef] [Google Scholar]
 J.F. Allard, C. Depollier, P. Guignouard, P. Rebillard: Effect of a resonance of the frame on the surface impedance of glass wool of high density and stiffness. Journal of the Acoustical Society of America 89 (1991) 999–1001. [CrossRef] [Google Scholar]
 P. Guignouard, M. Meisser, J.F. Allard, P. Rebillard, C. Depollier: Prediction and measurement of the acoustical impedance and absorption coefficient at oblique incidence of porous layers with perforated facings. Noise Control Engineering Journal 36 (1995) 129–135. [Google Scholar]
 J.S. Bolton, N.M. Shiau, Y.J. Kang: Sound transmission through multipanel structures lined with elastic porous materials. Journal of Sound and Vibration 191 (1996) 317–347. [CrossRef] [Google Scholar]
 D.L. Johnson, J. Koplik, D. Dashen: Theory of dynamic permeability and tortuosity in fluidsaturated porous media. Journal of Fluid Mechanics 176 (1987) 379–402. [CrossRef] [Google Scholar]
 A.N. Norris: On the viscodynamic operator in Biot’s equations of poroelasticity. Journal of WaveMaterial Interaction 1 (1986) 365–380. [Google Scholar]
 M.A. Biot, D.G. Willis: The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics 24 (1957) 594–601. [CrossRef] [Google Scholar]
 R.J.S. Brown: Connection between formation factor for electrical resistivity and fluidsolid coupling factor in Biot’s equations for acoustic waves in fuidfilled porous media. Geophysics 45 (1980) 1269. [CrossRef] [Google Scholar]
 D.L. Johnson, J. Koplik, L.M. Schwartz: New poresize parameter characterizing transport in porous media. Physical Review Letters 57 (1986) 2564–2567. [CrossRef] [PubMed] [Google Scholar]
 L.D. Landau, E. Lifshitz: Fluid mechanics. Pergamon Press, 1987. [Google Scholar]
 H.A. Lorentz: The theory of electrons and its applications to the phenomena of light and radiant heat, in: A course of lectures delivered in Columbia University, New York, in March and April 1906, Leipzig, B.G. Teubner, 1909. [Google Scholar]
 L.D. Landau, E. Lifshitz: Mechanics. ButterworthHeinemann, 2000. [Google Scholar]
 L.D. Landau, E. Lifshitz: Statistical physics. Pergamon Press, LondonParis, 1959. [Google Scholar]
 C. Boutin, C. Geindreau: Estimates and bounds of dynamic permeability of granular media, Journal of the Acoustical Society of America 124 (2009) 3576–3593. [Google Scholar]
 D. Lafarge: Porous and stratified porous media linear models of propagation, chapter 6, The equivalent fluid model, in: M. Bruneau, C. Potel, Eds., Materials and acoustics handbook, ISTE, Wiley, 2009, pp. 147–202; see pp. 178–179. [CrossRef] [Google Scholar]
 C. Perrot, F. Chevillotte, R. Panneton, J.F. Allard, D. Lafarge: On the dynamic viscous permeability tensor symmetry. Journal of the Acoustical Society of America 124 (2008) EL210–EL217. [CrossRef] [PubMed] [Google Scholar]
 L.H. Zheng, Y.C. Chiew: Computer simulation of diffusioncontrolled reactions in dispersions of spherical sinks. Journal of Chemical Physics 90 (1989) 322–327. [CrossRef] [Google Scholar]
 R. Roncen, Z.E.A. Fellah, D. Lafarge, E. Piot, F. Simon, E. Ogam, M. Fellah, C. Depollier: Acoustical modeling and bayesian inference for rigid porous media in the lowmid frequency regime. Journal of the Acoustical Society of America 144 (2018) 3084–3101. [CrossRef] [PubMed] [Google Scholar]
 J. Kergomard, D. Lafarge, J. Gilbert: Transients in porous media: Exact and modelled timedomain Green’s functions. Acta Acustica united with Acustica 99 (2013) 557–571. [CrossRef] [Google Scholar]
 G. Kirchhoff: Über des einfluss der wärmeleitung in einem gase auf die schallbewegung. Annalen der Physik 134 (1868) 177–193. English translation: On the influence of thermal conduction in a gas on sound propagation, in: R.B. Lindsay (Ed.), Physical Acoustics, in: Benchmark Papers in Acoustics, Vol. 4, Dowden, Hutchinson and Ross, Stroudsburg, PA, 1974, pp. 7–19. [CrossRef] [Google Scholar]
 N. Nemati, D. Lafarge: Check on a nonlocal Maxwellian theory of sound propagation in fluidsaturated rigid framed porous media. Wave Motion 51 (2014) 716–728. [CrossRef] [Google Scholar]
 C. Zwikker, C.W. Kosten: Sound absorbing materials. Elsevier Publishing Company Inc., New York, 1949. Reprinted 2012 by the NAG (Nederlands Akoestisch Genootschap), 1949. [Google Scholar]
 S.R. Pride, F.D. Morgan, A.F. Gangi: Drag forces of porousmedium acoustics. Physical Review B 47 (1993) 4964–4978. [CrossRef] [PubMed] [Google Scholar]
 A. Cortis: Dynamic acoustic parameters of porous media: a theoretical, numerical and experimental investigation. PhD thesis, Delft University of Technology, 2002. [Google Scholar]
 R. Burridge, J.B. Keller: Poroelastisity equations derived from microstructure. Journal of the Acoustical Society of America 70 (1981) 1140–1146. [CrossRef] [Google Scholar]
 P. Sheng, M.Y. Zhou: Dynamic permeability in porous media. Physical Review Letters 61 (1988) 1591–1594. [CrossRef] [PubMed] [Google Scholar]
 D.M.J. Smeulders, R.L.G.M. Eggels, M.E.H. van Dongen: Dynamic permeability: reformulation of theory and new experimental and numerical data. Journal of Fluid Mechanics 245 (1992) 211–227. [CrossRef] [Google Scholar]
 R. Venegas, C. Boutin: Enhancing sound attenuation in permeable heterogeneous materials via diffusion processes. Acta Acustica United with Acustica 104 (2018) 623–635. [CrossRef] [Google Scholar]
 A.D. Pierce: Acoustics: An introduction to its physical principles and applications. WileyBlackwell, 1989. [Google Scholar]
 E. SanchezPalencia: Nonhomogeneous media and vibration theory, Lecture Notes in Physics, vol. 127, SpringerVerlag, 1980. [Google Scholar]
 J.B. Keller: Effective behavior of heterogeneous media, in: U. Landman, Ed., Statistical mechanics and statistical methods in theory and application. Plenum, New York, 1977, pp. 631–644. [CrossRef] [Google Scholar]
 A. Bensoussan, J.B.L. Lions, C. Papanicolaou: Asymptotic analysis for periodic structures, in: Studies in mathematics and its applications, vol. 5. NorthHolland, Amsterdam, 1978. [Google Scholar]
 T.J. Plona: Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Applied Physics Letters 36 (1980) 259–261. [CrossRef] [Google Scholar]
 D.L. Johnson: Equivalence between fourth sound in liquid He II at low temperatures and the Biot slow wave in consolidated porous media. Applied Physics Letters 37 (1980) 1065–1067. [CrossRef] [Google Scholar]
 T. Levy, E. SanchezPalencia: Equations and interface conditions for acoustic phenomena in porous media. Journal of Mathematical Analysis and Applications 61 (1977) 813–834. [CrossRef] [Google Scholar]
 G. Russakoff: A derivation of the macroscopic Maxwell equations. American Journal of Physics 38 (1970) 1188–1195. [CrossRef] [Google Scholar]
 G. Klein: Mean firstpassage times of Brownian motion and related problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences A211 (1952) 431–443. [Google Scholar]
 N. Nemati, Y.E. Lee, D. Lafarge, A. Duclos, N. Fang: Nonlocal dynamics of dissipative phononic fluids. Physical Review B 95 (2017) 224304. [CrossRef] [Google Scholar]
Cite this article as: Lafarge D. Nemati N. & Vielpeau S. 2023. Brownian motion with radioactive decay to calculate the dynamic bulk modulus of gases saturating porous media according to Biot theory. Acta Acustica, xx, xx.
A forthcoming first work, “Acoustic waves in air or gasfilled structured porous media: new asymptotic tortuosity/compliability and characteristiclengths derived from nonlocal considerations”, will illustrate this by showing how nonlocalities in suitably structured rigidframed materials with smooth pore walls and for propagation along a macroscopic symmetry axis, can result in the highfrequency limit, or the small boundarylayer limit, in a new complex idealfluid tortuosity, an entirely new complex idealfluid “compliability”, and two new real viscous and thermal characteristic lengths; the second work will show how nonlocalities can lead, in the absence of particular symmetries of the material properties and of the type of field considered, to the appearance of even more complex quantities, the bulk modulus of the “effective medium” becoming, in particular, an elastic tensor of the solid fourindex type linking tensors of stresses (see Appendix A and ^{Footnote 13}) and strains with two indices.
A NavierStokesFourier model is used for the microscopic equations in the air, viscousinertial effects involving the NavierStokes equation and thermal effects the socalled Fourier thermal diffusion equation (24).
The interpretation in terms of diffusiondisintegration comes from the analysis in Sections 4 and 5.
Here, condensation significantly varies from point to point on the pore scale in the transition regime, where it distributes from isothermal regime values, b = p/K_{0}, at the pore walls, to values tending to the adiabatic regime, b = p/K_{a}, in the central pore region. These need to be ensembleaveraged or volumeaveraged, into a socalled macroscopic coarsegrained variable, B = 〈b〉, to yield an effective bulk modulus in macroscopic relationship, B = P/K_{f}(ω).
For simplicity of discussion, we express ourselves here as if the material is macroscopically isotropic and the quantities discussed are scalar; the general anisotropic case has been made in Section 3. For the heat exchange problem and the bulk modulus of the fluid that interests us here, K_{f}(ω), the quantities remain scalar in the general anisotropic case. However, we should be aware that this does not hold in the general case when spatial dispersion is introduced (see Appendix A and ^{Footnotes 1} and ^{13}).
To precise the meaning of “ball” or “bounded” see also the signal theory refinement of Lorentz’s notion of averaging brought by Russakoff [59].
If the material is already macroscopically homogeneous and unbounded, the extension is done ipsofacto: it consists in taking for the ${V}_{p}^{\left(\mathrm{ub}\right)}$ and ${S}_{p}^{\left(\mathrm{ub}\right)}$ the pore space V_{p} and the pore wall surface S_{p} of the material itself. If not, we have to perform an extension in the way indicated. This procedure reminds that used in conventional twoscale asymptotic homogenization, which introduces two set of space variables with extension of field quantities from 3 to 6dimensional space [48].
In the general case of arbitrary microgeometry lying outside the realm of Biot theory, we expect to see the notion of macroscopic pressure to be replaced by the notion of a macroscopic effective stress H_{ij}, with an indirect definition, −H_{ij}〈v_{j}〉_{p} = 〈pv_{i}〉_{p} (see Appendix A). When the fluid motion is almost incompressible, H_{ij} reduces to a diagonal form −Pδ_{ij} with P = 〈p〉_{p} and P〈v〉 = 〈pv〉, and this, whether the medium is isotropic or not. This corresponds to the present equivalence of two definitions. But when the microgeometry is such that the fluid motion is not incompressible the equivalence is not maintained. With arbitrary fields the nontrivial tensor character of H_{ij} will be created as a result of spatial dispersion even in media that would be isotropic in the ordinary sense.
We are not aware that this point has been previously clearly made. This would warrant verification and quantification, and is related to the limited physical consistency we expect, of the expansion principle used in conventional twoscale asymptotic homogenization [48], see Appendix B.
Let us note that, long before the introduction of these ideas by one of us in [4, 18], Levy and SanchezPalencia in a work using the twoscale homogenization [58], expressed the effective compressibility of the saturating gas in terms of a dynamic thermal function just proportional to our inverse dynamic thermal tortuosity. This had been done without noting the present parallel between viscous and thermal functions properly defined.
One can check that it would be useless to consider a more complete representation of the pressure in the starting equation (24), such as obtained with, p(t, r) = π(t, r) + [(r − x) · ∇P(t, x) + P(t, x)]. It would give the same relation as that implied in the problem (26)–(28), between the excess mean temperature, 〈τ〉(t), and the excitations, β_{0}T_{0}[∂〈p〉/∂t](t′), t′ ∈ [−∞, t]. Indeed, no mean excess temperature is created by the deviatoric pressure part π which is by construction of zero mean, 〈π〉 = 0, and the same is true of the term (r − x) · ∇P(t, x).
By definition the Laplace variable representation is, ${K}_{f}(t,\mathit{x})={\int}_{\gamma i.\mathrm{\infty}}^{\gamma +i.\mathrm{\infty}}\mathrm{}\frac{\mathrm{ds}}{2\pi i}{e}^{\mathrm{st}}{K}_{f}^{\mathcal{L}}(s,\mathit{x})$, ${K}_{f}^{\mathcal{L}}(s,\mathit{x})={\int}_{0}^{\mathrm{\infty}}\mathrm{}{K}_{f}(t,\mathit{x}){e}^{\mathrm{st}}\mathrm{d}t$, whereas the Fourier representation is, ${K}_{f}(t,\mathit{x})={\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}\mathrm{}\frac{\mathrm{d\omega}}{2\pi}{K}_{f}(\omega ,\mathit{x}){e}^{\mathrm{i\omega}t}$, ${K}_{f}(\omega ,\mathit{x})={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}\mathrm{}{K}_{f}(t,\mathit{x}){e}^{\mathrm{i\omega}t}\mathrm{d}t$; note that in the text here, to lighten the notation we do not make apparent the possible dependence on x, e.g. $k\mathrm{\prime}(\omega ,\mathit{x})=\varphi \left(\mathit{x}\right)D\langle \stackrel{\u0305}{t}\rangle \left(\mathit{x}\right)$, that would correspond to interpreting the V_{p}, S_{p}, in Section 4, as the ${V}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, ${S}_{p}^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, in Section 3.
Einstein’s relation gives the average of the square of the displacement associated with a given diffusion time. The inverse relation gives the average of the diffusion time required to reach a given displacement for the first time. That it has the inverse form of Einstein’s relation was demonstrated by Fürth (1917), see Klein [60].
Note that as an exact nonlocal description it is not limited to longwavelength, see [44]. In complex materials susceptible to be treated in statistical sense with a Gibbs statistical averaging, an exact statistical nonlocal description could also appear in principle that would not be limited to longwavelengths, see in this connection [8, 9, 61].
Note that, given a microgeometry, we can have a feeling whether or not it will fulfill the characteristic incompressibility assumption – such as not involving structures susceptible to resonate or lead to phenomena such as pressure diffusion [51]. But, in an intermediate situation, it may not be easy to make a quantitative estimate of this.
Cite this article as: Lafarge D. Nemati N. & Vielpeau S. 2023. Brownian motion with radioactive decay to calculate the dynamic bulk modulus of gases saturating porous media according to Biot theory. Acta Acustica, 7, 44.
All Figures
Figure 1 Schematic representation of the construction of Brownian motion paths r_{1}, r_{2}, …, r_{n} with succession of characteristic radii [R_{1}, R_{2}, …, R_{n}]. 

In the text 
Figure 2 Possible outcomes of the random walk with disintegration, associated probabilities and survival times. 

In the text 
Figure 3 Relaxation Function of cylindrical circular tubes: Theory, Random motion simulation with N = 10^{5} walkers, and Modeling. 

In the text 
Figure 4 Real part of K_{f}(ω) versus Reduced frequency Ω = ωR^{2}/8ν′: Theory, Random motion simulation, and Modeling. 

In the text 
Figure 5 Imaginary part of K_{f}(ω) versus Reduced frequency Ω = ωR^{2}/8ν′: Theory, Random motion simulation, and Modeling. 

In the text 
Figure 6 The two different, regular and “random”, arrangements of cylinders. There are 21 × 21 = 441 fibers in our representative cell. Each walker is released at the origin, where no fiber is present, and sees a new configuration – see the text for detail. N = 10^{7} trials have been taken. 

In the text 
Figure 7 Real part of the normalized bulk modulus K_{f}(ω)/K_{a} versus reduced frequency $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$ for the two different regular and “random” arrangements of cylinders. Brownian motion method is applied here with N = 10^{7} walkers, both for the computation of $\mathfrak{R}{K}_{f}\left(\omega \right)$, and the estimation of parameters ${k}_{0}^{\mathrm{\prime}}/\varphi $ and ${\alpha}_{0}^{\mathrm{\prime}}$ to plot the Model 1 and Model 2. 

In the text 
Figure 8 Imaginary part of the normalized bulk modulus K_{f}(ω)/K_{a} versus reduced frequency $\mathrm{\Omega}=\omega {{\mathrm{\Lambda}}^{\mathrm{\prime}}}^{2}/8{\nu}^{\mathrm{\prime}}$ for the two different regular and “random” arrangements of cylinders. Brownian motion method is applied here with N = 10^{7} walkers, both for the computation of $\mathfrak{I}{K}_{f}\left(\omega \right)$, and the estimation of parameters ${k}_{0}^{\mathrm{\prime}}/\varphi $ and ${\alpha}_{0}^{\mathrm{\prime}}$ to plot the Model 1 and Model 2. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.