© The Author(s), Published by EDP Sciences, 2023

1 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(ω), which accounts for heat exchange between the air and the solid constituents taking place at the microscopic level [4] in time harmonic regime. In case where the material can have macroscopic gradients of properties due to inhomogeneous microstructure, a dependence on macroscopic position x will also be present. The index f is for fluid, here the air: recall that in Biot’s theory one has a solid and fluid phase, each forming its own infinite cluster, resp. noted with indexes s and f, and whose coupled macroscopic displacement motions are followed separately [5]. This modulus is of concern in many airborne noise control applications when a precise prediction of sound reflection and transmission versus frequency is needed [3, 4]. Let us define it as the factor relating in harmonic regime e^−iωt the excess macroscopic pressure P, to the mean condensation: P/K_f(ω, x) = 〈ρ′〉_p/ρ₀, where 〈ρ′〉_p is the pore-scale-averaged excess density in the air, and ρ₀ is the air equilibrium density (ρ₀ + ρ′ the actual pore-scale variable air density ρ). See Appendix A and discussion in Section 3 for the definition of total-volume and pore-space-volume averages 〈·〉 and 〈·〉_p and the notion of macroscopic pressure. In the time domain, this implies an operator or integral relation of the form, $P(t,\mathit{x})={\widehat{K}}_f[\langle \rho \mathrm{\prime }{\rangle }_p/{\rho }_0](t,\mathit{x})={\int }_{-\mathrm{\infty }}^t\mathrm{}{K}_f(t-t\mathrm{\prime },\mathit{x}) [\langle \rho \mathrm{\prime }{\rangle }_p/{\rho }_0](t\mathrm{\prime },\mathit{x})\mathrm{d}t\mathrm{\prime }$, where, using by analogy a language familiar in electromagnetism [6–9], we see the presence of temporal dispersion effects but no spatial dispersion (nonlocalities in time but not in space). Actually, this absence indicates that some drastic simplifications of principle have been made in the Biot model. It is important to realize that the effective bulk modulus kernel, or the operator itself, owes its existence to these simplifications, as do other elements of the theory. Here, we will clearly establish what these simplifications are, and exploit them to devise a very simple method of calculation of the bulk modulus kernel in its Fourier K_f(ω) or Laplace $K_f^{\mathcal{L}}\left(s\right)$ representation. A general theory that would contradict conventional Biot’s theory would have to deal with underlying non-local aspects and could lead to very new quantities and operators¹. 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 (quasi-stationary 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(ω) will tend to the air isothermal bulk modulus K₀ (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 (2ν′/ω)^1/2, where ν′ = κ/ρ₀c_P is the thermal diffusivity [4], κ the thermal conductivity and c_P the specific heat coefficient at constant presssure, of air). The air will undergo adiabatic expansion and compression except in this ever-shrinking vicinity. In this limit K_f(ω) tends to the adiabatic bulk modulus K_a = γK₀ (with γ = c_P/c_V ≈ 1.4, the ratio of specific heat coefficients of air at constant pressure or constant volume). In between, at intermediate frequencies, a smooth relaxation transition is observed, which connects the “relaxed” state (low-frequency or “static” isothermal limit), to the “frozen” state (high-frequency adiabatic limit), and is determined by the geometry and dimensions of the pore space and the thermal diffusivity ν′ of air².

More generally, within Biot theory and when the saturating fluid is not air but any tri-variate fluid [12] with γ significantly departing from 1, meaning a gas as opposed to a liquid, the effective fluid bulk modulus is in principle entirely determined by the solution of a “Fourier”³ diffusion-disintegration-controlled problem⁴, which can be stated and solved in the pore space – see Appendix in [9] and the comprehensive presentation in Section 3, leading to equations (33)–(37), equations (29)–(30) and equations (42)–(43). In general, because of the complications in the geometries, no 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(ω) 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⁵. As noted by Torquato and Kim [13] in a larger context, this can be a wasteful calculation process as there is a significant amount of information that is lost in passing from the local to the averaged fields. To quote them: “such calculations can become exorbitant even when performed on a supercomputer”.

In this paper, we present a meshless stochastic numerical algorithm to perform the above operations and compute – with a possible shunting of the local field determination step – the effective modulus K_f(ω) that appears in Biot theory. The idea of the algorithm directly inspired by Torquato and Kim [13] is not new, see [14], but was given in the latter Ref. with one faulty reasoning (and a resulting faulty Eq. (9)). After due correction of the latter, it was, through personal communications, applied by Perrot et al. [15] to the description of thermal dissipation properties of 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⁶. 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 ω → ∞. 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^{\text{'}}$ (the inverse trapping constant [19], Lafarge [18]) (static) thermal tortuosity ${\alpha }_0^{\text{'}}$ (Lafarge [18]), and thermal characteristic length Λ′ (Allard 1991 [20]). Finally, we give conclusions in Section 7.

2 Special physical position of Biot theory and limited scope of the present work

Biot’s (1956) classical linear theory [1], complemented to account for thermal effects [2, 3] and losses in the solid [3], is a 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 motions – significantly distributed at the pore scale for the fluid because of the tortuous microgeometry, but not for the solid as it cannot slide like a fluid – are nearly 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⁷. It was extended and used by Attenborough [2], Allard and collaborators [20–26], and by others, e.g., Bolton and collaborators [27], as the basis for a general description of sound propagation in many common porous acoustic materials saturated by ambient air (polyurethane foams, glass wool, rock wool, etc.), used either singly or in multiple layers with other materials to effectively absorb sound, in noise control applications [3].

Here, we wish to highlight the cornerstone importance, in Biot theory, of the simplification of “quasi-incompressibility” (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 frequency-dependent 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₀⁸, the static tortuosity α₀ [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 Experiments – one in which a representative volume, jacketed and drained with small capillaries, is subjected to hydrostatic overpressure, and one in which, without a jacket, it is directly subjected to the same hydrostatic overpressure – the remaining macroscopic coefficients P, Q, R are expressed by the 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(ω) 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₀ and α₀ combine in a single notion, the frequency-dependent dynamic permeability/tortuosity function, k(ω) and α(ω) (see [18] and Appendix in [9]), well-studied in Johnson et al.’s landmark works, e.g. [28] (α(ω) ≡ νϕ/(−iωk(ω)), ν = η/ρ₀, with η dynamic viscosity, ν 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 α_∞, and, when the pore wall surface is smooth, to the important notion of characteristic viscous or electrical length Λ [28, 32]. As explained in [28], the dynamic tortuosity α(ω) directly determines the various effective densities (ρ₁₁, ρ₂₂, ρ₁₂), 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, ρ₁₁ + ρ₁₂ = (1 − ϕ)ρ_s, ρ₂₂ + ρ₁₂ = ϕρ₀, and, ρ₁₂ = −(α(ω) − 1)ϕρ₀, where ρ_s is the density of the elastic solid on which the frame is built; ρ₁₂ 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 α(ω) and k(ω) 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 ρ₁₁, ρ₁₂, ρ_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 frequency-dependent 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(ω) 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, β₀ ≡ −[∂ρ/(ρ∂P)]P₀, very close to 1/T₀ (T₀, absolute ambient temperature, P₀, atmospheric ambient pressure, ρ, the density). By the general thermodynamic relation, $\gamma -1=T_0{\beta }_0^2{c}_0^2/c_P$, ($c_0=\sqrt{K_a/{\rho }_0}$, adiabatic speed of sound), that can be found e.g. in Landau and Lifshitz [33], it means that the factor γ − 1 deviates significantly from zero (general tri-variate fluid [12]). For liquids, this is not the case because β₀ is very small and, as the thermodynamic identity shows, the deviation is a quadratic effect on β₀. Simple kinetic theory considerations give: γ – 1 = 2/n_d, which is not small, where $n_d$ is the number of excited degrees of freedom of a typical molecule (for air, mainly composed of diatomic molecules, n_d ≈ 5 and γ − 1 ≈ 0.4). In this case, the isothermal K₀ and adiabatic K_a = γK₀ bulk moduli differ significantly and there is a frequency dependence of the effective fluid bulk modulus that must be taken into account as already explained in the Introduction.

The frequency dependence of the effective fluid bulk modulus is directly obtained through the introduction of thermal dynamic permeability/tortuosity functions k′(ω) and α′(ω), kind of thermal counterparts of the notions of viscous dynamic permeability/tortuosity [4, 9, 18]. Johnson et al.’s entire analysis of the kind of frequency dependence that appears in these functions [28] translates directly for them, e.g. [4] (Appendix C) and [9] (Appendix). It is the model of K_f(ω) [9, 18] constructed using parameters k′₀/ϕ, α′₀, Λ′, 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(ω) directly on incompressibility, but we should not forget that this does not hold in the general case. We refer to the Appendices A and B for additional comments on the physical position of Biot theory and the present study.

3 Definition of the dynamic 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(ω) 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 observed – see 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⁹. 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¹⁰) 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_p^{\left(b\right)}$ (filled by the fluid) and a solid-fluid pore wall surface $S_p^{\left(b\right)}$. Underlying Biot’s local theoretical approach is the notion that the fluid velocity pattern v(t, r) in this region of central position x and thus its macroscopic mean 〈v〉(t, x), which is the quantity of interest to us, can be viewed as responding to the history of macroscopic pressure gradients applied (in the present and in the past) at this very position¹¹.

Before proceeding, let us distinguish two possible manifestations of the spatial derivative operator, $\nabla =\widehat{\mathbf{x}}\partial /\partial x+\widehat{\mathbf{y}}\partial /\partial y+\widehat{\mathbf{z}}\partial /\partial z$, which is also written in a classical shorthand notation, ∇ = ∂/∂x (e.g., Landau and Lifshitz, Mechanics [35], § 5, The Lagrangian for a system of particles), where $\mathbf{x}=x\widehat{\mathbf{x}}+y\widehat{\mathbf{y}}+z\widehat{\mathbf{z}}$. For the sake of clarity we adopt two different notations for a position vector x: r for a microscopic position vector indicating a position in the pore space, and x for a macroscopic position vector referring to a macroscopic quantity. With this in mind, in what follows, ∇ will mean ∂/∂r when acting on a microscopic variable, and ∇ = ∂/∂x when acting on a coarse-grained variable such as P or 〈v〉, where x denotes the “central position” of the averaging sphere that is used to define the macroscopic quantity in question (with, for the specific case of macroscopic pressure in the fluid, either of the definitions that we mentioned in Appendix A, ϕP = 〈p〉 or 〈pv〉 = P〈v〉).

Since the actual microscopic air motion in the spherical spatial averaging region satisfies, at any time, the equation,

$${\rho }_0\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}=-\nabla p+\eta \Delta \mathit{v}+\left(\frac{\eta }{3}+\zeta \right)\nabla \left(\nabla ·\mathit{v}\right), \hspace{1em}\mathrm{in} V_p^{\left(b\right)},$$(1)

(which is one among the complete set of governing coupled equations mentioned above (7.1-6) in [9], η and ζ, first and second viscosities, and p a thermodynamic variable as explained in [8]), with boundary conditions:

$$\mathit{v}=\mathbf{0}, \hspace{1em}\mathrm{on} S_p^{\left(b\right)},$$(2)

and since in Biot theory we apply a 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:

$${\rho }_0\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}=-\nabla p+\eta \Delta \mathit{v},\hspace{1em}\mathrm{in} V_p^{\left(b\right)},$$(3)

$$\nabla ·\mathit{v}=0,\hspace{1em}\mathrm{in} V_p^{\left(b\right)},$$(4)

with boundary conditions:

$$\mathit{v}=\mathbf{0}, \hspace{1em}\mathrm{on} S_p^{\left(b\right)}.$$(5)

Because of the lack of definition of variable p in (3) (which is not the same variable as p in (1), and in particular is no longer a thermodynamic variable), and boundary conditions on the remaining surfaces bounding the averaging sphere in the fluid, these equations do not correspond to a 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_p^{\left(b\right)}$ and $S_p^{\left(b\right)}$ into the space outside the sphere in an appropriately defined “uniform manner”, as if the material were macroscopically homogeneous. Using the periodic or stationary-random 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_p^{\left(\mathrm{ub}\right)}$ and $S_p^{\left(\mathrm{ub}\right)}$, where “(ub)” stands for “unbounded”. When the material is macroscopically inhomogeneous, these regions should be called $V_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and $S_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ because their construction refers specifically to the macroscopic position x at which a certain type of local structure of the material is assumed. By construction, the regions $V_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and $S_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ are macroscopically homogeneous in nature¹².

Let −∇P(t′, x) for t′ ∈ [−∞, t] define the values that the physical macroscopic pressure gradients have taken at the x position, and recall that we can define these gradients using either of the definitions of macroscopic pressure that we mentioned in Appendix A. Stated and solved in the construct unrestricted domains $V_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$ and $S_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, the equations (3)–(5) have a unique abstract mathematical solution for the velocity v(t, r), that will be a valid representation of the physical velocities in a representative volume around x, if we constrain the abstract pressure, p(t′, r) = P + π, to be split in a uniformly variable abstract macroscopic part P defined from the values −∇P(t′, x) by: P(t′, r) = (r − x) · ∇P(t′, x) + P(t′, x) (note that the knowledge of −∇P(t′, x) does not give that of P(t′, x) but the latter value will play no role here), and a deviatoric part, π(t′, r), that will automatically be 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′, r), (−∂P(t′, r)/∂r = −∇P(t′, x)), which constancy is itself the result of the incompressibility of the considered motion. Obviously this condition that π must be 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 = 〈p〉 or the definition P〈v〉 = 〈pv〉, to obtain our representative macroscopic pressure P. Finally, with one or the other definition, we solve the same problem where the spatially uniform, arbitrarily time variable forcing −∇P, is given by specified values. With one or the other definition, the equations of the problem are at all times t′ ∈ [−∞, t]:

