Brownian motion with radioactive decay to calculate the dynamic bulk modulus of gases saturating porous media according to Biot theory

– We present a new stochastic simulation method for determining the long-wavelength 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 simpli ﬁ cation that the ﬂ uid and solid motions are quasi-incompressible 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 diffusion-disintegration-controlled problem that considers Brownian motion of diffusing particles undergoing radioactive-type decay in the pore volume and instant decay at the pore walls. By appropriately choosing the decay time and the diffusion coef ﬁ cient, the stochastic algorithm we develop to determine the average lifetime of the diffusing particles, directly gives the effective apparent modulus of the saturating ﬂ uid. 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 ﬁ brous 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.


Introduction
One of the elements of Biot's macroscopic linear theory of sound propagation [1], as extended for air-saturated porous materials by Attenborough [2] and Allard and Atalla [3], is the complex effective modulus of air, K f (x), 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 Àixt the excess macro-scopic pressure P, to the mean condensation: P/K f (x, x) = hq 0 i p /q 0 , where hq 0 i p is the pore-scale-averaged excess density in the air, and q 0 is the air equilibrium density (q 0 + q 0 the actual pore-scale variable air density q).See Appendix A and discussion in Section 3 for the definition of total-volume and pore-space-volume averages hÁi and hÁi p and the notion of macroscopic pressure.In the time domain, this implies an operator or integral relation of the form, P ðt; xÞ ¼ Kf ½hq 0 i p =q 0 ðt; xÞ ¼ Z t À1 K f ðtÀt 0 ; xÞ ½hq 0 i p =q 0 ðt 0 ; xÞdt 0 , where, using by analogy a language familiar in electromagnetism [6][7][8][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 (x) or Laplace K L f ðsÞ representation.A general theory that would contradict conventional Biot's theory would have to deal with underlying non-local 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 non-local 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 well-known 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 non-local 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 low-frequency 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 quasi-isothermal expansion and compression [11].Therefore in this limit K f (x) will tend to the air isothermal bulk modulus K 0 (the atmospheric pressure).
Conversely, in a high-frequency limit, heat exchange does not has time to occur, except in the immediate vicinity of the pore walls (inside a so-called thermal boundary layer of characteristic thickness (2m 0 /x) 1/2 , where m 0 = j/q 0 c P is the thermal diffusivity [4], j 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 ever-shrinking vicinity.In this limit K f (x) tends to the adiabatic bulk modulus K a = cK 0 (with c = 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 (low-frequency 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 m 0 of air 2 .
More generally, within Biot theory and when the saturating fluid is not air but any tri-variate fluid [12] with c 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 diffusion-disintegration-controlled problem 4 , which can be stated and solved in the pore spacesee 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 closed-form 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 (x) can become significantly cumbersome: first, the local governing equations for the Fourier diffusion-disintegration 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 computewith a possible shunting of the local field determination stepthe effective modulus K f (x) 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 three-dimensional reconstructed unit cells of aluminum foams saturated by the ambient air.Here we describe precisely the Biot context of the effective fluid bulk-modulus 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 dynamic-viscous and dynamic-thermal tortuosities/ permeabilities response functions emerge within Biot's simplifications, and show how the thermal ones determine the effective fluid bulk-modulus.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 low-frequency 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 high-frequency limit can be computed in another recursive process taking a limit x ? 1.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 well-known closed-form 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 long-proposed frequency-dependent model [18], recently re-exposed ( [9], Appendix), that includes geometric parameters, in addition to porosity /, called (static) thermal permeability k 0 0 (the inverse trapping constant [19], Lafarge [18]) (static) thermal tortuosity a 0 0 (Lafarge [18]), and thermal characteristic length K 0 (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 long-wavelength semi-empirical model allowing to describe low-amplitude 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 averaging-volume, 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 motionssignificantly distributed at the pore scale for the fluid because of the tortuous microgeometry, but not for the solid as it cannot slide like a fluidare nearly divergence-free 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 fluid-saturated 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][21][22][23][24][25][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 divergence-free motion for the fluid and quasi uniform motion of the solid, at the pore scale, i.e. in a representative averaging-volume.
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 low-frequency limit where Biot theory was first formulated, the important parameters are porosity /, Darcy permeability k 0 8 , the static tortuosity a 0 [9,18,29], and the abstract Biot-Willis elastic constants P, Q, R, N [5,30].The latter elastic constants are identified through the three quasi-static Biot-Willis 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 Experimentsone 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 overpressurethe remaining macroscopic coefficients P, Q, R are expressed by the so-called Biot-Willis 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 quasi-static 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 quasi-static thought experiments.The Biot-Willis 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 frequency-dependent values (in practice complex amplitude values with so-called 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 (x) for the original real K f .The fact that the incompressibility hypothesis is the reason for the transfer of the static Biot-Willis experiments into the harmonic domain with the above-mentioned substitutions makes the importance of this hypothesis clear.
At arbitrary frequencies it has been described how the parameters k 0 and a 0 combine in a single notion, the frequency-dependent dynamic permeability/tortuosity function, k(x) and a(x) (see [18] and Appendix in [9]), well-studied in Johnson et al.'s landmark works, e.g.[28] (a(x) m//(Àixk(x)), m = g/q 0 , with g dynamic viscosity, m kinematic viscosity).In particular, Johnson et al., 1986, clarified the high-frequency limit properties of these functions which gave rise to the classical [31] notion of ideal fluid tortuosity a 1 , and, when the pore wall surface is smooth, to the important notion of characteristic viscous or electrical length K [28,32].As explained in [28], the dynamic tortuosity a(x) directly determines the various effective densities (q 11 , q 22 , q 12 ), that appear in Biot theory.It is a complex frequency-dependent function to account in detail for the inertial-viscous interaction between the solid and the fluid.Explicitly, the densities verify, q 11 + q 12 = (1 À /)q s , q 22 + q 12 = /q 0 , and, q 12 = À(a(x) À 1)/q 0 , where q s is the density of the elastic solid on which the frame is built; q 12 describes the inertial-viscous 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 a(x) and k(x) are strongly constrained by the condition that the fluid motion at the pore scale is divergence-free 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 low-frequency and high-frequency geometric parameters (see Appendix A in [28] and also, Appendix in [9]).Actually, the above relations on the Biot effective densities q 11 , q 12 , q 22, are obtained on the condition that, not only the fluid motion is divergence-free, 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 divergence-free 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 rigid-frame motionless structure of same microgeometry as the actual, undeformed, material.
Like the densities and for the same reasons the function K f (x) is unique, whether the frame oscillations characteristic of Biot theory manifest themselves or not.It can be likewise studied for a rigid-frame 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 non-negligible thermal expansion.This is the case of gases in general, such as air here.Gases have an expansion-compression behavior well described by the perfect gas law, and that means they have a coefficient of thermal expansion, b 0 À[oq/(qoP)]P 0 , very close to 1/T 0 (T 0 , absolute ambient temperature, P 0 , atmospheric ambient pressure, q, the density).By the general thermodynamic relation, c , adiabatic speed of sound), that can be found e.g. in Landau and Lifshitz [33], it means that the factor c À 1 deviates significantly from zero (general tri-variate fluid [12]).For liquids, this is not the case because b 0 is very small and, as the thermodynamic relation shows, the deviation is a quadratic effect on b 0 .Simple kinetic theory considerations give: c -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 c À 1 % 0.4).In this case, the isothermal K 0 and adiabatic K a = cK 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 0 (x) and a 0 (x), 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 (x) [9,18] constructed using parameters k 0 0 //, a 0 0 , K 0 , and based on these general considerations, that will be checked in our numerical simulations in Section 6.
Historically, the low-frequency 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 inertial-viscous 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 (x) 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.

Definition of the dynamic bulk-modulus 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 (x) 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 non-local 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.1-6) 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, coarse-grained 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 point-to-point changes that can be meaningfully observedsee the general remarks of Lorentz [34] pp.133-134, § § 113 and 114.That is, in addition to a long-wavelength 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 solid-fluid averaging ball (b) be a bounded spherical region with a connected pore space V ðbÞ p (filled by the fluid) and a solid-fluid pore wall surface S ðbÞ p .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 hvi(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, r ¼ x@=@xþ ŷ@=@y þ ẑ@=@z, which is also written in a classical shorthand notation, r = o/ox (e.g., Landau and Lifshitz, Mechanics [35], § 5, The Lagrangian for a system of particles), where x ¼ xx þ y ŷ þ 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, r will mean o/or when acting on a microscopic variable, and r = o/ox when acting on a coarse-grained variable such as P or hvi, 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 = hpi or hpvi = Phvi).
Since the actual microscopic air motion in the spherical spatial averaging region satisfies, at any time, the equation, (which is one among the complete set of governing coupled equations mentioned above (7.1-6) in [9], g and f, first and second viscosities, and p a thermodynamic variable as explained in [8]), with boundary conditions: and since in Biot theory we apply a divergence-free 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: with boundary conditions: 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 well-defined problem.However, a well-defined 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 ðbÞ p and S ðbÞ p 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 stationary-random extension of the internal structure over the entire space, which generates the regions denoted below V ðubÞ p and S ðubÞ p , where "(ub)" stands for "unbounded".When the material is macroscopically inhomogeneous, these regions should be called V ðubÞ p ðxÞ and S ðubÞ p ðxÞ 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 ðubÞ p ðxÞ and S ðubÞ p ðxÞ are macroscopically homogeneous in nature 12 .unbounded, the extension is done ipso-facto: it consists in taking for the V ðubÞ p and S ðubÞ p 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 two-scale asymptotic homogenization, which introduces two set of space variables with extension of field quantities from 3-to 6-dimensional space [48].
Let ÀrP(t 0 , x) for t 0 2 [À1, 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 ðubÞ p ðxÞ and S ðubÞ p ðxÞ, 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 0 , r) = P + p, to be split in a uniformly variable abstract macroscopic part P defined from the values ÀrP(t 0 , x) by: P(t 0 , r) = (r À x) Á rP(t 0 , x) + P(t 0 , x) (note that the knowledge of ÀrP (t 0 , x) does not give that of P(t 0 , x) but the latter value will play no role here), and a deviatoric part, p(t 0 , r), that will automatically be stationary-random or periodic as a function of r, to be coherent with the spatial constancy of the gradient of the above-defined abstract P(t 0 , r), (ÀoP (t 0 , r)/or = ÀrP(t 0 , x)), which constancy is itself the result of the incompressibility of the considered motion.Obviously this condition that p must be stationary-random or periodic must hold in same manner whether we take, on the abstract field p or abstract fields p, v, the definition /P = hpi or the definition Phvi = hpvi, 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 ÀrP, is given by specified values.With one or the other definition, the equations of the problem are at all times t 0 2 [À1, t]: with boundary conditions: with only the difference that the stationary-random p(t 0 , r) solution, that is unique up to an additive integration constant C(t 0 ), can be fixed in one or the other case either with hpi = 0 or with hpvi = 0.However, as these two ways to remove the additive constant are obviously equivalent we must conclude that the same solution, v(t 0 , r), p(t 0 , r), is obtained, and thus the same macroscopic pressure P is obtained, in both the direct and indirect averaging-conceptions 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., ÀrP ¼ f ðtÞêðtÞ, with arbitrary amplitude function f(t) and arbitrary time-dependent orientation of the unit vector ê, ê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 half-axis [28] (Appendix A), as opposed to the full half-plane that is allowed in pure causality considerations [6,36].This becomes clear in Avellaneda and Torquato's, 1991, relaxation-times representation of the viscous dynamical permeability [19]: their relaxation times Â n correspond to purely imaginary singular angular frequencies, x 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 p, will appear in another abstract problem extended as above: with boundary conditions: 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 0 2 [À1, t], numerically equal to the previously considered pressure gradient: In this second problem, the quantity p, 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, p, 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 + p in the first problem, p = 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 divergence-free 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 p, 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 hvi 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 Dirac-delta stirring force, set along a given direction ê 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 âij ðxÞ, or response factors a ij (t, x), such that, or equivalently, local "permeability" macroscopic operators kij ðxÞ, or response factors k ij (t, x), such that, (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 hi p performed in the averaging-volume of central position x where the values f ðtÞêðtÞ 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 a 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 frequency-dependent functions, dynamic inertial/viscous tortuosity: and dynamic viscous/inertial permeability: At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem: with boundary conditions: and ÀrP, a given spatial vector constant.From the three different unique solutions v (l) associated with the source term ÀrP ¼ À dP d' êðlÞ oriented along the three different mutually orthogonal directions êðlÞ of the coordinate axes (ê ðlÞ i ¼ d il , Kronecker symbol; dP variation of pressure along distance d' taken along êðlÞ ), we get, after averaging in the fluid and using the relations/definitions ( 16) and ( 17) the values of the desired functions k ij (x, x) or a À1 ij ðx; xÞ: The obtained functions a ij (x, x) or k ij (x, x) are automatically symmetric [38,39], expressing an Onsager-Casimir symmetry property, and they are related by: The frequency-dependent dynamic permeability relation (17) between the macroscopic flow velocity hvi = /hvi 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 excess-temperature-pattern s in air, appearing in the above finite coarse-graining-ball of volume V ðbÞ p and pore surface S ðbÞ p centered at x.It obeys the following Fourier equation in V ðbÞ p (which is one among the complete set of governing coupled equations mentioned above (7.1-6)[9]), with boundary conditions: 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 divergence-free 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 averaging-volume can be viewed as the sum of a macroscopic part P (defined with /P = hpi or equivalently hpvi = Phvi), which is uniformly variable (P(t, r) = (r À x) Á rP(t, x) + P(t, x), where rP(t, x) is a vector constantit has null derivative o/or), and a deviatoric part p À P, with zero mean value (in either sense, h(p À P)i = 0 or h(p À P)vi = 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 s(t, r) in the region of central position x and thus its macroscopic mean hsi(t, x), which is the quantity of interest to us, can be viewed as responding to the history of macroscopic terms b 0 T 0 oP/ot 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 ðbÞ p and S ðbÞ p 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: 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 ðubÞ p and S ðubÞ p , periodic or stationary random), with boundary conditions: 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 p field did not exactly mimic the deviatoric overpressure, p À P, i.e. the pore-scale 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 s (and thus its macroscopic mean) in the original sphere, and the course of the macroscopic term b 0 T 0 oP/ot 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 16 Let us note that, long before the introduction of these ideas by one of us in [4,18], Levy and Sanchez-Palencia in a work using the two-scale 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. 17One 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) = p(t, r) + [(r À x) Á rP(t, x) + P(t, x)].It would give the same relation as that implied in the problem ( 26)- (28), between the excess mean temperature, hsi(t), and the excitations, b 0 T 0 [ohpi/ot](t 0 ), t 0 2 [À1, t].Indeed, no mean excess temperature is created by the deviatoric pressure part p which is by construction of zero mean, hpi = 0, and the same is true of the term (r À x) Á rP(t, x).
damped normal modes (Avellaneda and Torquato's viscous relaxation-times analysis is easily transposed in a thermal relaxation-times analysis, see [18]).
By writing the convolution we will learn that there are local macroscopic operators of "thermal tortuosity" â0 ðxÞ, associated with response, kernel factors a 0 (t, x), such that, or equivalently, local macroscopic operators of "thermal permeability" k0 ðxÞ, associated with response, kernel factors k 0 (t, x), such that, The above relations are written with h i p performed in the averaging-volume of the central position x where the f(t) values of the macroscopic source term b 0 T 0 oP/ot have been recorded.Recall that the values of the material parameters, porosity / and kernels a 0 and k 0 , 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 frequency-dependent functions, dynamic thermal tortuosity a 0 (x, x): and dynamic thermal permeability k 0 (x, x): At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem: with boundary conditions: From the unique solution s of this diffusion-controlled 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 a 0 (x, x) or k 0 (x, x): which are related by: The frequency-dependent dynamic thermal permeability relation (32) between the macroscopic excess temperature hsi = /hsi 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 q 0 /q 0 where q 0 is the excess density, and excess temperature: where the notation v 0 is employed for the adiabatic compressibility 1/K a (setting, s = 0, in (38), the resulting relation, cv 0 p = b, must describe the isothermal excess pressure-condensation relation, so that, obviously, cv 0 = 1/K 0 ).Taking the average and the time derivative it gives the relation, since hpi p = P, In harmonic regime e Àixt and working with the complex amplitudes we can use (31) to express the right-hand side in (39).Using the general thermodynamic identity seen in Section 2 and suppressing the common factor (Àix) this gives after straightforward calculation the relation, Therefore in harmonic regime we find a relation having the form, with It gives the Fourier coefficients of the kernel functions v À1 f ðt; xÞ or K f (t, x) which appear in the arbitrary timevariable-regime, operator or integral equation, We recognize here, the relation announced in Section 1.The function v f (x, x), resp.its inverse K f (x, 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, b(x, x), such that, in harmonic regime, It is given by, Note that in (46) we have directly written that, ohbi/ot = ÀrÁhvi, using the identity hr Á vi = r Á hvi, 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 = o/or, and in the right, r = o/ox).
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 s(r) in thermal problem defined in (33)- (35) and the mean survival time tðrÞ 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 tðrÞ 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 tðrÞ 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 ðubÞ p ðxÞ, S ðubÞ p ðxÞ, with the dependence (x) not indicated, for simplicity.Then, with the substitution of V ðubÞ p , S ðubÞ p , 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 s b , from 1 (no disintegration, low frequency limit) to 0 (instant disintegration, high frequency limit), will allow us to explore the effective bulk-modulus function in its Laplace-Transform variable representation, K L f ðsÞ, from s = 0 to s = 1, where s is related to x according to s = Àix using the convention e Àixt .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 s b in the volume V p , i.e., if present (in bulk) at time t the particle has probability j b dt = dt/s 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 diffusion-disintegration process which is assumed to mimick Brownian motion and radioactive decay laws.We define the cumulative probability function Dðr; tÞ (replacing the Lifson and Jackson cumulative probability function Wð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 "D"for "Disappearance" (or "Decay", "Disintegration", "Delete" or Death")qualification for this probability is appropriate if any particle touching the surface S p , or the excitation it carries, is instantly removed; Dðr; tÞ only differs from Wð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 Wðr; tÞ in [16], still hold for Dðr; tÞ.Evidently, for a point s on the surface S p , Dðs; tÞ ¼ 1 for all t.For an interior point we have Dðr; 1Þ ¼ 1, while Dðr; 0Þ ¼ 0. A difference with [16] is that the condition, lim now becomes, after accounting for the finite rate of radioactive decay, lim To see this, recall that Dðr; sÞ is the probability that a particle (or its excitation), starting from r at t = 0, dies before t = s, either by touching S p , or by spontaneously decaying in the bulk.So to calculate D, we launch N ? 1 particles at r at t = 0, and, after a time interval s, we see N 00 paths of particles that have not touched S p or died in the bulk.The remaining N 0 = N À N 00 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, Dðr; sÞ ¼ lim N !1 ðN 0 =N Þ.When s ?0 there is not enough time to reach S p , so the non-dead path amounts to N 00 predicted on the only basis of mass "radioactive" absorption: N 00 = N(1 À j b s); consequently N 0 = Nj b s and D ¼ j b s, as stated in (49).
We can now relate (as Lifson and Jackson did for probability W) the probability D at time s, to its value at the next time t + s.For this purpose we define the function i(r, t, q), such that, i(r, t, q)dq, is the transition probability for a particle, to go from point r, in V p , at time t = 0, to a volume element dq about another point, q, 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 t(r, t, q)dq in [16], because, as our particle now decays in a manner typical of radioactive decay, i is therefore reduced compared to its counterpart t without radioactive decay: We may call Lðr; tÞ Z V p iðr; t; qÞdq the "Life 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 i(r, t, q) is normalized in the sense that, From ( 51) and ( 49) we also note that, It is now clear that the following integral equation holds: The meaning of the above equation is as follows: If a particle has reached the surface S p during the time interval (0, t + s) 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, s), or else, it gets to the point q at time s 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 (s, s + t) 18 .We now recall known results for the successive moments of the ordinary transition probability t 0 (r, t, q) in an infinite medium, which will immediately translate into identical results for the successive moments of i(r, t, q).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: where f is the hydrodynamic friction coefficient and t 0 (r, s, q) 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, lim s!0 To see this we note that the effect of the, force-driven, constant drift velocity F i /f, which brings in (q i À r i ) a term sF i /f which superposes to a random, non-driven one, will bring a cross term, s À1 (F i F j )s 2 /f 2 , which does not contribute in the limit s ?0, and two non-cross-terms, that disappear by the absence of direction of the random nondriven motion.Thus, no force-driven effect appears in the second moment, and the transition probability t 0 (r, t, q) can be replaced by the one without force, Gaussian in each direction, ¼ e ÀjqÀrj 2 =2dDt =ð2pdDtÞ, where d is the dimension.The result then expresses Einstein's translational Brownian motion formula, q i À r i ð Þ q j À r j À Á ¼ 2Dsd ij .For all higher moments it can be checked that, ð58Þ 18 Recall that in all reasonings we can imagine replacing the actual V p , S p , by the V ðbÞ p , S ðbÞ p , or more precisely, the V ðubÞ p , S ðubÞ p .
Now, in view of the fact that the function i(r, t, q) must approach the ordinary (i.e.infinite medium) transition probability as s approaches to zero, the moments of i(r, t, q) must approach the moments of the ordinary transition probability.We then have, from equations ( 54), (55), and (57), We now use the integral equation for Dðr; tÞ (Eq.( 53)), to derive a partial differential equation.Expanding Dðq; tÞ about the point r we have Utilizing the various limiting properties of i and D presented in ( 49), ( 52), ( 60), ( 61), (62), we get a differential equation for the cumulative probability Dð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: This interprets as the result obtained for W in Lifson and Jackson [16] equation (9), modified by the presence of a term Àj b L added in the left-hand side to account for an additional rate of change coming from the radioactive disintegration.Now observe that the equation, Dðr; t þ dtÞ ¼ Dðr; tÞ þ dtoDðr; tÞ=ot; 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 dtoDðr; tÞ=ot 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: From this we finally obtain the following desired differential equation satisfied by tðrÞ: with boundary condition: When the field F is zero one may directly write the solutions for a line, cylinder, and sphere.One has, with notation l a a ffiffiffiffiffiffiffiffiffiffi ffi Three-dimensional case: radius of sphere R, r 2 [0, R], 5 Efficient algorithm to compute K f ðxÞ 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: 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 and next, substituting in the calculations that are made to evaluate the product Dhti, the values D = m 0 and j b = Àix of the coefficients D and j b .Likewise, we can compute the bulk modulus K f (x) by, first, writing the relation (e.g.combine (43) with (75) and the correspondence (72)-( 74)): and next, substituting the above values of D and j b in the expressions obtained for the product j b hti.In fact, it means that in the Laplace-variable representation of the kernel K f (t) of the bulk modulus operator Kf , (45) 19 , we simply have, , with this time, D = m 0 and j 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 t and its average hti over position or over configuration (see in Appendix A, the respective viewpoints of Lorentz and Gibbs).
For the sake of simplicity of Figure 1 and the discussion here we take our material to be an air-saturated 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 tðrÞ of random walkers released in one configuration at a given place, or to calculate its average hti 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 zig-zag 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 (t ðsÞ R , s for 19 By definition the Laplace variable representation is, in the text here, to lighten the notation we do not make apparent the possible dependence on x, e.g.k 0 ðx; xÞ ¼ /ðxÞDhtiðxÞ, that would correspond to interpreting the V p , S p , in Section 4, as the V ðubÞ p ðxÞ, S ðubÞ p ðxÞ, in Section 3. surface), to reach R, will be given by the inverse of Einstein's relation, i.e. 20 , (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 t ðbÞ R (b for bulk).Note that t ðbÞ 1 is the average survival time s b of the walker in the infinite medium.
Let 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: The above two types of outcomes allow us to write it, also, in the following form: To get another relationship between the same two, but yet unknown quantities, p(R) and t ðbÞ R , 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 s b .A proportion 1 À p(R) of them do not reach radius R and live the average time t ðbÞ R .The remaining proportion p(R) reach R and then will later disappear at some instant.The latter particles therefore live first the average time t ðsÞ R to reach R, and then take the additional average time s b to disappear.Thus, we have the obvious relation, Combined with the preceding, it reads, 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: 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 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: 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 t Ri , a systematic two by two cancellation of all intermediate terms occurs, leaving only the two ends, and giving a result in simpler form: The exact result without truncation will be with Q 1 i¼1 pðR i Þ, for the probability product.
In a given realization of the medium, the final mean survival time t 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: 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: The microscopic excess temperature field at position r 1 , (see Eq. ( 72)), will be directly obtained geometrically from the construction of N ? 1 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 ? 1 the mean survival time or equivalently the microscopic field s at this position is determined.
To determine the dynamic thermal permeability we only need to have a macroscopic average hsi(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 hsi(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 pre-determination of the microscopic field.Of course, as the number of walkers increases and the result converges (in 1= ffiffiffiffi N p ), 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 pre-determination 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= ffiffiffiffi N p , 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 t as used in equation (86)21 : 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 Q 1 i¼1 pðR i Þ.Now recall the correspondence (75)-(77).There follows immediately from it and the equations (91) and (78) that the dynamic thermal permeability is given by: In general, whatever the geometry tortuous or not, the high-frequency limit a 0 1 of a 0 (x) is automatically equal to 1.It is customary to express the dynamic thermal tortuosity a 0 (x) in terms of its high-frequency limit, 1, and a viscous relaxation-function, v 0 (x), 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): A comparison of this definition (90) with ( 89) and (37), immediately shows us that the averaged infinite product, Q 1 i¼1 pðR i Þ, is nothing but this thermal relaxation function: 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 In the low frequency limit, it is a known result [9,41], that a 0 (x) can be expanded in the following Laurent series: where k 0 0 is the static thermal permeability or inverse trapping constant, a 0 0 is the static thermal tortuosity, and a 0 1 , a 0 2 , etc., are some higher order purely geometrical dimensionless parameters.Comparing this expansion with the expansions that follow from using equations ( 91) and (82) with l R ¼ R ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Àix=m 0 p , 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 0 and a 0 0 the following purely geometric expressions: and, We let the reader derive higher-order expressions for higherorder geometric parameters.
With respect to the high-frequency 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], with thermal characteristic length K 0 , given by, 2/K 0 = 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, K 0 = R, R 0 = R 2 and V 0 = R 3 ).To our knowledge there are no known general expressions for the higher order parameters, R 0 , V 0 , 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 v 0 (x) calculated using (91), we will have the means to estimate, sequentially, the values of the high frequency parameters, considering the x ? 1 limits, e.g. for K 0 , R 0 and V 0 : and, etc., where the frequency-dependent complete expressions of the p(R i ) are left.Provided the statistics is sufficient to accurately estimate the averaged infinite product Q 1 i¼1 pðR i Þ the above high-frequency parameters, K 0 , R 0 , V 0 , etc., will be directly computable by the simulations from estimation of the above recursive limits.

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 air-filled pores aligned in identical cylinders, we can choose a situation where the cross-section 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 z-axis 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 non-local 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): where oP/oz is a given constant over the cross-section, (ii) equations determining the excess temperature pattern s(x, y) = s(r): where b 0 T 0 oP/ot is again a spatial constant over the cross-section, which will be combined with the following additional relations: (iii) equation of continuity: (iv) thermodynamic equation of state: In harmonic regime e Àixt the solution fields v(r) and s(r) are of the following form: For the case of other cross-sectional 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 h i denoting the average at the cross-section of the tube, we can define macroscopic variables such as average axial velocity hvi, average excess temperature hsi, average over-density hq 0 i, or average condensation B = hbi = hq 0 i/q 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: where k(x) can be directly computed by using this definition and the average at the cross-section, 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: The bracket is just the inverse dynamic viscous tortuosity function 1/a(x).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 high-frequency nontrivial limit, a 1 , and a viscous relaxation-function, v(x), 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): For all cylindrical tubes (whatever the form of the perimeter including fractal) a 1 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 relaxation-function v(x) 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: where k 0 (x) can be directly computed by using this definition and the average at the cross-section, of the harmonic regime solution s(x, y).For cylindrical circular tubes of fixed radius R, s(x, y) = s(r) is given by (106), and the integration shows that: The bracket is just the inverse of the so-called dynamic thermal tortuosity function 1/a 0 (x), 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: Note that for the special geometry of aligned cylindrical tubes having any type of cross-section, the viscous and thermal problems are essentially the same and the frequency-dependences 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, v 0 (x) = v(xPr), or, k(x) = k 0 (x/Pr), or a(x) = a 0 (x/Pr), where Pr m/m 0 gc P /j 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: 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 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 (pore-wall 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.All calculations were performed for 51 reduced frequencies 23 , spanning the low, intermediate, and high frequencies.As mentioned, a convergence in 1= ffiffiffiffi N p 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 0 ¼ k 0 , whose theoretical value is, R 2 /8 = 0.125R 2 , and the value, 1.3336, for the parameter a 0 0 ¼ a 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 tortuosity-permeability.For the effective bulk modulus it reads: The thermal permeability k 0 0 and thermal characteristic length K 0 are independent parameters.However, the following dimensionless ratio 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 0 ¼ R 2 =8/, and, K 0 = R.The expression (114) simplifies and happens to coincide with the so-called Allard-Champoux simpler model [20]: Thus here, Model 1 is not different from the usual widely used Allard-Champoux model that considers the approximation, k 0 0 ¼ /K 0 2 =8.Making apparent the ratio M 0 in (114) it can be seen that the magnitude of k 0 0 , or K 0 2 , determines a rollover frequency where a transition occurs between low-and high-frequencies.The shape of this transition is affected by the value of the "form factor" M 0 .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 a 0 0 , which determines an additional form factor, The refined model replaces the square root in (114) by the expression: 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 X ¼ xK 0 2 =8m 0 .(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.
What we can conclude from these figures and the numerical findings on k 0 0 and a 0 0 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 slit-like 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 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð1 À /Þ=p p were accomodated.
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 K 0 in both geometries, random and regular, are the same (for instance K 0 is here automatically R//(1 À /)).Thus, using the reduced frequency, X ¼ xK 0 2 =8m 0 , 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 0 and a 0 0 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).
We note that we have at our disposal prior determinations of the parameters k 0 0 and a 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 0 =/ ¼ 0:17141L 2 and a 0 0 ¼ 1:0660.Those of Cortis were, k 0 0 =/ ¼ 0:17144L 2 and a 0 0 ¼ 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 a 0 0 which is a measure of disorder (we have the known formula [18] and Appendix in [9], a 0 0 ¼ hs 2 0 i p =hs 0 i 2 p , where s 0 is the relaxed excess-temperature pattern), increases significantly in the random case (1.3833 versus 1.0660), leading to a rather small value of the shape factor m 0 (0.1671 versus 0.5205).This is why model 1, which assumes m 0 = 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 low-frequency information a 0 0 leads to an improvement of the description of the relaxation transition.Clearly, the proposed random walk method provides an efficient way to numerically calculate quasi-exactly the dynamic bulk modulus in a given geometry.Model 2 was found to be very close to the exact result.

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 (x) highlight a simple relaxation function model ( [9], Appendix) in terms of the following three pore space geometric parameters: static thermal permeability k 0 0 divided by porosity /, thermal tortuosity a 0 0 , and thermal characteristic length K 0 .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 diffusion-controlled 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 mesh-free method by quickly determining the thermal functions in pores with fractal boundaries.In this case, the high-frequency 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 Curie-von-Schweidler law with exponent 1/2 (see Appendix in [9]), the revised model will correspond to a general Curie-von-Schweidler 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 divergence-free 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.definition of P the average is made indifferently with totalvolume normalization, h i, or with air-phase normalization, h i p (both give same result for P, as we have the general identity, hÁi(x) = /(x)hÁi 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 h i(x) or h i p (x) as volume averages performed at some central position x in a given sample (the so-called 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 so-called Gibbs point of view in Refs.[8,9]), the averages h i(x) or h i 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 fluid-volume normalization.Our reasonings on quantities such as the effective bulk modulus or the thermal permeability/tortuosity, will be the same whether our averages h i, h i p , are conceived in terms of volume average, ensemble average, or a mixture of the two.We note also that the averaging-volume 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: Àhpv i i = H ij hv j i, that will provide a relevant general definition of a macroscopic symmetric stress field H ij in the fluid, in the case of non-local theory capable to describe arbitrary microgeometries.This last comment refers to the mentioned general theory non-local in time and space that still has to be developed in its full details as recalled in Sections 1 and 2.

Figure 3 .
Figure 3. Relaxation Function of cylindrical circular tubes: Theory, Random motion simulation with N = 10 5 walkers, and Modeling.

Figure 6 .
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 configurationsee the text for detail.N = 10 7 trials have been taken.

Figure 7 .
Figure 7. Real part of the normalized bulk modulus K f (x)/K a versus reduced frequency X ¼ xK 0 2 =8m 0 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 RK f ðxÞ, and the estimation of parameters k 0 0 =/ and a 0 0 to plot the Model 1 and Model 2.

Figure 8 .
Figure 8. Imaginary part of the normalized bulk modulus K f (x)/K a versus reduced frequency X ¼ xK 0 2 =8m 0 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 IK f ðxÞ, and the estimation of parameters k 0 0 =/ and a 0 0 to plot the Model 1 and Model 2.