$${\rho }_0\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}=-\nabla \pi +\eta \Delta \mathit{v}-\nabla P,\hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(6)

$$\nabla ·\mathit{v}=0, \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(7)

$$\pi : \mathrm{stationary} \mathrm{random}/\mathrm{periodic}, \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(8)

with boundary conditions:

$$\mathit{v}=\mathbf{0}, \hspace{1em}\mathrm{on} S_p^{\left(\mathrm{ub}\right)},$$(9)

with only the difference that the stationary-random π(t′, r) solution, that is unique up to an additive integration constant C(t′), can be fixed in one or the other case either with 〈π〉 = 0 or with 〈πv〉 = 0. However, as these two ways to remove the additive constant are obviously equivalent we must conclude that the same solution, v(t′, r), π(t′, r), is obtained, and thus the same macroscopic pressure P is obtained, in both the direct and indirect averaging-conceptions of Appendix A¹³. We note also that, the spatial constancy of the considered macroscopic pressure gradient is a signature that there can be no account of spatial dispersion in the description to be developed. Since this constancy is directly related to the incompressibility condition in the posed mathematical problem, it shows the already mentioned connection between the assumption of incompressibility at the pore level and the rejection of spatial dispersion at macroscopic level.

As for the time variations of the macroscopic pressure gradient, they can be assumed to be arbitrary (as long as this is compatible with the long wavelength considerations), i.e., $-\nabla P=f\left(t\right)\widehat{\mathit{e}}\left(t\right)$, with arbitrary amplitude function f(t) and arbitrary time-dependent orientation of the unit vector $\widehat{\mathit{e}}$, ${\widehat{\mathit{e}}}^2=1$. It should be clear that, if we allow for arbitrary time variation here, this in no way means that the effects of temporal dispersion fully appear in the description. Rather, as mentioned in Section 2, the allowed temporal dispersion effect is very strictly limited. The use of the incompressibility assumption not only precludes the occurrence of spatial dispersion effects, but also severely restricts the temporal dispersion effects: as previously mentioned in Section 2 it causes the zeros and singularities of the response functions in the complex frequency plane to lie on the imaginary 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, ω_n = −i/Θ_n. It follows that these functions are constructed with basis elements which have strictly monotone behavior on the real frequency axis (see [18, 37] and Appendix in [9]).

Note that, exactly the same response fields, v and π, will appear in another abstract problem extended as above:

$${\rho }_0\frac{\mathrm{\partial }\mathit{v}}{\mathrm{\partial }t}=-\nabla \pi +\eta \Delta \mathit{v}+\mathit{f},\hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(10)

$$\nabla ·\mathit{v}=0, \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(11)

$$\pi : \mathrm{stationary} \mathrm{random}/\mathrm{periodic}, \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(12)

with boundary conditions:

$$\mathit{v}=\mathbf{0}, \hspace{1em}\mathrm{on} S_p^{\left(\mathrm{ub}\right)},$$(13)

and in which the imposed spatial constant vector f, instead of representing a macroscopic pressure gradient, now represents an imposed external, bulk, body force, directly acting in the fluid. The same motion as before will result when this spatial constant force will be taken, at the different instants t′ ∈ [−∞, t], numerically equal to the previously considered pressure gradient: $\mathit{f}\left(t\mathrm{\prime }\right)=f\left(t\mathrm{\prime }\right)\widehat{\mathit{e}}\left(t\mathrm{\prime }\right)=-\nabla P(t\mathrm{\prime },\mathit{x})$. In this second problem, the quantity π, in a way, represents the pressure itself, which appears in the fluid in the given nontrivial microgeometry, in reaction to the imposed, past and present, forces f. In the preceding problem, π, in a way, was the local deviation to the macroscopic excess pressure. We say “in a way”, because, in these mathematical problems, since the fluid is viewed as completely incompressible which is a purely mental view, the quantities that seem to act as actual physical overpressures (p = P + π in the first problem, p = π in the second) have only symbolic meaning, not physical meaning (remember that the thermodynamic meaning has been lost). We cannot expect these abstract mathematical quantities to yield a truly quantitatively meaningful representation of the physical excess pressure in the actual physical problem of sound propagation. Recall that the unphysical spatial uniformity of the source terms, related to the 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 π, small, abstract pressure variations, will necessarily be truncated representations of the physical ones; they will be sufficient for our purpose of accurately representing the velocity pattern and computing its spatial average 〈v〉 in a representative volume but we believe that they will not very accurately capture the pattern of the pressure¹⁴.

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 $\widehat{\mathit{e}}$ in the stationary random or periodic medium). Of course, this solution exhibits the limitations indicated in [28] (Appendix A), in that it is constructed with purely damped normal modes¹⁵. By explicitly writing the convolution we will learn that there are local “tortuosity” macroscopic operators ${\widehat{\alpha }}_{\mathrm{ij}}\left(\mathit{x}\right)$, or response factors α_ij(t, x), such that,

$${\rho }_0\frac{\mathrm{\partial }\langle {\mathit{v}}_i{\rangle }_p}{\mathrm{\partial }t}\left(t\right)=-{\widehat{\alpha }}_{\mathrm{ij}}^{-1}\left(\mathit{x}\right){\nabla }_{j}P\left(t\right)=-{\int }_{-\mathrm{\infty }}^t\mathrm{}\mathrm{d}{t}^{\mathrm{\prime }}{\alpha }_{\mathrm{ij}}^{-1}\left(t-t^{\mathrm{\prime }},\mathit{x}\right){\nabla }_{j}P\left(t^{\mathrm{\prime }}\right),$$(14)

or equivalently, local “permeability” macroscopic operators ${\widehat{k}}_{\mathrm{ij}}\left(\mathit{x}\right)$, or response factors k_ij(t, x), such that,

$$\varphi \left(\mathit{x}\right)\langle {\mathit{v}}_i{\rangle }_p\left(t\right)=-\frac{{\widehat{k}}_{\mathrm{ij}}\left(\mathit{x}\right)}{\eta }{\nabla }_{j}P\left(t\right)=-\frac{1}{\eta }{\int }_{-\mathrm{\infty }}^t\mathrm{}\mathrm{d}{t}^{\mathrm{\prime }}{k}_{\mathrm{ij}}(t-t^{\mathrm{\prime }},\mathit{x}){\nabla }_{j}P\left(t^{\mathrm{\prime }}\right),$$(15)

(with symmetry on the exchange of i = 1, 2, 3 and j = 1, 2, 3, see below; the general case of an anisotropic medium is considered, the indices i and j refer to a particular arbitrary choice of three orthogonal coordinate axes; the summation convention on repeated indices is used). The above relations are written with 〈 〉_p performed in the averaging-volume of central position x where the values $f\left(t\right)\widehat{\mathit{e}}\left(t\right)$ of the macroscopic pressure gradient were recorded. Recall that this volume served to construct a representative periodic or stationary random homogeneous medium, so that the values of the material parameters, porosity ϕ and kernels α_ij or k_ij, depend on x if the existing starting material has gradients of properties.

In harmonic regime the operator relations (14) and (15) become multiplications with the hat quantities becoming the usual frequency-dependent functions, dynamic inertial/viscous tortuosity:

$$-\mathrm{i\omega }{\rho }_0{\alpha }_{\mathrm{ij}}(\omega ,\mathit{x})\langle {\mathit{v}}_j{\rangle }_p=-{\nabla }_{i}P,$$(16)

and dynamic viscous/inertial permeability:

$$\varphi \left(\mathit{x}\right)\langle {\mathit{v}}_i{\rangle }_p=-\frac{k_{\mathrm{ij}}\left(\omega ,\mathit{x}\right)}{\eta }{\nabla }_{j}P.$$(17)

At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem:

$$-\mathrm{i\omega }{\rho }_0\mathit{v}=-\nabla \pi +\eta \Delta \mathit{v}-\nabla P,\hspace{1em} \mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(18)

$$\nabla ·\mathit{v}=0,\hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(19)

$$\pi : \mathrm{stationary} \mathrm{random}/\mathrm{periodic},\hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(20)

with boundary conditions:

$$\mathit{v}=\mathbf{0},\hspace{1em}\mathrm{on} S_p^{\left(\mathrm{ub}\right)},$$(21)

and −∇P, a given spatial vector constant. From the three different unique solutions v^(l) associated with the source term $-\nabla P=-\frac{\delta P}{\delta \mathcal{l}}{\widehat{\mathit{e}}}^{\left(l\right)}$ oriented along the three different mutually orthogonal directions ${\widehat{\mathit{e}}}^{\left(l\right)}$ of the coordinate axes (${\widehat{\mathit{e}}}_i^{\left(l\right)}={\delta }_{\mathrm{il}}$, Kronecker symbol; δP variation of pressure along distance $\delta \mathcal{l}$ taken along ${\widehat{\mathit{e}}}^{\left(l\right)}$), we get, after averaging in the fluid and using the relations/definitions (16) and (17) the values of the desired functions k_ij(ω, x) or ${\alpha }_{\mathrm{ij}}^{-1}(\omega ,\mathit{x})$:

$$k_{\mathrm{ij}}(\omega ,\mathit{x})=-\frac{\eta \varphi \left(\mathit{x}\right)}{\delta P/\delta \mathcal{l}}\langle {\mathit{v}}^{\left(j\right)}{\rangle }_p·{\widehat{\mathit{e}}}^{\left(i\right)}, \hspace{1em}{\alpha }_{\mathrm{ij}}^{-1}(\omega ,\mathit{x})=-\frac{\left(-\mathrm{i\omega }{\rho }_0\right)}{\delta P/\delta \mathcal{l}}\langle {\mathit{v}}^{\left(j\right)}{\rangle }_p·{\widehat{\mathit{e}}}^{\left(i\right)}.$$(22)

The obtained functions α_ij(ω, x) or k_ij(ω, x) are automatically symmetric [38, 39], expressing an Onsager-Casimir symmetry property, and they are related by:

$$k_{\mathrm{ij}}\left(\omega ,\mathit{x}\right)=\frac{\nu \varphi \left(\mathit{x}\right)}{-\mathrm{i\omega }}{\alpha }_{\mathrm{ij}}^{-1}\left(\omega ,\mathit{x}\right).$$(23)

The frequency-dependent dynamic permeability relation (17) between the macroscopic flow velocity 〈v〉 = ϕ〈v〉_p and the macroscopic pressure gradient, is referred to as a dynamic regime Darcy’s law. The other relation (16) is a dynamic regime Newton’s law.

To introduce the analogous notions of “dynamic thermal tortuosity and permeability”¹⁶, we now proceed in a somewhat parallel manner. For heat exchange, we are interested in the description of the excess-temperature-pattern τ in air, appearing in the above finite coarse-graining-ball of volume $V_p^{\left(b\right)}$ and pore surface $S_p^{\left(b\right)}$ centered at x. It obeys the following Fourier equation in $V_p^{\left(b\right)}$ (which is one among the complete set of governing coupled equations mentioned above (7.1-6) [9]),

$${\rho }_0{c}_P\frac{\mathrm{\partial }\tau }{\mathrm{\partial }t}={\beta }_0{T}_0\frac{\mathrm{\partial }p}{\mathrm{\partial }t}+\kappa \Delta \tau ,\hspace{1em} \mathrm{in} V_p^{\left(b\right)},$$(24)

with boundary conditions:

$$\tau =0,\hspace{1em}\mathrm{on} S_p^{\left(b\right)}.$$(25)

In the full set of equations of which (24) and (25) is an excerpt, p is a thermodynamic variable; it will now lose this status in the simplified way in which we will consider (24), isolated from the rest. Since in Biot theory we apply a 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 = 〈p〉 or equivalently 〈pv〉 = P〈v〉), which is uniformly variable (P(t, r) = (r − x) · ∇P(t, x) + P(t, x), where ∇P(t, x) is a vector constant – it has null derivative ∂/∂r), and a deviatoric part p − P, with zero mean value (in either sense, 〈(p − P)〉 = 0 or 〈(p − P)v〉 = 0).

Parallel to what was said earlier for the velocity field, we now claim that, underlying Biot’s local theoretical approach is the notion that the excess temperature pattern τ(t, r) in the region of central position x and thus its macroscopic mean 〈τ〉(t, x), which is the quantity of interest to us, can be viewed as responding to the history of macroscopic terms β₀T₀∂P/∂t applied (in the present and in the past) at this very position. To define this response, we first expand, as before, in a suitable “uniform manner” (periodic or stationary random), the domains $V_p^{\left(b\right)}$ and $S_p^{\left(b\right)}$ so that they encompass the entire space outside the sphere, as if the material were macroscopically homogeneous. Next, we assume that the excess temperature pattern we seek, is approximately the same as that (periodic or stationary random) determined in the following well-posed problem:

$${\rho }_0{c}_P\frac{\mathrm{\partial }\tau }{\mathrm{\partial }t}={\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \Delta \tau , \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(26)

$$P=f\left(t\right), \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(27)

where f(t) is a spatial constant which represents the history of the macroscopic pressure at the central position x of the original averaging ball that served to construct the abstract homogeneous unbounded medium representation ($V_p^{\left(\mathrm{ub}\right)}$ and $S_p^{\left(\mathrm{ub}\right)}$, periodic or stationary random), with boundary conditions:

$$\tau =0, \hspace{1em}\mathrm{on} S_p^{\left(\mathrm{ub}\right)}.$$(28)

We can justify the way this problem is written as follows. Since P represents now the time variations of the macroscopic pressure at the central position x of the original averaging sphere, the equation (26) is reminiscent of the Stokes equation (6) or (10) used to define the velocity pattern after incompressibility was introduced. It should be recalled that the π field did not exactly mimic the deviatoric overpressure, p − P, i.e. the 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 τ (and thus its macroscopic mean) in the original sphere, and the course of the macroscopic term β₀T₀∂P/∂t with the spatial constant P as the averaged central variable¹⁷. Finally, the problem (24) and (25) becomes well-posed if it is simplified and extended over the whole space in the way indicated to give (26)–(28).

As for the time variations of the uniform excitation P, they can be assumed to be arbitrary (as long as they are compatible with long wavelengths), but again, this does not mean that the intervening temporal dispersion will be general. For arbitrary time variations, the solution is a convolution of the elementary solution of the impulse Fourier problem considered in [18] (response to a Dirac delta heat source in time). This solution has the constraints given in Appendix C in [4], since it is constructed with purely damped normal modes (Avellaneda and Torquato’s viscous 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” $\widehat{\alpha }\mathrm{\prime }\left(\mathit{x}\right)$, associated with response, kernel factors α′(t, x), such that,

$${\rho }_0{c}_P\frac{\mathrm{\partial }\langle \tau {\rangle }_p}{\mathrm{\partial }t}\left(t\right)={\widehat{\alpha }}^{{\mathrm{\prime }}^{-1}}\left(\mathit{x}\right){\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}\left(t\right)={\int }_{-\mathrm{\infty }}^t\mathrm{}\mathrm{d}{t}^{\mathrm{\prime }}{\alpha }^{{\mathrm{\prime }}^{-1}}\left(t-t^{\mathrm{\prime }},\mathit{x}\right){\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}\left(t^{\mathrm{\prime }}\right),$$(29)

or equivalently, local macroscopic operators of “thermal permeability” $\widehat{k}\mathrm{\prime }\left(\mathit{x}\right)$, associated with response, kernel factors k′(t, x), such that,

$$\varphi \left(\mathit{x}\right)\langle \tau {\rangle }_p\left(t\right)=\frac{{\widehat{k}}^{\mathrm{\prime }}\left(\mathit{x}\right)}{\kappa }{\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}\left(t\right)=\frac{1}{\kappa }{\int }_{-\mathrm{\infty }}^t\mathrm{}\mathrm{d}t\mathrm{\prime }k\mathrm{\prime }(t-t\mathrm{\prime },\mathit{x}){\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}\left(t\mathrm{\prime }\right).$$(30)

The above relations are written with 〈 〉_p performed in the averaging-volume of the central position x where the f(t) values of the macroscopic source term β₀T₀∂P/∂t have been recorded. Recall that the values of the material parameters, porosity ϕ and kernels α′ and k′, depend on x if there are gradients in macroscopic properties; note that for shortness this argument x is omitted in the fields in the equations ((14)–(15), (29)–(30)).

In harmonic regime the operator relations (29) and (30) become multiplications with the hat quantities becoming the usual frequency-dependent functions, dynamic thermal tortuosity α′(ω, x):

$$-\mathrm{i\omega }{\rho }_0{c}_P\alpha \mathrm{\prime }(\omega ,\mathit{x})\langle \tau {\rangle }_p={\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t},$$(31)

and dynamic thermal permeability k′(ω, x):

$$\varphi \left(\mathit{x}\right)\langle \tau {\rangle }_p=\frac{k^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)}{\kappa }{\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}.$$(32)

At each position they are computed by solving in the associated unbounded (ub) homogeneous periodic or stationary random medium the following problem:

$$-\mathrm{i\omega }{\rho }_0{c}_P\tau =-\mathrm{i\omega }{\beta }_0{T}_{0}P+\kappa \Delta \tau , \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(33)

$$P=\mathrm{a} \mathrm{given} \mathrm{spatial} \mathrm{constant}, \hspace{1em}\mathrm{in} V_p^{\left(\mathrm{ub}\right)},$$(34)

with boundary conditions:

$$\tau =0,\hspace{1em}\mathrm{on} S_p^{\left(\mathrm{ub}\right)}.$$(35)

From the unique solution τ of this 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 α′(ω, x) or k′(ω, x):

$$k\mathrm{\prime }(\omega ,\mathit{x})=\frac{\kappa \varphi \left(\mathit{x}\right)}{{\beta }_0{T}_0\mathrm{\partial }P/\mathrm{\partial }t}\langle \tau {\rangle }_p, {\alpha }^{\mathrm{\prime }-1}(\omega ,\mathit{x})=\frac{\left(-\mathrm{i\omega }{\rho }_0{c}_P\right)}{{\beta }_0{T}_0\mathrm{\partial }P/\mathrm{\partial }t}\langle \tau {\rangle }_p,$$(36)

which are related by:

$$k^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)=\frac{{\nu }^{\mathrm{\prime }}\varphi \left(\mathit{x}\right)}{-\mathrm{i\omega }}{\alpha }^{\mathrm{\prime }-1}\left(\omega ,\mathit{x}\right).$$(37)

The frequency-dependent dynamic thermal permeability relation (32) between the macroscopic excess temperature 〈τ〉 = ϕ〈τ〉_p and the macroscopic source term, is referred to as a dynamic regime “thermal Darcy’s law”. The other relation (31) may be referred to as a dynamic regime “thermal Newton’s law”. It remains here to explain how these notions of thermal dynamic permeability and tortuosity determine the effective bulk modulus of the air. The link between thermal dynamic permeability/tortuosity and effective bulk modulus has been presented e.g. in [4, 9, 18]. For completeness we briefly recall here the relation.

As derived in [8, 9], the thermodynamic equation of state of the saturating fluid leads, after linearization, to the following relation between excess pressure, condensation b ≡ ρ′/ρ₀ where ρ′ is the excess density, and excess temperature:

$$\gamma {\chi }_{0}p=b+{\beta }_0\tau , \hspace{1em}\mathrm{in} V_p^{\left(b\right)},$$(38)

where the notation χ₀ is employed for the adiabatic compressibility 1/K_a (setting, τ = 0, in (38), the resulting relation, γχ₀p = b, must describe the isothermal excess pressure-condensation relation, so that, obviously, γχ₀ = 1/K₀). Taking the average and the time derivative it gives the relation, since 〈p〉_p = P,

$$\gamma {\chi }_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t}=\frac{\mathrm{\partial }\langle b{\rangle }_p}{\mathrm{\partial }t}+{\beta }_0\frac{\mathrm{\partial }\langle \tau {\rangle }_p}{\mathrm{\partial }t}.$$(39)

In harmonic regime e^−iωt and working with the complex amplitudes we can use (31) to express the right-hand side in (39). Using the general thermodynamic identity seen in Section 2 and suppressing the common factor (−iω) this gives after straightforward calculation the relation,

$${\chi }_0\left[\gamma -\frac{\gamma -1}{{\alpha }^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)}\right]P=\langle b{\rangle }_p.$$(40)

Therefore in harmonic regime we find a relation having the form,

$$P=K_f(\omega ,\mathit{x})\langle b{\rangle }_p={\chi }_f^{-1}(\omega ,\mathit{x})\langle b{\rangle }_p,$$(41)

with

$$K_f^{-1}\left(\omega ,\mathit{x}\right)={\chi }_f\left(\omega ,\mathit{x}\right)={\chi }_0\left[\gamma -\frac{\gamma -1}{{\alpha }^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)}\right],$$(42)

$$={\chi }_0\left[\gamma -\left(\gamma -1\right)\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}\varphi }{k}^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)\right].$$(43)

It gives the Fourier coefficients of the kernel functions ${\chi }_f^{-1}(t,\mathit{x})$ or K_f(t, x) which appear in the arbitrary time-variable-regime, operator or integral equation,

$$P(t,\mathit{x})={\widehat{\chi }}_f^{-1}\langle b{\rangle }_p(t,\mathit{x})={\int }_{-\mathrm{\infty }}^t\mathrm{}{\chi }_f^{-1}(t-t\mathrm{\prime },\mathit{x})\langle b{\rangle }_p(t\mathrm{\prime },\mathit{x})\mathrm{d}t\mathrm{\prime },$$(44)

$$={\widehat{K}}_f\langle b{\rangle }_p(t,\mathit{x})={\int }_{-\mathrm{\infty }}^t\mathrm{}{K}_f(t-t\mathrm{\prime },\mathit{x})\langle b{\rangle }_p(t\mathrm{\prime },\mathit{x})\mathrm{d}t\mathrm{\prime }.$$(45)

We recognize here, the relation announced in Section 1. The function χ_f(ω, x), resp. its inverse K_f(ω, x), represents a local dynamic compressibility, resp. dynamic bulk modulus, of the saturating fluid, that depends on frequency because of the thermal exchanges between fluid and solid.

Finally, we note that it is customary to define a normalized dynamic compressibility, β(ω, x), such that, in harmonic regime,

$${\chi }_0\beta \left(\omega ,\mathit{x}\right)\varphi \left(\mathit{x}\right)\frac{\mathrm{\partial }P}{\mathrm{\partial }t}={\chi }_0\beta \left(\omega ,\mathit{x}\right)\frac{\mathrm{\partial }\langle p\rangle }{\mathrm{\partial }t}=\frac{\mathrm{\partial }\langle b\rangle }{\mathrm{\partial }t}=-\nabla ·\langle \mathit{v}\rangle .$$(46)

It is given by,

$$\beta \left(\omega ,\mathit{x}\right)=\frac{{\chi }_f\left(\omega ,\mathit{x}\right)}{{\chi }_0}=\gamma -\frac{\gamma -1}{{\alpha }^{\mathrm{\prime }}\left(\omega ,\mathit{x}\right)}.$$(47)

Note that in (46) we have directly written that, ∂〈b〉/∂t = −∇ · 〈v〉, using the identity 〈∇ · v〉 = ∇ · 〈v〉, which is valid (with no indice p) even if the porosity varies spatially. We refer the reader to [9] where the explanations are given on how this identity is obtained (where in the left, ∇ = ∂/∂r, and in the right, ∇ = ∂/∂x).

4 Generalization of the method of Pontrjagin, Andronow and Witt, in the presence of radioactive decay of diffusing species

As will be seen in Section 5, there is a direct proportionality between the distribution of excess temperature field τ(r) in thermal problem defined in (33)–(35) and the mean survival time $\bar{t}\left(\mathit{r}\right)$ of a tiny particle released at position r, moving and decaying under the combined action of thermal agitation (Brownian motion) and radioactive disintegration, and instantly absorbed by reaching a specified boundary (that of the pore walls). (The overbar designates the average over realizations of the random walk in the limit of large numbers). In this section, we will first study how to compute the mean survival time $\bar{t}\left(\mathit{r}\right)$ of the particle with instant absorption (instant decay), at the pore walls. In a paper on the statistical treatment of dynamic systems, Pontrjagin et al. [17], have discussed a general method, allowing in particular to establish a differential equation for the average time it takes a Brownian motion particle moving under the combined effect of thermal agitation and a stationary field of force, to reach a defined boundary. Focusing on this particular problem, an helpful detailed account of Pontrjagin, Andronow and Witt’s general method has been presented by Lifson and Jackson in [16]. Here, closely following their exposition, we will establish a differential equation for the average time it takes for a Brownian motion particle to reach the pore walls, moving under the combined effects of thermal agitation and superimposed radioactive decay of the diffusing species. The method of calculation of $\bar{t}\left(\mathit{r}\right)$ in a complex geometry will be based on the known solutions of the differential equation in the most simple open (disk and ball in 2D and 3D, see Eqs. (73) and (74)).

In their presentation, Lifson and Jakson [16], as in the general considerations in [17], started with a closed region V bounded by a surface S. Let us substitute for volume V our pore space V_p, occupied by air, and for enclosing surface S, the pore wall surface S_p between air and porous solid. Our pore space V_p is not a closed region but this is not a significant difference. What is important is that this space is only bounded by the surface S_p. If it happens in our following reasonings that V_p and S_p are not macroscopically homogeneous regions, we will understand to identify them with the homogeneous domains we referred previously by the notations $V_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, $S_p^{\left(\mathrm{ub}\right)}\left(\mathit{x}\right)$, with the dependence (x) not indicated, for simplicity. Then, with the substitution of $V_p^{\left(\mathrm{ub}\right)}$, $S_p^{\left(\mathrm{ub}\right)}$, for V, S, in [16], only minimal change in the considerations will occur. This is why we take the liberty, in the following, to paraphrase almost completely parts of Lifson and Jakson’s exposition, mutatis mutandis only what is necessary in the presence of the superimposed radioactive decay of the diffusing particles.

As we will see, varying the radioactive decay time τ_b, from ∞ (no disintegration, low frequency limit) to 0 (instant disintegration, high frequency limit), will allow us to explore the effective bulk-modulus function in its Laplace-Transform variable representation, $K_f^{\mathcal{L}}\left(s\right)$, from s = 0 to s = ∞, where s is related to ω according to s = −iω using the convention e^−iωt. In [16], the conservative force field, F(r), acts on the diffusing species and is required to be finite everywhere in the V_p region. For the sake of generality in this section we retain it, though we need not consider it in the sequel. We will later (Sect. 5) set, F(r) = 0.

We consider Brownian motion particles, moving, in the region V_p, under the combined influence of thermal agitation and the force F; the new element compared to [16] is that each diffusing particle undergoes “radioactive decay” with time characteristic τ_b in the volume V_p, i.e., if present (in bulk) at time t the particle has probability κ_b dt = dt/τ_b of being destroyed in the infinitesimal time interval (t, t +dt).

“Disintegration” may refer to the particle itself that can be removed e.g. by chemical reactions, or to some physical excitation carried by the particle that can be deexcited by interaction in the bulk. In all cases we speak of “the particle” and of its disappearance. We are not concerned here with the physical implementation but only with the posed mathematics of the diffusion-disintegration process which is assumed to mimick Brownian motion and radioactive decay laws. We define the cumulative probability function $\mathcal{D}(\mathit{r},t)$ (replacing the Lifson and Jackson cumulative probability function $\mathcal{W}(\mathit{r},t)$), as the cumulative probability function that a particle located at point r at time t = 0 will reach the bounding surface S_p during the time interval (0, t), or else, will disappear in the bulk V_p without reaching S_p during the same interval.

The “$\mathcal{D}$” – for “$\mathcal{D}$isappearance” (or “$\mathcal{D}$ecay”, “$\mathcal{D}$isintegration”, “$\mathcal{D}$elete” or $\mathcal{D}$eath”) – qualification for this probability is appropriate if any particle touching the surface S_p, or the excitation it carries, is instantly removed; $\mathcal{D}(\mathit{r},t)$ only differs from $\mathcal{W}(\mathit{r},t)$ in [16] in that the superimposed radioactive decay (of diffusing particles or their excitation), provides an additional spontaneous means of removal (in bulk). A number of simple properties previously listed for $\mathcal{W}(\mathit{r},t)$ in [16], still hold for $\mathcal{D}(\mathit{r},t)$. Evidently, for a point s on the surface S_p, $\mathcal{D}(\mathit{s},t)=1$ for all t. For an interior point we have $\mathcal{D}(\mathit{r},\mathrm{\infty })=1$, while $\mathcal{D}(\mathit{r},0)=0$. A difference with [16] is that the condition,

$$\underset{\tau \to 0}{\lim }\left[\mathcal{W}\left(\mathit{r},\tau \right)/\tau \right]=0,$$(48)

now becomes, after accounting for the finite rate of radioactive decay,

$$\underset{\tau \to 0}{\lim }\left[\mathcal{D}\left(\mathit{r},\tau \right)/\tau \right]=\frac{1}{{\tau }_b}={\kappa }_b.$$(49)

To see this, recall that $\mathcal{D}(\mathit{r},\tau )$ is the probability that a particle (or its excitation), starting from r at t = 0, dies before t = τ, either by touching S_p, or by spontaneously decaying in the bulk. So to calculate $\mathcal{D}$, we launch N → ∞ particles at r at t = 0, and, after a time interval $\tau$, we see N″ paths of particles that have not touched S_p or died in the bulk. The remaining N′ = N − N″ paths, are dead paths, either stopped in the bulk when particles died spontaneously, or stopped at the boundary S_p where instant decay occurs. Then, $\mathcal{D}(\mathit{r},\tau )={\lim }_{N\to \mathrm{\infty }}(N\mathrm{\prime }/N)$. When τ → 0 there is not enough time to reach S_p, so the non-dead path amounts to N″ predicted on the only basis of mass “radioactive” absorption: N″ = N(1 − κ_bτ); consequently N′ = Nκ_bτ and $\mathcal{D}={\kappa }_b\tau$, as stated in (49).

We can now relate (as Lifson and Jackson did for probability $\mathcal{W}$) the probability $\mathcal{D}$ at time τ, to its value at the next time t + τ. For this purpose we define the function ι(r, t, ρ), such that, ι(r, t, ρ)dρ, is the transition probability for a particle, to go from point r, in V_p, at time t = 0, to a volume element dρ about another point, ρ, in V_p, at time t, without having touched the surface S_p during the time interval (0, t). This is not the same as υ(r, t, ρ)dρ in [16], because, as our particle now decays in a manner typical of radioactive decay, ι is therefore reduced compared to its counterpart υ without radioactive decay:

$$\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)=e^{-\tau /{\tau }_b}\upsilon \left(\mathit{r},\tau ,\mathit{\rho }\right).$$(50)

We may call $\mathcal{L}(\mathit{r},t)\equiv {\int }_{V_p}\mathrm{}\iota (\mathit{r},t,\mathit{\rho })\mathrm{d}\mathit{\rho }$ the “$\mathcal{L}$ife expectancy” from going to r, at t = 0, to anywhere else in the pore space, at time t. As a particle must be counted as either deceased or living, the function ι(r, t, ρ) is normalized in the sense that,

$$\mathcal{L}\left(\mathit{r},t\right)\equiv {\int }_{V_p}\mathrm{}\iota \left(\mathit{r},t,\mathit{\rho }\right)d\mathit{\rho }=1-\mathcal{D}\left(\mathit{r},t\right).$$(51)

From (51) and (49) we also note that,

$$\underset{\tau \to 0}{\lim }{\int }_{V_p}\mathrm{}\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)d\mathit{\rho }=1-{\kappa }_b\tau +\dots$$(52)

It is now clear that the following integral equation holds:

$$\mathcal{D}\left(\mathit{r},t+\tau \right)=\mathcal{D}\left(\mathit{r},\tau \right)+{\int }_{V_p}\mathrm{}\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)\mathcal{D}\left(\mathit{\rho },t\right)\mathrm{d}\mathit{\rho }.$$(53)

The meaning of the above equation is as follows: If a particle has reached the surface S_p during the time interval (0, t + τ) or has died in the volume V_p during this time interval without touching it, then, either it has reached S_p or died in V_p during the time interval (0, τ), or else, it gets to the point ρ at time τ without touching the surface S_p or decaying in V_p, and then, reaches from this point the surface S_p or decays in V_p without touching it, in the time interval (τ, τ + t)¹⁸.

We now recall known results for the successive moments of the ordinary transition probability υ₀(r, t, ρ) in an infinite medium, which will immediately translate into identical results for the successive moments of ι(r, t, ρ). For Brownian motion without disintegration in an infinite medium in the presence of a field of force, the particle has an average velocity at a point, proportional to the force, i.e., the following limit exists:

$$\underset{\tau \to 0}{\lim }\frac{\langle {\rho }_i-r_i\rangle }{\tau }=\frac{F_i}{f}=\underset{\tau \to 0}{\lim } {\tau }^{-1}{\int }_{V_p}\mathrm{}\left({\rho }_i-r_i\right){\upsilon }_0\left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho },$$(54)

where f is the hydrodynamic friction coefficient and υ₀(r, τ, ρ) the transition probability in the infinite medium. Furthermore, the diffusion coefficient, D = kT/f, is related to the second moment of the displacement, according to,

$$\underset{\tau \to 0}{\lim }\frac{\langle \left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right)\rangle }{\tau }=2D{\delta }_{\mathrm{ij}},$$(55)

$$=\underset{\tau \to 0}{\lim } {\tau }^{-1}{\int }_{V_p}\mathrm{}\left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right){\upsilon }_0\left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho }.$$(56)

To see this we note that the effect of the, force-driven, constant drift velocity F_i/f, which brings in (ρ_i − r_i) a term τF_i/f which superposes to a random, non-driven one, will bring a cross term, τ⁻¹(F_iF_j)τ²/f², which does not contribute in the limit τ → 0, and two non-cross-terms, that disappear by the absence of direction of the random non-driven motion. Thus, no force-driven effect appears in the second moment, and the transition probability υ₀(r, t, ρ) can be replaced by the one without force, Gaussian in each direction, $=e^{-|\mathbf{\rho }-\mathbf{r}{|}^2/2\mathrm{dDt}}/\left(2\pi \mathrm{dDt}\right)$, where d is the dimension. The result then expresses Einstein’s translational Brownian motion formula, $\overline{\left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right)}=2\mathrm{D\tau }{\delta }_{\mathrm{ij}}$. For all higher moments it can be checked that,

$$\underset{\tau \to 0}{\lim }{\tau }^{-1}\langle {\left({\rho }_i-r_i\right)}^l{\left({\rho }_j-r_j\right)}^m{\left({\rho }_k-r_k\right)}^n\rangle =0,$$(57)

$$=\underset{\tau \to 0}{\lim }{\tau }^{-1}{\int }_{V_p}\mathrm{}{\left({\rho }_i-r_i\right)}^l{\left({\rho }_j-r_j\right)}^m{\left({\rho }_k-r_k\right)}^n{\upsilon }_0\left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho },$$(58)

$$\left(l+m+n\right)>2.$$(59)

Now, in view of the fact that the function ι(r, t, ρ) must approach the ordinary (i.e. infinite medium) transition probability as τ approaches to zero, the moments of ι(r, t, ρ) must approach the moments of the ordinary transition probability. We then have, from equations (54), (55), and (57),

$$\underset{\tau \to 0}{\lim }\frac{\langle {\rho }_i-r_i\rangle }{\tau }=\frac{F_i}{f}=\underset{\tau \to 0}{\lim }{\tau }^{-1}{\int }_{V_p}\mathrm{}\left({\rho }_i-r_i\right)\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho },$$(60)

$$2D{\delta }_{\mathrm{ij}}=\underset{\tau \to 0}{\lim }{\tau }^{-1}{\int }_{V_p}\mathrm{}\left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right)\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho },$$(61)

$$0=\underset{\tau \to 0}{\lim }{\tau }^{-1}{\int }_{V_p}\mathrm{}{\left({\rho }_i-r_i\right)}^l{\left({\rho }_j-r_j\right)}^m{\left({\rho }_k-r_k\right)}^n\iota \left(\mathit{r},\tau ,\mathit{\rho }\right)\mathrm{d}\mathit{\rho },\left(l+m+n\right)>2.$$(62)

We now use the integral equation for $\mathcal{D}(\mathit{r},t)$ (Eq. (53)), to derive a partial differential equation. Expanding $\mathcal{D}(\mathit{\rho },t)$ about the point r we have

$$\frac{\mathrm{\partial }\mathcal{D}\left(\mathit{r},t\right)}{\mathrm{\partial }t}=\underset{\tau \to 0}{\lim }{\tau }^{-1}\left[\mathcal{D}\left(\mathit{r},t+\tau \right)-\mathcal{D}\left(\mathit{r},t\right)\right],$$(63)

$$=\underset{\tau \to 0}{\lim }{\tau }^{-1}\{-\mathcal{D}(\mathit{r},t)+\mathcal{D}(\mathit{r},\tau )+{\int }_{V_p}\mathrm{}d\mathit{\rho }\iota (\mathit{r},\tau ,\mathit{\rho }) [\mathcal{D}(\mathit{r},t)+\sum _i\mathrm{}({\rho }_i-r_i)\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }{r}_i}+\frac{1}{2!}\sum _{i,j}\mathrm{}\left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right)\frac{{\mathrm{\partial }}^2\mathcal{D}}{\mathrm{\partial }{r}_i\mathrm{\partial }{r}_j}+\frac{1}{3!}\sum _{i,j,k}\mathrm{}\left({\rho }_i-r_i\right)\left({\rho }_j-r_j\right)\times ({\rho }_k-r_k)\frac{{\mathrm{\partial }}^3\mathcal{D}}{\mathrm{\partial }{r}_i\mathrm{\partial }{r}_j\mathrm{\partial }{r}_k}+\dots \left]\right\}.$$(64)

Utilizing the various limiting properties of ι and $\mathcal{D}$ presented in (49), (52), (60), (61), (62), we get a differential equation for the cumulative probability $\mathcal{D}(\mathit{r},t)$ that a particle launched at r in V_p at t = 0, vanishes before t, either by touching S_p or decaying in the bulk:

$$\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}=f^{-1}\mathit{F}·\nabla \mathcal{D}+D\Delta \mathcal{D}-{\kappa }_b\mathcal{D}+{\kappa }_b.$$(65)

This interprets as the result obtained for $\mathcal{W}$ in Lifson and Jackson [16] equation (9), modified by the presence of a term $-{\kappa }_b\mathcal{L}$ added in the left-hand side to account for an additional rate of change coming from the radioactive disintegration.

Now observe that the equation, $\mathcal{D}(\mathit{r},t+\mathrm{d}t)=\mathcal{D}(\mathit{r},t)+\mathrm{d}t\mathrm{\partial }\mathcal{D}(\mathit{r},t)/\mathrm{\partial }t,$ should be interpreted as follows: If a particle starting from r at t = 0 vanishes during the time interval (0, t + dt) (either by reaching S_p or decaying in the bulk V_p), then, either it vanishes during the time interval (0, t) (by reaching S_p or decaying in the bulk V_p), or it is still alive at t, and then, during the time interval (t, t + dt), reaches for the first time S_p or vanishes in the bulk. This shows that $\mathrm{d}t\mathrm{\partial }\mathcal{D}(\mathit{r},t)/\mathrm{\partial }t$ is the probability of a particle to vanish (reach for the first time S_p or vanish in the bulk) during (t, t + dt), when starting from r at time t = 0. Therefore, the average survival time for a particle is:

$$\bar{t}\left(\mathit{r}\right)={\int }_0^{\mathrm{\infty }}\mathrm{}t\frac{\mathrm{\partial }\mathcal{D}}{\mathrm{\partial }t}\left(\mathit{r},t\right)\mathrm{d}t.$$(66)

From this we finally obtain the following desired differential equation satisfied by $\bar{t}\left(\mathit{r}\right)$:

$$\left(1/f\right)\mathit{F}·\nabla \bar{t}+D\Delta \bar{t}-{\kappa }_b\bar{t}=-1, \hspace{1em}\mathrm{in} V_p,$$(67)

with boundary condition:

$$\bar{t}=0,\hspace{1em}\mathrm{on} S_p.$$(68)

When the field F is zero one may directly write the solutions for a line, cylinder, and sphere. One has, with notation ${\mu }_a\equiv a\sqrt{{\kappa }_b/D}$,

One dimensional case: length of line 2L, x ∈ [−L, L],

$$\bar{t}\left(x\right)={\tau }_b\left[1-\frac{\cosh \left({\mu }_x\right)}{\cosh \left({\mu }_L\right)}\right];$$(69)

Two-dimensional case: radius of cylinder R, r ∈ [0, R],

$$\bar{t}\left(r\right)={\tau }_b\left[1-\frac{I_0\left({\mu }_r\right)}{I_0\left({\mu }_R\right)}\right];$$(70)

Three-dimensional case: radius of sphere R, r ∈ [0, R],

$$\bar{t}\left(r\right)={\tau }_b\left[1-\frac{{\mu }_R\sinh ({\mu }_r)}{{\mu }_r\sinh \left({\mu }_R\right)}\right].$$(71)

5 Efficient algorithm to compute $K_f\left(\omega \right)$ and other functions and intervening parameters

With F = 0 we can directly compare the average survival time problem (67)–(68) with the excess temperature pattern problem (33)–(34). Examination of both shows that the following quantities, with the same dimensions, can be put into correspondence:

$$\bar{t}\leftrightarrow \frac{1}{-\mathrm{i\omega }\frac{{\beta }_0{T}_0}{{\rho }_0{c}_P}P}\tau ;$$(72)

$$D\leftrightarrow \frac{\kappa }{{\rho }_0{c}_P}\equiv {\nu }^{\text{'}};$$(73)

$${\kappa }_b\equiv \frac{1}{{\tau }_b}\leftrightarrow -\mathrm{i\omega }.$$(74)

Therefore, when we have solved the mean survival time problem, we have solved the thermal problem. For example we can compute the thermal dynamic permeability by, first, writing the relation (e.g. combine (32) with the correspondence (72)–(74)):

$$k^{\mathrm{\prime }}\left(\omega \right)=\varphi D\langle \bar{t}\rangle ,$$(75)

and next, substituting in the calculations that are made to evaluate the product $D\langle \bar{t}\rangle$, the values D = ν′ and κ_b = −iω of the coefficients D and κ_b.

Likewise, we can compute the bulk modulus K_f(ω) by, first, writing the relation (e.g. combine (43) with (75) and the correspondence (72)–(74)):

$$K_f^{-1}\left(\omega \right)=K_a^{-1}\left[\gamma -\left(\gamma -1\right){\kappa }_b\langle \bar{t}\rangle \right],$$(76)

and next, substituting the above values of D and κ_b in the expressions obtained for the product ${\kappa }_b\langle \bar{t}\rangle$. In fact, it means that in the Laplace-variable representation of the kernel K_f(t) of the bulk modulus operator ${\widehat{K}}_f$, (45)¹⁹, we simply have, ${\left(K_f^{\mathcal{L}}\right)}^{-1}\left(s\right)=K_a^{-1}\left[\gamma -(\gamma -1){\kappa }_b\langle \bar{t}\rangle \right]$, with this time, D = ν′ and κ_b = s: The problem of mean survival time, most naturally translates in the description of thermal response factors expressed in Laplace representation.

Having established the connection between the mean survival time problem and the thermal problem, let us now focus on the former and show, how a simple geometric construction of the Brownian random motion path considered in Torquato and Kim [13], see Figure 1, leads us to a simple geometric algorithm to calculate the mean survival time $\bar{t}$ and its average $\langle \bar{t}\rangle$ 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 $\bar{t}\left(\mathit{r}\right)$ of random walkers released in one configuration at a given place, or to calculate its average $\langle \bar{t}\rangle$ either viewed as a volume average (in a given configuration) when walkers are released at any place with equal probability (Lorentz viewpoint), or as an ensemble average when they are released at the same position but the configuration is varied appropriately (Gibbs viewpoint), is now described.

Let r₁ 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₁ which does not overlap any solid trap. As the random motion establishes no preference on the directions, the next position r₁ of the walker is chosen at random on this circle/sphere. In this way, Brownian motion paths, r₁, r₂, …, 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 (${\bar{t}}_R^{\left(s\right)}$, s for surface), to reach R, will be given by the inverse of Einstein’s relation, i.e.²⁰,

$${\bar{t}}_R^{\left(s\right)}=\frac{R^2}{4\mathrm{D}}, \left(2\mathrm{D}\right); {\bar{t}}_R^{\left(s\right)}=\frac{R^2}{6\mathrm{D}}, \left(3\mathrm{D}\right),$$(77)

(ii) the decay occurs before radius R is reached, the associated probability is 1 − p(R) and we assume that these aborted random walks are performed in an average time ${\bar{t}}_R^{\left(b\right)}$ (b for bulk). Note that ${\bar{t}}_{\mathrm{\infty }}^{\left(b\right)}$ is the average survival time τ_b of the walker in the infinite medium.

Let ${\bar{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:

$${\bar{t}}_R={\tau }_b\left[1-\frac{1}{I_0\left({\mu }_R\right)}\right], \left(2\mathrm{D}\right);\hspace{1em}{\bar{t}}_R={\tau }_b\left[1-\frac{{\mu }_R}{\sinh \left({\mu }_R\right)}\right], \left(3\mathrm{D}\right).$$(78)

The above two types of outcomes allow us to write it, also, in the following form:

$${\bar{t}}_R=\left[1-p\left(R\right)\right]{\bar{t}}_R^{\left(b\right)}+p\left(R\right){\bar{t}}_R^{\left(s\right)}.$$(79)

To get another relationship between the same two, but yet unknown quantities, p(R) and ${\bar{t}}_R^{\left(b\right)}$, let us now imagine that the boundary at radius R is fictitious and drawn in an infinite medium. Let us again follow the particles released in the center O. We know that on average they live the time τ_b. A proportion 1 − p(R) of them do not reach radius R and live the average time ${\bar{t}}_R^{\left(b\right)}$. The remaining proportion p(R) reach R and then will later disappear at some instant. The latter particles therefore live first the average time ${\bar{t}}_R^{\left(s\right)}$ to reach R, and then take the additional averagetime τ_b to disappear. Thus, we have the obvious relation,

$${\tau }_b=\left[1-p\left(R\right)\right]{\bar{t}}_R^{\left(b\right)}+p\left(R\right)\left[{\bar{t}}_R^{\left(s\right)}+{\tau }_b\right].$$(80)

Combined with the preceding, it reads, ${\tau }_b={\bar{t}}_R+p\left(R\right){\tau }_b$ (irrespective of the dimension 2 or 3), that is:

$${\bar{t}}_R={\tau }_b\left[1-p\left(R\right)\right].$$(81)

Inserting in (81) the mean survival time expressions (78) we therefore find for the hitting probability p(R) the following expressions in the 2D and 3D cases:

$$p\left(R\right)=\frac{1}{I_0\left({\mu }_R\right)}, \left(2\mathrm{D}\right);\hspace{1em}p\left(R\right)=\frac{{\mu }_R}{\sinh \left({\mu }_R\right)}, \left(3\mathrm{D}\right);\hspace{1em}{\mu }_R\equiv R\sqrt{\frac{{\kappa }_b}{D}}.$$(82)

Having determined the probability p(R) (for having survived till crossing for the first time the radius r = R, after having started at the center r = 0) and on the other hand the mean time ${\bar{t}}_R$ (for crossing for the first time the radius r = R or decaying before touching it, after having started at r = 0) we now can write the following simple expression:

$$\begin{array}{c}\bar{t}\left({\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_n\right)={\bar{t}}_{R_1}+p\left(R_1\right){\bar{t}}_{R_2}+p\left(R_1\right)p\left(R_2\right){\bar{t}}_{R_3}+\dots \\ \dots +p\left(R_1\right)\dots p\left(R_{n-1}\right){\bar{t}}_{R_n},\end{array}$$(83)

for the mean time for trapping, associated with the path construction, r₁, r₂, …, 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₁ at t = 0. Finally, by inserting in (83) the expression (81) of times ${\bar{t}}_{R_i}$, a systematic two by two cancellation of all intermediate terms occurs, leaving only the two ends, and giving a result in simpler form:

$$\bar{t}\left({\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_n\right)={\tau }_b\left[1-p\left(R_1\right)p\left(R_2\right)\dots p\left(R_n\right)\right].$$(84)

The exact result without truncation will be with ${\prod }_{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)$, for the probability product.

In a given realization of the medium, the final mean survival time ${\bar{t}}_{{\mathit{r}}_1}$ at position r₁, will be obtained by repeating the path constructions and taking the average of the above mean trapping times, over all these constructions:

$${\bar{t}}_{{\mathbf{r}}_1}=\underset{N\to \mathrm{\infty }}{\lim }\frac{\sum \mathrm{}\bar{t}\left({\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_n\right)}{N}.$$(85)

The sum corresponds to the systematic repetition of the constructions and N is the number of times the constructions are repeated. The r₂, …, r_n including the final index value n vary from one construction to another. Likewise in the vectors [R₁, R₂, …, R_n] specifying one construction, only the first value R₁ repeats itself, all the others, R₂, …, 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:

$${\bar{t}}_{{\mathit{r}}_1}={\tau }_b\left[1-\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}\right].$$(86)

The microscopic excess temperature field at position r₁,

$$\tau \left({\mathit{r}}_1\right)=-\mathrm{i\omega }\left({\beta }_0{T}_0/{\rho }_0{c}_P\right)P{\bar{t}}_{{\mathit{r}}_1},$$(87)

(see Eq. (72)), will be directly obtained geometrically from the construction of N → ∞ vectors [R₁, R₂, …, 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₁ and making the N different constructions for the given same sample of material. In the limit N → ∞ the mean survival time or equivalently the microscopic field τ at this position is determined.

To determine the dynamic thermal permeability we only need to have a macroscopic average 〈τ〉(x). To have it we can for example proceed as follows.

First consider Lorentz’s conception of the average. We work with one sample of the material and must consider the result of an average in a representative volume. To calculate this volume average 〈τ〉(x), where x is the cental position of the representative volume, we can randomly select an initial position r₁, accept it and increment the number of trials N if it lies in the fluid part of this volume, then perform the construction, [R₁, R₂, …, 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₁ increases, without pre-determination of the microscopic field. Of course, as the number of walkers increases and the result converges (in $1/\sqrt{N}$), the microscopic field is also automatically extracted at any point with increasing accuracy; but this is an output, not an input, in the averaging process.

With a Gibbs’ conception of the average, the initial position x is fixed and we can randomly select, each time, a new configuration of the traps. We accept it and increment the number of trials N if x lies in the fluid, each time giving rise to one single new construction [R₁, R₂, …, 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/\sqrt{N}$, with N the used number of constructions. Before moving on to validation and method examples, it is useful now to restate, further detail, and summarize our findings.

Let us return to the macroscopic average mean life time of diffusing radioactive species that are directly absorbed in the pore wall. We can consider it in the Lorentz conception when it is calculated in a single sample by means of constructions derived from random points of a representative fluid volume, or in the Gibbs conception when the point of origin is invariable and the construction each time is carried out in a new configuration including certain random variations. (Also, mixed conceptions that make partial use of both ways of averaging are possible). By indicating symbolically in all cases, with an overline (which must be interpreted appropriately), the exact way in which the average is performed, one can always designate the mean survival time in question $\bar{t}$ as used in equation (86)²¹:

$$\bar{t}={\tau }_b\left[1-\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}\right].$$(88)

The important thing is that whatever conception is considered, for the overline it will correspond to a certain vision of the macroscopic medium and a definite procedure to generate the construction [R₁, R₂, …, R_n,…] and calculate the average $\overline{{\prod }_{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}$.

Now recall the correspondence (72)–(74). There follows immediately from it and the equations (88) and (75) that the dynamic thermal permeability is given by:

$$k^{\mathrm{\prime }}\left(\omega \right)=\frac{{\nu }^{\mathrm{\prime }}\varphi }{-\mathrm{i\omega }}\left[1-\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}\right].$$(89)

In general, whatever the geometry tortuous or not, the high-frequency limit ${\alpha }_{\mathrm{\infty }}^{\text{'}}$ of α′(ω) is automatically equal to 1. It is customary to express the dynamic thermal tortuosity α′(ω) in terms of its high-frequency limit, 1, and a viscous relaxation-function, χ′(ω), taking a value of 1 when the thermal effects are in relaxed state (isothermal motions originating from the pore walls having time to fully penetrate the fluid during the cycle time), and a value of 0 when they are in frozen state (isothermal motions having no time to penetrate the fluid during the cycle time):

$$\frac{1}{{\alpha }^{\mathrm{\prime }}\left(\omega \right)}=1-{\chi }^{\mathrm{\prime }}\left(\omega \right).$$(90)

A comparison of this definition (90) with (89) and (37), immediately shows us that the averaged infinite product, $\overline{{\prod }_{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}$, is nothing but this thermal relaxation function:

$${\chi }^{\mathrm{\prime }}\left(\omega \right)=\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}.$$(91)

Once the relaxation function is conveniently directly determined through the average (91), the other thermal functions derive immediately. The effective bulk modulus, in particular, is given by

$$K_f^{-1}\left(\omega \right)=K_a^{-1}\left[1+\left(\gamma -1\right){\chi }^{\mathrm{\prime }}\left(\omega \right)\right].$$(92)

In the low frequency limit, it is a known result [9, 41], that α′(ω) can be expanded in the following Laurent series:

$${\alpha }^{\mathrm{\prime }}\left(\omega \right)=\frac{{\nu }^{\mathrm{\prime }}\varphi }{-\mathrm{i\omega }{k}_0^{\mathrm{\prime }}}+{\alpha }_0^{\mathrm{\prime }}+{{\alpha }^{\mathrm{\prime }}}_1^2\frac{\mathrm{i\omega }{k}_0^{\mathrm{\prime }}}{{\nu }^{\mathrm{\prime }}\varphi }+{{\alpha }^{\mathrm{\prime }}}_2^3{\left(\frac{\mathrm{i\omega }{k}_0^{\mathrm{\prime }}}{{\nu }^{\mathrm{\prime }}\varphi }\right)}^2+\dots ,$$(93)

where $k_0^{\mathrm{\prime }}$ is the static thermal permeability or inverse trapping constant, ${\alpha }_0^{\mathrm{\prime }}$ is the static thermal tortuosity, and ${\alpha }_1^{\mathrm{\prime }}$, ${\alpha }_2^{\mathrm{\prime }}$, etc., are some higher order purely geometrical dimensionless parameters. Comparing this expansion with the expansions that follow from using equations (91) and (82) with ${\mu }_R=R\sqrt{-\mathrm{i\omega }/\nu \mathrm{\prime }}$, as given by (73) and (74), and the known small argument expansion of the function I₀ or sinh, it is easy to derive explicit geometric expressions for the above geometric parameters. For example in 2D, we find for $k_0^{\mathrm{\prime }}$ and ${\alpha }_0^{\mathrm{\prime }}$ the following purely geometric expressions:

$$k_0^{\mathrm{\prime }}/\varphi =\frac{1}{4}\overline{\sum _{i=1}^{\mathrm{\infty }}\mathrm{}{R}_i^2},$$(94)

and,

$$\frac{{k_0^{\mathrm{\prime }}}^2{\alpha }_0^{\mathrm{\prime }}}{{\varphi }^2}=\overline{\frac{3}{64}\sum _{i=1}^{\mathrm{\infty }}\mathrm{}{R}_i^4+\frac{1}{16}\sum _{i=1}^{\mathrm{\infty }}\mathrm{}\sum _{j=i}^{\mathrm{\infty }}\mathrm{}{R}_i^2{R}_j^2}.$$(95)

We let the reader derive higher-order expressions for higher-order 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],

$${\chi }^{\mathrm{\prime }}\left(\omega \right)={\left(\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}+\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\frac{1}{{\mathrm{\Sigma }}^{\mathrm{\prime }}}+{\left(\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\right)}^{3/2}\frac{1}{4V^{\mathrm{\prime }}}+\dots ,$$(96)

with thermal characteristic length Λ′, given by, 2/Λ′ = S_p/V_p, where S_p and V_p are pore wall surface and pore volume. (The numerical factors have been chosen so that for cylindrical circular pores of radius R, one finds, Λ′ = R, Σ′ = R² and V′ = R³). To our knowledge there are no known general expressions for the higher order parameters, Σ′, V′, etc., which have surface, volume, etc., dimensions. In this limit the expression (91) cannot be expanded because for the first large radii in the product we are in the large argument limit, but for the last very small ones near the pore wall we are in the small argument limit. Nevertheless, given the HF data of the function χ′(ω) calculated using (91), we will have the means to estimate, sequentially, the values of the high frequency parameters, considering the ω → ∞ limits, e.g. for Λ′, Σ′ and V′:

$$\frac{2}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}=\underset{\omega \to \mathrm{\infty }}{\lim }{\left(\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}\right)}^{1/2}\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)},$$(97)

$$\frac{1}{{\mathrm{\Sigma }}^{\mathrm{\prime }}}=\underset{\omega \to \mathrm{\infty }}{\lim }\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}\left[\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}-{\left(\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}\right],$$(98)

and,

$$\frac{1}{4V^{\mathrm{\prime }}}=\underset{\omega \to \mathrm{\infty }}{\lim }{\left(\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}\right)}^{3/2}\left[\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}-{\left(\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\right)}^{1/2}\frac{2}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}-\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\frac{1}{{\mathrm{\Sigma }}^{\mathrm{\prime }}}\right],$$(99)

etc., where the frequency-dependent complete expressions of the p(R_i) are left. Provided the statistics is sufficient to accurately estimate the averaged infinite product $\overline{{\prod }_{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}$ the above high-frequency parameters, Λ′, Σ′, V′, etc., will be directly computable by the simulations from estimation of the above recursive limits.

6 Validation and examples of computation

We can now first validate our algorithm (Sect. 5) by considering a simple case where we have a known analytical solution within the Biot framework of the simplifications made. In particular, and for the special case where the material has 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²² 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):

$${\rho }_0\frac{\partial v}{\partial t}=-\frac{\partial P}{\partial z}+\eta \left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)v,\hspace{1em}\left(\mathrm{cross}‐\mathrm{section}\right),$$(100)

$$v=\mathbf{0},\hspace{1em}\left(\mathrm{pore}‐\mathrm{wall}\right),$$(101)

where ∂P/∂z is a given constant over the cross-section,

(ii) equations determining the excess temperature pattern τ(x, y) = τ(r):

$${\rho }_0{c}_P\frac{\partial \tau }{\partial t}={\beta }_0{T}_0\frac{\partial P}{\partial t}+\kappa (\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2})\tau , \hspace{1em}\left(\mathrm{cross}‐\mathrm{section}\right),$$(102)

$$\tau =0,\hspace{1em}\left(\mathrm{pore}‐\mathrm{wall}\right),$$(103)

where β₀T₀∂P/∂t is again a spatial constant over the cross-section, which will be combined with the following additional relations:

(iii) equation of continuity:

$$\frac{\partial b}{\partial t}+\frac{\partial v}{\partial z}=\frac{1}{{\rho }_0}\frac{\partial \rho \text{'}}{\partial t}+\frac{\partial v}{\partial z}=0, \hspace{1em}\left(\mathrm{cross}‐\mathrm{section}\right),$$(104)

and,(iv) thermodynamic equation of state:

$$\gamma {\chi }_{0}P=b+{\beta }_0\tau ,\hspace{1em}\left(\mathrm{cross}‐\mathrm{section}\right).$$(105)

In harmonic regime e^−iωt the solution fields v(r) and τ(r) are of the following form:

$$v=\frac{\nu }{-\mathrm{i\omega }}\left[-\frac{\mathrm{\partial }P}{\eta \mathrm{\partial }z}\right]\left[1-\frac{I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{\nu }}r\right)}{I_0(\sqrt{\frac{-\mathrm{i\omega }}{\nu }}R})\right];\hspace{1em}\tau =\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\left[{\beta }_0{T}_0\frac{\mathrm{\partial }P}{\kappa \mathrm{\partial }t}\right]\left[1-\frac{I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}r\right)}{I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)}\right].$$(106)

For the case of other 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 〈 〉 denoting the average at the cross-section of the tube, we can define macroscopic variables such as average axial velocity 〈v〉, average excess temperature 〈τ〉, average over-density 〈ρ′〉, or average condensation B = 〈b〉 = 〈ρ′〉/ρ₀.

Therefore, comparing the above simplified governing equations (possibly written in the general cylindrical case) with those discussed in Section 3, we see that in the harmonic regime we have, at the macroscopic, averaged level, a viscous, dynamic Darcy law, stating that:

$$\langle v\rangle =-\frac{k\left(\omega \right)}{\eta }\frac{\mathrm{\partial }P}{\mathrm{\partial }z},$$(107)

where k(ω) can be directly computed by using this definition and the average at the 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:

$$k\left(\omega \right)=\frac{\nu }{-\mathrm{i\omega }}\left[1-2\frac{I_1\left(\sqrt{\frac{-\mathrm{i\omega }}{\nu }}R\right)}{\left(\sqrt{\frac{-\mathrm{i\omega }}{\nu }}R\right)I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{\nu }}R\right)}\right].$$(108)

The bracket is just the inverse dynamic viscous tortuosity function 1/α(ω). In general, when the geometry is not trivial like here in cylindrical tubes, but tortuous, it is convenient to express the dynamic viscous tortuosity in terms of its high-frequency nontrivial limit, α_∞, and a viscous relaxation-function, χ(ω), taking a value of 1 when the viscous effects are in relaxed state (vortical motions originating from the pore walls having time to fully penetrate the fluid during the cycle time), and a value of 0 when they are in frozen state (vortical motions having no time to penetrate the fluid during the cycle time):

$$\frac{1}{\alpha \left(\omega \right)}=\frac{1}{{\alpha }_{\mathrm{\infty }}}\left[1-\chi \left(\omega \right)\right].$$(109)

For all cylindrical tubes (whatever the form of the perimeter including fractal) α_∞ is trivially equal to 1. Here, for cylindrical circular tubes of fixed radius, we see, comparing (108) and (109), that the complex expression in the brackets (108) is just the above relaxation-function χ(ω) for this cylindrical circular case.

In the same way, the comparison of the above simplified governing equations with those of Section 3, also leads us to the analogous thermal and dynamic “Darcy” law, which states that:

$$\langle \tau \rangle =\frac{k^{\mathrm{\prime }}\left(\omega \right)}{\kappa }{\beta }_0{T}_0\frac{\mathrm{\partial }P}{\mathrm{\partial }t},$$(110)

where k′(ω) can be directly computed by using this definition and the average at the cross-section, of the harmonic regime solution τ(x, y). For cylindrical circular tubes of fixed radius R, τ(x, y) = τ(r) is given by (106), and the integration shows that:

$$k^{\mathrm{\prime }}\left(\omega \right)=\frac{{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }}\left[1-2\frac{I_1\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)}{\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)}\right].$$(111)

The bracket is just the inverse of the so-called dynamic thermal tortuosity function 1/α′(ω), for this cylindrical circular case. Comparing (111) and (90) we see that the complex expression in the brackets (111) is just the thermal relaxation function for the case of cylindrical circular tubes of radius R:

$${\chi }_{\mathrm{circ}}^{\mathrm{\prime }}\left(\omega \right)=2\frac{I_1\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)}{\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)I_0\left(\sqrt{\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}}R\right)}.$$(112)

Note that for the special geometry of aligned cylindrical tubes having any type of 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, χ′(ω) = χ(ωPr), or, k(ω) = k′(ω/Pr), or α(ω) = α′(ω/Pr), where Pr ≡ ν/ν′ ≡ ηc_P/κ is the Prandtl number. In arbitrary geometries, on the contrary, there is no general relationship between the viscous and thermal response functions.

Recalling the general identification (91) of the thermal relaxation function, we see that, to validate our simulation technique in the case of a cylindrical circular tube, we only need to check the following. When we randomly launch the initial position of a random walker inside a disk of radius R, then construct and average the probability product over a large number of trials, it should be observed that:

$$\overline{\prod _{i=1}^{\mathrm{\infty }}\mathrm{}p\left(R_i\right)}={\chi }_{\mathrm{circ}}^{\mathrm{\prime }}\left(\omega \right).$$(113)

In Figure 3 we present an Argand plot of this relaxation function, as given by the theoretical closed form solution (112), and as estimated with N = 10⁵ walkers by our random motion simulation technique. For our programming in Matlab here, Lorentz volume averaging was used, where the position of the circular section is considered fixed, while the initial position of the walker is uniformly random within it. The consecutive absolute positions r_i are dummy except for their differences, that define the construction steps R_i = |r_i+1 − r_i|. To ensure uniform filling of the section required for volume averaging, it is sufficient to define the first radius R₁ 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²³, spanning the low, intermediate, and high frequencies. As mentioned, a convergence in $1/\sqrt{N}$ is slow. With N = 10⁵ 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⁷ walkers not shown, the eye no longer perceives the differences. Indeed, with N = 10⁷ walkers, the equations (93) and (95) were found to give respectively the value, 0.1249996R², for the parameter $k_0^{\mathrm{\prime }}=k_0$, whose theoretical value is, R²/8 = 0.125R², and the value, 1.3336, for the parameter ${\alpha }_0^{\mathrm{\prime }}={\alpha }_0$, whose theoretical value is, 4/3 = 1.33333….

Also plotted in the Figure 3, are the values (for the present case) of the following two general models applicable in arbitrary geometries.

Model 1, was derived by Lafarge [4, 18], in clarification of previous work by Allard and Champoux [20]. It is the thermal pure counterpart of Johnson et al. general model [28] of the viscous dynamic tortuosity-permeability. For the effective bulk modulus it reads:

$$\frac{K_a}{K_f\left(\omega \right)}=\gamma -\left(\gamma -1\right){\left[1+\frac{{\nu }^{\mathrm{\prime }}\varphi }{-\mathrm{i\omega }{k}_0^{\mathrm{\prime }}}{\left(1-\frac{4\mathrm{i\omega }{k^{\mathrm{\prime }}}_0^2}{{\nu }^{\mathrm{\prime }}{{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2{\varphi }^2}\right)}^{1/2}\right]}^{-1}.$$(114)

The thermal permeability $k_0^{\mathrm{\prime }}$ and thermal characteristic length Λ′ are independent parameters. However, the following dimensionless ratio

$$M^{\mathrm{\prime }}=\frac{8k_0^{\mathrm{\prime }}}{\varphi {{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2},$$(115)

is generally found to depart from unity by less than one order of magnitude. For the circular tubes here, M is exactly equal to 1, as it is known results that, $k_0^{\mathrm{\prime }}=R^2/8\varphi$, and, Λ′ = R. The expression (114) simplifies and happens to coincide with the so-called Allard-Champoux simpler model [20]:

$$\frac{K_a}{K_f\left(\omega \right)}=\gamma -\left(\gamma -1\right){\left[1+\frac{8{\nu }^{\mathrm{\prime }}}{-\mathrm{i\omega }{{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2}{\left(1-\frac{\mathrm{i\omega }{{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2}{16{\nu }^{\mathrm{\prime }}}\right)}^{1/2}\right]}^{-1}.$$(116)

Thus here, Model 1 is not different from the usual widely used Allard-Champoux model that considers the approximation, $k_0^{\mathrm{\prime }}=\varphi {{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2/8$. Making apparent the ratio M′ in (114) it can be seen that the magnitude of $k_0^{\mathrm{\prime }}$, or ${{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2$, determines a rollover frequency where a transition occurs between low- and high-frequencies. The shape of this transition is affected by the value of the “form factor” M′.

Model 2, is a refined model, also derived by Lafarge [9, 18] in clarification of previous work by Pride et al. [46]. It essentially gives an improved definition of the shape of the transition frequency region, by adding the information provided by the thermal static tortuosity parameter ${\alpha }_0^{\mathrm{\prime }}$, which determines an additional form factor,

$$m^{\mathrm{\prime }}=\frac{M^{\mathrm{\prime }}}{4\left({\alpha }_0^{\mathrm{\prime }}-1\right)}.$$(117)

The refined model replaces the square root in (114) by the expression:

$$1-m^{\mathrm{\prime }}+m^{\mathrm{\prime }}{\left(1+\frac{M^{\mathrm{\prime }}}{2{m^{\mathrm{\prime }}}^2}\frac{-\mathrm{i\omega }}{{\nu }^{\mathrm{\prime }}}\frac{k_0^{\mathrm{\prime }}}{\varphi }\right)}^{1/2}.$$(118)

In doing so it is capable to retrieve the exact two first terms in (93). For an understanding of these models we refer to [4, 18, 28, 42], and more recently [9] (Appendix).

As we can see, Model 2 performs better than Model 1, however, small discrepancies subsist. The discrepancies are magnified in the Argand Plot representation, which particularly emphasizes the transition region. A direct look at the quantity of interest, the effective dynamic apparent modulus as a function of frequency, is given in Figures 4 and 5 where we plot the real and imaginary parts of the normalized bulk modulus versus reduced frequency $\mathrm{\Omega }=\omega {{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2/8{\nu }^{\mathrm{\prime }}$. (This choice of the reduced frequency is because of the next illustration, see below). On these figures the number of walkers is N = 10⁷ 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^{\mathrm{\prime }}$ and ${\alpha }_0^{\mathrm{\prime }}$ is that the Brownian motion method is validated. When the number of random walkers is not too high, it is possible to have very fast but not very precise estimates. The results presented above are prepared for the cylindrical circular pores, which is a very simple shape. One could check similarly that for other or more complex shapes for which we have the analytical solution, such as 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\sqrt{(1-\varphi )/\pi }$ 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₁, R₂, …, 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₂, …, r_n, …] were found to never go outside the 21 L × 21 L region, with typical drift, before reaching the pore walls, observed to be on the order of L. Very rarely, a drift as large as 8 L was observed in the random geometry. Note that the drift values give an idea of the representative cell dimensions of the material’s microstructure. Note also that, by eliminating the overlap between fibers, we automatically ensure that the porosity ϕ and the thermal characteristic length Λ′ in both geometries, random and regular, are the same (for instance Λ′ is here automatically Rϕ/(1 − ϕ)). Thus, using the reduced frequency, $\mathrm{\Omega }=\omega {{\mathrm{\Lambda }}^{\mathrm{\prime }}}^2/8{\nu }^{\mathrm{\prime }}$, we know that we should get the exact same frequency dependence for both geometries in the HF asymptotic limit. All calculations are done for the same 51 reduced frequencies as before.

It will be shown in this illustration that the irregular position of the fibers in real glass wool has a direct influence on the thermal relaxation of the effective bulk modulus, i.e. the way of moving from a relaxed isothermal state to a frozen adiabatic state. Note that due to the large density of glass compared to air, a porosity as high as ϕ = 0.97 is still a safe value to ensure that the fibers essentially remain at ambient temperature, justifying the assumed boundary conditions (28). The real part of the dynamic bulk modulus for the regular and “random” settings at porosity ϕ = 0.97 is plotted in Figure 7 versus the reduced frequency. The corresponding imaginary part is plotted in Figure 8. The predictions of Model 1 and Model 2 are also indicated, using as input parameters for $k_0^{\mathrm{\prime }}$ and ${\alpha }_0^{\mathrm{\prime }}$ the values inferred from the random walk algorithm, that were calculated beforehand with the same number of trials, N = 10⁷, using (94) and (95).

We note that we have at our disposal prior determinations of the parameters $k_0^{\mathrm{\prime }}$ and α₀, 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⁷ random walkers and by systematically taking 150 steps for one construction (rare events were observed where the previous value of 100 steps was insufficient to reach the boundary), our values are, $k_0^{\mathrm{\prime }}/\varphi =0.17141L^2$ and ${\alpha }_0^{\mathrm{\prime }}=1.0660$. Those of Cortis were, $k_0^{\mathrm{\prime }}/\varphi =0.17144L^2$ and ${\alpha }_0^{\mathrm{\prime }}=1.0654$. The differences concern digits that are not significant within the given FreeFM finite elements calculations.

In the random case the value of the static thermal permability is found to be approximately two times (1.87) that of the regular arrangement (namely 0.3197 L² versus 0.1714 L²). Thus the transition frequency is shifted to lower frequencies. The transition is steeper in the regular case, corresponding to a geometry where microscopic scales are less distributed than in the random case. This effect is clearly better described by an enhanced model where a second form factor (117) allows accounting for the broadening of the transition. Indeed, it can be noted that the value of the thermal tortuosity ${\alpha }_0^{\mathrm{\prime }}$ which is a measure of disorder (we have the known formula, ${\alpha }_0^{\mathrm{\prime }}=\langle {\tau }_0^2{\rangle }_p/\langle {\tau }_0{\rangle }_p^2,$ see [18] or Appendix in [9], where τ₀ 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.1671 versus 0.5205). This is why model 1, which assumes m′ = 1, is significantly less precise in the random case. Model 2 is also less precise, but to a much lesser extent. The introduction of the low-frequency information ${\alpha }_0^{\mathrm{\prime }}$ leads to an improvement of the description of the relaxation transition. Clearly, the proposed random walk method provides an efficient way to numerically calculate quasi-exactly the dynamic bulk modulus in a given geometry. Model 2 was found to be very close to the exact result.

7 Conclusion

A very simple simulation technique utilizing the known properties of Brownian motion with radioactive decay has been proposed to estimate, within the framework of Biot theory, the dynamic bulk modulus of a gas filling a porous medium. This technique has been validated on the simplest case of a cylindrical circular tube, and shown to provide a convenient way to calculate thermal relaxation in arbitrary geometries. The case of regular and random parallel fiber arrangement is considered because it corresponds to two extreme cases of glass wool microgeometry. The results of the quasi exact Brownian motion determination of the relaxation pattern of the function K_f(ω) highlight a simple relaxation function model ([9], Appendix) in terms of the following three pore space geometric parameters: static thermal permeability $k_0^{\mathrm{\prime }}$ divided by porosity ϕ, thermal tortuosity ${\alpha }_0^{\mathrm{\prime }}$, and thermal characteristic length Λ′. They show how the fiber arrangement is reflected on the shape of the relaxation transition.

The proposed simulation technique can provide useful clues on questions relative to the notion of macroscopic averaging. It is easy to implement and can be directly applied to all problems related to diffusion-controlled reactions between sinks, where the diffusing species are added in volume by using a spatially uniform, arbitrarily time-variable, 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.

Conflict of interest

The authors declare that they have no conflicts of interest in relation to this article.

Acknowledgments

The authors would like to thank Professors Evgeny A. Pchelintsev and Claude Depollier for sending them a copy of the article by Pontrjagin et al. [17] on which the method in this paper is based. N. Nemati acknowledges financial support from a French government Grant managed by ANR within the framework of the National Program Investments for the Future ANR-11-LABX-0022-01.

Data availability statement

The research data associated with this article are included within the article.

Appendix A

Additional comments on incompressibility, pressure, mean pressure, and averages

As we are concerned here only with rigid-framed materials and motion in the air, we introduce a macroscopic averaging operation 〈 〉(x) in the air, as a total–volume coarse graining average, $\langle ·\rangle \left(\mathit{x}\right)=\frac{1}{V}{\int }_{V_p}\mathrm{}·\mathrm{d}{V}_p$, given by the local air–phase integral of the considered quantities, divided by the corresponding total (air–phase + solid–phase) volume, where V_p(x) is the pore-space volume (the fluid volume) and V the total volume, in the representative volume (averaging sphere) of central position x. An averaging with normalization taken with respect to air-volume V_p instead of total-volume V, say $\langle ·{\rangle }_p\left(\mathit{x}\right)=\frac{1}{V_p\left(\mathit{x}\right)}{\int }_{V_p}\mathrm{}·\mathrm{d}{V}_p$, where index p over 〈 〉 is for pore-space, may be preferable when willing to have averaged values representative of the values in the fluid phase only. (Note that by abuse of language we use the same notation for the fluid domain V_p(x) and its corresponding volume; See also the two paragraphs before equation (1) for more details, and our two different notations of microscopic (r) and macroscopic (x) positions: The fluid volume element dV_p is d³r whereas x is reserved for the central position of an averaging sphere).

As we assume that Biot theory is a relevant description of the actual physics, the fluid and solid motions are, as insisted in the text, long-wavelength quasi-incompressible motions at the pore scale (i.e. in an averaging sphere). In particular it means here that, the microscopic excess-pressure acoustic field in the saturating fluid (the ambient air, here), p, is very nearly time-variable pore-scale constant. No important difference is, quantitatively, to be made between the microscopic, p, and spatial averaged macroscopic, P, excess-pressure. The quantitative double-scale asymptotic homogenization-theory justification for this [4, 48–50] is pertinent but not fully exact as explained in Sections 2 and 3. The latter sections provide qualitative justification of the incompressibility and approximate correspondence between p and P (not ∇p and ∇P!), as mere consequences of the absence of very different pore sizes (which excludes local resonance and pressure diffusion process [51]); see also the discussion around equations (7.160)–(7.161) in [9] (Appendix). The very small variations p − P that occur genuinely at the pore scale are removed in the macroscopic average P. Not explicitly mentioned in literature, within Biot’s theory the latter can be synonymously defined as: (i) the local direct coarse-grained average, P = 〈p〉_p (here local refers to making the average around a chosen central position x and is not to be confused with the notion of local and non-local theory), or, (ii) the local indirect coarse-grained average given by, 〈pv_i〉 = P〈v_i〉, where v_i is the ith component of the microscopic velocity in the fluid, i.e. this time, the Heaviside-Poynting or Umov averaged pressure (see Refs. [8, 9]) that can also be termed as a lumped pressure (see [52]). In this indirect definition of P the average is made indifferently with total-volume normalization, 〈 〉, or with air-phase normalization, 〈 〉_p (both give same result for P, as we have the general identity, 〈·〉(x) = ϕ(x)〈·〉_p(x), with ϕ(x) = V_p(x)/V the porosity). Here and in the text, for shortness, we sometimes omit to mention “excess” or “over” before “pressure” for p or P, simply called pressures. We assert that as long as Biot theory remains relevant, i.e. the fluid/solid incompressibility applies, both definitions of macroscopic pressure P(x), direct and indirect, are completely equivalent and indistinguishable. We note that, in what has been said above we have interpreted our averages 〈 〉(x) or 〈 〉_p(x) as volume averages performed at some central position x in a given sample (the 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 〈 〉(x) or 〈 〉_p(x) are more simply viewed as the expectation value of the considered quantity at the position x, in an ensemble of realizations of the macroscopic medium; the first is when the quantities are extended to zero in the solid and all realizations are taken into consideration, obtaining a total volume normalization, the second is when we discard realizations with x lying in the solid, obtaining a 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 〈 〉, 〈 〉_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: −〈pv_i〉 = H_ij〈v_j〉, that will provide a relevant general definition of a macroscopic symmetric stress field H_ij in the fluid, in the case of 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.

Appendix B

Additional comments on the constraints of Biot’s theory

Loosing sight of the semi-empirical origin of Biot theory it is sometimes written (e.g. [48]) that this theory is rigorously derived from the microscopic equations by the method of two-scale asymptotic homogenization developed and studied by Sanchez-Palencia [53], Keller [54], Bensoussan et al. [55], and others. However, this view is a misinterpretation of the character of this homogenization method, which automatically introduces the assumption of fluid/solid incompressibility as a result of its own procedure that hypothesizes the existence of one characteristic pore size, in microgeometry. The two-scale homogenization method does not derive Biot theory from microscopic equations only, but from the microgeometric constraints on which it is based.

Biot theory, as shown by the previously mentioned absence of spatial dispersion, and the related crucial role played by the incompressibility hypothesis, is not general in some potentially important aspects, such as the definitions with which it can be expressed (like the systematic assimilation of macroscopic values with direct volume averages, see comments in Appendix A), and the quantities involved. It constitutes a special (though widely applicable) simplified limit (non-local in time but local in space) of a much more general theory (applicable to much wider classes of microgeometries) non-local in both time and space. This more general theory will possibly lead, in other simplified limits, to more general local (local in space) descriptions. The fully general non-local theory or its more general local simplifications would not satisfy the characteristic assumption of fluid/solid incompressibility²⁴ and would lead to much richer behavior than conventional Biot theory allows.

The characteristic reduction in the possibility of temporal dispersion in Biot model is seen in the fact highlighted by Johnson [28] (Appendix A) and Lafarge [18] and [4] (Appendix C), that the singularities of the response functions α(ω), k(ω), α′(ω), k′(ω), in the complex ω-plane, lie only on the imaginary half-axis, as opposed to the full half-plane allowed by the general causality principle; this, also, expresses in the monotonic behaviors mentioned in Section 3. Using a divergence-free representation of all fluid displacements combined with a uniform representation of all solid displacements, at the pore scale, Biot’s model cannot account for spatial dispersion, at all. The latter remains entirely outside its framework. With spatial dispersion present as soon as significant divergent motions appear at pore scale [9], the usual simple basis for the definition of local response factors (that is detailed in Sect. 3), simply vanishes. In these circumstances the operators and field quantities must be thoroughly generalized, possibly enhancing even their tensor rank (see Appendix A and ^{Footnotes 1} and ¹⁶). In a forthcoming work (the first mentioned in ^{Footnote 1}), we will illustrate how an imbalanced divergent fluid motion at the small scale can generate a much more general temporal dispersion in connection with underlying presence of spatial dispersion.

In general it is the relatively good satisfaction of the incompressibility of solid and fluid motion at long wavelengths in a variety of microgeometries that explains the large scope and experimental success of Biot theory. These simplifications are expected to hold true in the real propagation problem and to justify Biot theory, not only when there is a well-defined pore size, but more generally as long as the complications in the geometries do not lead to locally resonant structures [8, 9] or other effects related to fluid flow divergence, such as pressure diffusion [51]. In all these cases, both the long-wavelength spatial dispersion (related to the non-uniform or divergent parts of the solid or the fluid motions on the pore scale) and the mentioned additional temporal dispersion (also triggered by these motions) are practically absent. This is an expression of the fact that, compared to the nearly uniform (solid) or non-divergent (fluid) motions on the pore scale in Biot theory, the actual non-uniform/divergent distributed motions of the solid/fluid on the pore scale usually have a very weak amplitude. Their complete absence in Biot theory is usually not a problem, but a reasonable simplification.

Therefore, instead of assuming a very well-defined pore size, as is the case with two-scale asymptotic homogenization, it would be more appropriate to characterize the physical position of Biot theory by recognizing that it corresponds to situations in which one can neglect the long-wavelength spatial dispersion phenomena and a corresponding part of the temporal dispersion, which allows considerable simplifications. We note that the fact that the asymptotic two-scale homogenization method automatically preserves, to a first approximation, the incompressibility of the solid and fluid motion at the pore scale is an indication that this method is not suitable to describe the spatial dispersion and part of the temporal dispersion. For this reason, the higher order terms in this method do not accurately describe the onset of the spatial dispersion corrections that occur at decreasing wavelengths, regardless of the microgeometries, simple or complex, and the part of the temporal dispersion effects associated with this onset. This indicates the limited physical consistency of the expansion principle used in this method. Without the unnecessary detour via this method, as made in [4], we can base our description of the function K_f(ω, x) directly on the incompressibility of fluid and solid motion, which we introduce here as an explicit hypothesis/simplification justified by the type of geometry considered.

References

---

All Figures