## Abstract

The presence of a particulate phase defines the effective response of a suspension to changes of the average heat flux. Using the random-point approximation we show that within the first order in the concentration, one needs to solve the problem for the temperature field created by a single inclusion in a matrix subject to an unsteady temperature gradient at infinity. We solve this problem by means of Laplace transform and use the solution as the first-order kernel in the functional expansion. From this kernel, we find the statistical average for the heat flux, which turns out to be a memory integral of the spatially averaged time-dependent temperature gradient. Thus, we discover that the constructive relationship between the average flux and averaged temperature gradient is not local in time, but rather involves a convolution integral that represents the memory due to the heterogeneity of the system. This is a novel result, which *inter alia* gives a rigorous justification to the usage of generalizations of the heat conduction law involving fractional time derivatives. The decay of the kernel is very close to *t*^{−1/2} for dimensionless times lesser than one and abruptly changes to *t*^{−3/2} for larger times.

## 1. Introduction

Estimating the effective transport coefficients of heterogeneous media is important for technology, and has attracted vigorous investigation (see the detailed accounts in Beran (1968) and Miton (2002), and especially the compilation of review articles in Markov & Preziosi (2000)). The most typical examples of such media are the particulate materials (suspensions) in which the second (particulate) phase comprises spherical particles (the *filler*) that are randomly dispersed throughout the continuous phase (the *matrix*). The different transport problems that can be considered for a suspension are the effective electric or heat conductivity, effective viscosity and effective elasticity. The effective thermal conductivity is currently attracting a considerable attention in connection with the nanofluids and their use as coolants (Vadasz *et al.* 2005; Vadasz 2006). Respectively, the mathematically analogous problem for the effective diffusivity of a particulate medium is of importance for modelling the suspensions of swollen microgel particles (Routh & Zimmerman 2003).

Different methodologies and techniques have been developed for evaluating the effective heat transport properties of composites. They range from purely empirical models to very detailed description of the field around individual particles and consecutive average. The averaging over volume and the notion of effective modulus was first introduced by Maxwell (1873: §314), who compared the potential created by *n* spheres, each of radius *a*, to the potential of a sphere that encompasses the swarm of small spheres and has an equivalent electric conductivity. Maxwell obtained the contribution to the effective conductivity of first order with respect to the volume fraction of the particulate phase. The same idea was applied by Einstein (1906) for computing the first order in volume fraction contribution to the effective viscosity of a suspension. For the elastic moduli of a suspension, the same approach was corroborated in Walpole (1972). Jeffrey (1973) argued that the method proposed by Maxwell can give correctly only the first order in the volume fraction and went on to discuss the statistical properties of the centres of spheres. He extended the arguments from Batchelor (1972) and Batchelor & Green (1972) and reached the conclusion that for the second-order approximation, the solution around two spheres is needed. The rigorous inclusion of the probability distributions of the centres of spheres for the case of effective heat conduction was introduced in Keller & Feuillebois (1993). A comprehensive review of the works on viscosity of suspensions can be found in Herczyński & Pieńkowska (1980). The effective constitutive properties of heterogeneous elastic materials with various shapes of the inclusions are exhaustively treated in Ostoja-Starzewski (1998).

Because of the complexity of the mathematical problems involved, different approximations are still being sought, e.g. limits on the moduli, etc. (e.g. Lipton & Vernescu 1996; Markov 1999; Miton 2002 for the current stage of the research). Usually, the temperature field is considered to be stationary with constant mean gradient, and the heat conduction problem is solved approximately. The obtained flux is averaged and as a result the overall relationship between heat flux and temperature gradient is established.

However, the problem for the average response of the material to a time varying spatially constant gradient does not seem to have attracted the same attention as the stationary case. Even a brief glance at the unsteady problem reveals that it is qualitatively different from the steady one, in the sense that the presence of an interface between the phases provides resistance to the propagation of heat and the average of the flux may not be immediately proportional to the current value of the average temperature gradient. Clearly, the described resistance will introduce some kind of a memory in the constitutive relation. The memory effect is well documented starting with the groundbreaking work of Basset (1888) on the resistance on a sphere in a viscous liquid. The propagation of heat through a dilute suspension of small particles is investigated in Coimbra & Rangel (2000) using the Duhamel’s integral, which approach appears to be more straightforward than the Laplace transform method for the case when there is no average gradient of the temperature field. The appearance of memory in heterogeneous materials has also been argued on the basis of the homogenization method in Tartar (1989). It seems important to investigate the memory effect in case of particulate materials subject to temperature field with constant gradient, which is the objective of the present paper.

A rigorous approach to constructing successive approximations to the random fields that describe a mechanical property of a composite (such as viscosity, elasticity, heat diffusivity, etc.) has been proposed in Christov (1981). The approach is based on an idea by Wiener in the late 1940s (see Wiener 1958) to expand the unknown random fields in functional series (see Volterra 1982 for the definition) in which the kernels are non-random functions and only the basis function contains the information about the randomness of the field. The original Volterra–Wiener Expansion (VWE) employs the Gaussian white noise as the basis function. The VWE can be made orthogonal by using the Hermite polynomials of the white-noise function, and the series is alternatively called ‘Wiener–Hermite Expansion’ (WHE). The WHE was applied to some turbulence models (e.g. Siegel *et al.* 1963) but faced formidable difficulties for nonlinear problems and required frequent renormalization (Meecham 1974). This shortcoming was explained in Christov (1981) by making the observation that the normally distributed Gaussian white noise has zero third moments, while most physical systems are far from normality, exhibiting significant third cumulants. Thus, the general theorem from Cameron & Martin (1947) that any random function can be expanded into WHE is of not much help, because WHE is practically divergent when the kernels are identified from nonlinear hierarchical systems.

The most natural conjecture for a particulate material is that the centres of the inclusions are represented by a system of random points, while the effect of each inclusion is defined by the mechanical properties and the geometry of the inclusion, i.e. the local effect is deterministic. The superposition of non-random structures that appear in random positions of space or moments of time is called random point function/process (RPF; see Snyder 1975). This led one of the present authors to suggest that VWE with a RPF basis opens the way for rigorous construction of consecutive approximations to the problem of mechanical behaviour of composites. Naturally, the simplest RPF function is the Poisson process (Snyder 1975), and the proposal forwarded in Christov (1981) was to use a Poisson basis. Actually, the first paper in which the Poisson process was used as the basis function of the functional expansion was Ogura (1972), where it was demonstrated that the Charlier polynomials render the functional series orthogonal. However, no applications of the Poisson–Wiener expansion were attempted in Ogura (1972), and the exclusive role of these series for nonlinear problems was argued in Christov (1981).

It turned out that the plain Poisson basis function was a very good guess for the basis in modelling the chaotic behaviour of nonlinear dynamical systems, e.g. turbulence. It goes beyond the scope of the present paper to give an overview of the research in that direction. We refer the reader to Christov & Nartov (1992) for a comprehensive exposition or to the English language reviews (Christov 1990; Nartov & Christov 1995). The application to suspensions/emulsion of the random-point approximation (RPA) was sketched first in Christov (1983*b*, 1984*b*) (see also Christov (1983*a*) for treating the interaction between the phases in suspensions/emulsions).

The Poisson RPF is not very well suited for composite materials, because if the random centres are assumed to be completely independent statistically, then the inclusions can consequently overlap. However, the RPA is so well tailored to this kind of problem that the (anti)correlation of the random points needed to prevent the overlapping only contributes to the second order with respect to the volume fraction. This was the reason why in Christov & Markov (1985*b*) the RPFs of perfect disorder were introduced, for which the multi-point probability distributions are defined to be zero if the centres of the respective inclusions are closer than the sum of their radii. If the inclusions do not overlap, then the multi-point probability density is equal to unity, which corresponds to a Poisson distribution. The effective heat conduction of a monodisperse suspension was treated in Christov & Markov (1985*b*) and shown that the linear contribution of the filler to the effective modulus is equal exactly to Maxwell’s formula for the effective electric conductivity. Later on, the elastic problem was solved in Christov & Markov (1985*c*). More insight on the mechanisms that define the two-point probability distribution for spherical suspensions can be found in Torquato (1986). The role of the two-point correlation of arbitrary-shaped inclusions on the bounds for the effective coefficient is elucidated in Markov (1999).

After the RPA was generalized for the case of marked RPF in Christov (1985), the most general case of a polydisperse suspension of perfect disorder type became amenable to the VWE method (see Christov & Markov 1985*a*). It has been rigorously established in the cited works that the knowledge of the first kernel allows one to obtain the respective effective modulus with a first order of approximation in volume fraction *c* of the filler. Respectively, it has been proved that only if one has the two-sphere solution, then one can rigorously obtain the second-order approximation in the volume fraction, as the two-sphere solution gives precisely the second-order kernel in the formal VWE.

While, the one-sphere solution was known for many years, the two-sphere solution could not be obtained analytically. The numerical technique from Christov (1984*a*) has been recently implemented in Chowdhury & Christov (2009*a*), and the actual temperature distribution around two spheres has been found for a constant gradient at infinity. The contribution due to the ‘binary interaction’ (two-sphere solution) is of orders of magnitude smaller than the contribution of a single sphere. The case of very closely situated spheres has been interrogated in Chowdhury & Christov (2009*b*) and the two-sphere contribution has been shown to be less than 20 per cent of the one-sphere one. It is reasonable to expect that the relative importance of the second kernel will be similar in the case when the temperature gradient is function of time. This motivates the goal of the present work: to compute the contribution of the first kernel (one-sphere solution) for unsteady temperature gradient and to apply the RPA to finding the constitutive relation between the average heat flux and temperature gradient.

## 2. Medium with random discontinuous transport coefficient

We consider the two-phase medium that consists of an array of perfectly disordered spheres of random radii with thermal conductivity ϰ_{f} distributed in an unbounded matrix of conductivity ϰ_{m}. Thus, the suspension is a medium with a discontinuous random coefficient of thermal conductivity, which can be expressed in the following way:
2.1
where *h* is the shape function of a single spherical inclusion, and is the four-dimensional space which is the cartesian product of the three dimensional configurational (geometric) space and the set of the possible values of the random mark (so-called ‘mark space’). Without loosing the generality, we can assume that because if the set is any finite interval in , then we can augment the probability density with zero outside of the finite interval. The randomness shows up through the function
which is a random superposition of Dirac delta functions. Here *ξ*_{α} is the random position at which appears a sphere of radius *a*_{α}. For the case of spherical suspensions, the random mark has a clear physical meaning of the radius of the sphere that happens to be centred at the random position *ξ*_{α}. Function *ω*(** x**,

*a*) is called generalized random density function. Let

*γ*denote the number density of the system of random points, and

*P*(

*a*) is the probability density of the mark. When

*P*(

*a*)=

*δ*(

*a*−

*a*

_{0}), one encounters the case of monodisperse suspension containing spheres of the same radius

*a*

_{0}.

To forbid the spheres from intersecting, the multipoint probability distributions for the system of random points were stipulated in Christov & Markov (1985*a*) to obey the law:
2.2
The above expression for *Q* ensures that when the centres of two spheres are farther than the sum of their radii, they are statistically independent, and the basis function behaves as a Poissonian random function (as sometimes called the ‘perfectly random function’). In Christov & Markov (1985*a*), random functions of the type considered here were called ‘random functions of perfect disorder type’. Note that *Q*_{ik} is symmetric with respect to interchanging the arguments within the first and second pairs.

We observe here that the average volume of a single inclusion is the following mathematical expectation:
and then the concentration of the particulate phase is by definition *c*=*γ*〈*V* 〉. One should not confuse the notation for the average volume 〈*V* 〉 with the notation *V*_{a}=4*πa*^{3}/3, which stands for the volume of geometric domain that is encompassed by a single spherical inclusion of radius *a* centred at the origin.

In this paper, we will consider the equation of the heat diffusion
2.3
where the coefficient of thermal diffusivity ϰ(** x**) is a random function, given in equation (2.1). We will seek a solution of the problem in functional series.

## 3. Volterra–Wiener functional expansions

For the continuum considered here, the randomness of the solution is a result of the randomness of the coefficient. The natural approach to this kind of problem with known source of randomness is to consider a functional series with basis function that is related to the random coefficient (the random density function that defines the coefficient). This means that the temperature field can be considered as a Volterra–Wiener functional series in the products of the basis function with non-random kernels *T*_{n}. The direct application involves somewhat unwieldy expressions for the scalar products of the different products of function *ω*. To streamline the derivations, polynomials of *ω* can be used, for which the different terms are orthogonal to each other.

As already mentioned, to orthogonalize the functionals, one needs to use multivariate stochastic polynomials of the random density function *ω*. In Christov (1985), the Charlier polynomials were generalized for the case of the random density functions of perfect-disorder type. The basis functions are as follows (the first three are presented here)
3.1a
3.1b
We mention here that from equation (3.1a) we find . Then upon introducing the above expression for *ω* and for the concentration in equation (2.1), we find an alternative expression for the discontinuous transport coefficient, namely
3.2
where 〈ϰ〉=ϰ_{m}+*c*[[ϰ]] denotes the average conductivity of the medium. Note that this average is not equal to the effective heat conductivity coefficient. Respectively, the Volterra–Wiener series for the temperature field adopt the form
3.3
where *T*_{0}(** x**)≡〈

*T*(

**,**

*x**t*)〉 is the ensemble-averaged temperature field. We note here that for finding the effective coefficient, it is enough to consider the case when ∇

*T*

_{0}(

**)=〈∇**

*x**T*(

**)〉=**

*x***(**

*G**t*), where

**(**

*G**t*) is some spatially constant vector gradient at infinity. It is supposed to be a ‘slow’ function of time in the sense that it changes slowly compared to the ‘diffusion time’

*a*

^{2}/ϰ

_{f}. The meaning of this requirement is paramount to any treatment of the effective properties of a heterogeneous material: the overall gradients and fluxes should not change rapidly in time (or abruptly in space) around a single inclusion.

Let us note that it was shown in Christov & Markov (1985*a*,*b*) that the functional series is *virial* in the sense that the integral containing the *n*th kernel contributes to the averaged characteristics a quantity of order *γ*^{n}, where *γ* is the above-mentioned number density. Actually, after acknowledging the connection of *γ* to *c*, we can reformulate the property of ‘viriality’ and to claim that each term contributes a quantity of order of *c*^{n}. Since *c*<1, the virial property of the functional series ensures their convergence. Despite the fact that for very concentrated suspensions the convergence may be slow, the viriality is an indispensable property, because it provides for a rigorous perturbation methodology.

Combining equations (3.2) and (3.3), we find the expression for the heat flux as a random function within the adopted order of approximation *O*(*γ*^{3}):
3.4
Here the notation ∇_{x} stands for the gradient of *T*_{2} with respect to ** x**, which gradient is obtained from the gradients with respect to its two spatial vector arguments.

In order to obtain an expression for the averaged heat flux, we need the expressions for the average values of the products from the above cited works. For the first few we have within the adopted order *O*(*γ*^{3}) the following:
3.5a
3.5b
3.5c
3.5d
3.5e

Equation (3.5e) is obtained from equation (3.5d) after acknowledging the connection between *C*^{(2)} and *C*^{(1)} from equation (3.1b), orthogonality of *C*^{(2)} and *C*^{(1)}, as well as the fact that *C*^{(2)} is a centred random function.

## 4. Governing equation for kernels—asymptotic uncoupling

The basic problem when employing functional expansions is the identification of the respective kernels. We follow Christov & Markov (1985*a*), which is a less accessible publication. We mention that the gist of the method can be found also in Christov & Markov (1985*b*), but for the limiting case of monodisperse suspensions. We assume that the average gradient does not depend on the spatial coordinate, namely ∇*T*_{0}(** x**)=

**(**

*G**t*). Taking the average of equation (2.3) with equation (3.4) acknowledged, gives At this junction, in order to simplify the treatment of the double integral, we define an auxiliary function as follows (see Chowdhury & Christov in press): 4.1 where the subscript

*P*refers to the fact that the average is taken with respect to the probability distribution

*P*(

*b*). Note that is introduced mostly for simplifying the notation in the general formulas. The integration over

*b*is not an obvious procedure, and very often it is better deferred to the stage where the integration with respect to

**has already been performed (see Chowdhury & Christov in press). The meaning of function**

*y**S*is that when used as a kernel in an integral, it parametrizes the effect of the spheres being restricted from overlapping. In what follows we call

*S*the ‘interaction function’. As it will become evident later, the statistical average of the interaction function plays an important role. For this reason we introduce the special notation .

In addition we introduce the ‘volume average’ of the gradient of *T*_{1}
4.2
over the volume of a single sphere. In some cases one may find it helpful to use the divergence theorem and to render the above integral to a surface integral, but this goes beyond the scope of the present paper. In a similar fashion we define
4.3
where the meaning is that the average is taken over all different volumes of sphere whose radii are distributed with probability *P*. Then the first equation for the kernel can be rewritten as (note that ∇*T*_{0}(** x**)≡

**(**

*G**t*),

*c*=

*γ*〈

*V*〉): 4.4a

Following the procedure outlined in Christov (1981) and developed further in Christov & Markov (1985*b*), Christov (1985) and Christov & Markov (1985*a*), we multiply equations (2.3) and (3.4) by the first-order Charlier polynomial , and obtain within the order of approximation *O*(*γ*^{3}) the following equation for the non-random kernel *T*_{1}(** x**;

*a*): Before proceeding further, we need to manipulate the integral on the fourth line of the above equation taking into account the property of function

*Q*(

*ξ*_{1},

*ξ*_{1};

*a*

_{1},

*a*

_{1}), namely where the localization properties of function

*h*are taken into account. Now, we are equipped to derive the most convenient form of the equation for the kernel

*T*

_{1}: 4.4b

Following the general procedure, we multiply equations (2.3) and (3.4) by the second-order Charlier polynomial to obtain 4.4c

As one should have expected, the method of functional expansion—as all the other methods—leads to an infinite hierarchy of equations. What makes the difference here, however, is that the expansion under consideration is virial and the higher-order functionals contribute to the averaged characteristics quantity proportional to the respective power of *γ* (or which is the same – to *c*). Respectively the *n*th order kernel appears in the hierarchy always multiplied by *γ*^{n} (or *c*^{n}). As far as *c*<1 one can effectively decouple the hierarchy by expanding the kernels into power series with respect to the concentration *c*.
In order to display the gist of the method it is enough to consider the first two kernels. So we discard terms proportionals to *γ*^{n},*n*≥3, and by combining the terms of the like powers we arrive to the following systems of equations (with 〈ϰ〉=ϰ_{m}+*c*[[ϰ]] acknowledged):
4.5a
4.5b
4.5c
The equations contain discontinuous coefficients and their solutions have derivatives that are generalized functions. In what follows we concern ourselves only with the first-order solution. Equation (4.5a) for the kernel *T*_{10} is nothing else but the equation for the heat propagation.

## 5. First-order time-dependent problem

Equation (2.3) contains a discontinuous coefficient, and as such can be considered as two equations for the functions and for the distribution outside and inside of the sphere, respectively: 5.1a 5.1b

The problem that is relevant to the theory of effective heat conduction of suspensions requires that the temperature field at infinity is a linear function of the spatial coordinates (constant temperature gradient). However, in the present work, the gradient at infinity can be a function of time, i.e. the temperature is proportional to [** G**(

*t*)+

*G*_{0}]⋅

**for . We choose to separate the ‘pure’ unsteady part of the temperature gradient from the steady part**

*x*

*G*_{0}. Then we can consider function

*T*

_{10}(

**,**

*x**t*) from equation (2.3) as a perturbation, i.e. it decays at infinity 5.2

The fact that the sought function is considered as a perturbation, requires that in the boundary conditions we ensure that both the temperature and the total flux (perturbation added to the constant gradient) are continuous across the sphere’s boundary, |** x**|=

*a*: 5.3

Thus, equations (5.1) and the conditions—equations (5.2), (5.3)—form the boundary value problem under investigation. Before proceeding further, we mention that the stationary problem (when ** G**(

*t*)=0) has the following well-known solution 5.4 Then, we can introduce a new function

*u*(

**,**

*x**t*)=

*T*

_{10}(

**,**

*x**t*)−

*u*

_{0}(

**), for which we have the same equations (5.1), but with the following jump conditions: 5.5 We assume here that a constant gradient has acted long enough before the initial moment**

*x**t*=

*t*

_{0}to have been able to create the distribution presented in equation (5.4). The first moment of time for which the temperature gradient will be allowed to change is

*t*=

*t*

_{0}. Thus,

*u*(

**,**

*x**t*

_{0})=0 (and

**(**

*G**t*

_{0})=0) is the initial condition for the new function. Without loosing the generality, we can set

*t*

_{0}=0.

A convenient method for solving equation (2.3) is to take the Laplace transform with respect to *t*. We denote the Laplace transform of function *u*(** x**,

*t*) and its derivative by which allows us to rewrite equation (2.3) as follows 5.6 Applying the Laplace transform to the gradient, we find

**=∫**

*G*^{∞}

_{0}

*e*

^{−st}

**(**

*G**t*)

*dt*. Thus, we arrive at a problem for the Laplace transform of the sought function, which is very similar to the steady one, with the difference that it contains also a Helmholtz term. To solve equation (5.6), we make the same stipulation as for the stationary case, namely

*U*(

**;**

*x**s*)=

*f*(

*r*;

*s*)

**⋅**

*G***, where**

*x**r*=|

**|. Then from the Laplace operator we find the following ODE for**

*x**f*: 5.7 where the prime stands for a differentiation of function

*f*with respect to the argument

*r*. Note that the dependence of

*f*on

*s*is as of on a parameter (no derivatives are considered). The standard way to solve the equation with radial symmetry equation (5.7) is to introduce a new function

*W*=

*fr*, which upon substitution in equation (5.7) renders the latter to an equation with constant coefficients, whose solutions are , where .

Now we return to function *f* and retain in the two regions (inside and outside of the sphere) only the fundamental solutions that do not diverge. The expression for *f*(*r*,*s*) gives the solution of the homogeneous equation inside and outside of the sphere:
5.8a
5.8b
where *A* and *B* are two unknown constants to be defined by the boundary conditions.

Now, by taking the Laplace transform of the boundary conditions equation (5.3) we find the conditions on *U* to be satisfied at |** x**|=

*a*, namely

After introducing the expressions for the general solution (5.8) we find:
5.9a
5.9b
taken at *r*=*a*. The first of these equations, equation (5.9a), is satisfied if
5.10
To satisfy the second equation (5.9b), we first show the following identities
5.11a
5.11b
Here *ϕ*(*r*) is an arbitrary function of the variable *r*=|** x**|, and we have used the fact that does not depend on the spatial variables. Now we can use the last result and to plug for

*ϕ*each of the function on the left- and right-hand sides of equation (5.9b) in order to get the second condition for

*A*and

*B*. Thus for

*r*=0 we have 5.12 where a term has been cancelled on both sides of the equation. Upon applying the operators on both sides of equation (5.12), we obtain the second condition for

*A*and

*B*. 5.13

Now by using the system of equations (5.10), (5.13), we can solve for *A* and *B*
5.14a
5.14b

Thus, the Laplace transform of the temperature field has been found.

## 6. The effective heat flux

Taking the average of Laplace transform of equation (3.4) we find 6.1

In order to find the effective heat flux, we need the integral over the inclusion from the gradient of the temperature, where *D*_{a} is the sphere of radius *a*. Respectively, is given by equation (5.8a) in which the equation (5.14a) is acknowledged. Now, by substituting or *ϕ*(*r*)=−*β*/*s**G*_{0}⋅** x**, we find with the help of equation (5.11b) that
6.2

We introduce spherical coordinates as usual: ** x**=

*r*sin

*θ*cos

*φ*

**+**

*i**r*sin

*θ*sin

*φ*

**+**

*j**r*cos

*θ*

**, where**

*k***,**

*i***,**

*j***are the unit vectors of a Cartesian coordinate system. Then where some trigonometric integrals have been acknowledged. Respectively, and By using all the results we can obtain the final expression for the integral over the sphere: 6.3 where 6.4 Thus, we have found the following constitutive relation for the overall heat conduction of the composite in terms of the Laplace transforms: 6.5 where ϰ**

*k*_{eff}is the Maxwell effective modulus for stationary heat conduction thorough a suspension (see Christov & Markov 1985

*b*). After applying the inverse Laplace transform, the last equality gives the effective constitutive law for the heat conduction of a suspension: 6.6 where

*K*(

*t*)=ℒ

^{−1}[

*H*(

*s*)] is the inverse Laplace transform of the function

*H*(

*s*). Here is the unsteady part of the flux, generated by the unsteady part of the temperature gradient. The law (6.6) naturally encompasses the Maxwell result for the steady heat conduction throughout a suspension, which is recovered by setting

**≡**

*G***0**. Naturally, the deviation from the Fourier law is proportional to the jump of the heat coefficient [[ϰ]].

There was a benefit in deriving the above expressions in dimensional form in order to keep in the front the physical meaning of the different terms, but it is more convenient to present the parametric study of the memory relation in dimensionless variables. We introduce dimensionless time and Laplace parameter using the so-called ‘diffusion time’ *τ*=*a*^{2}/ϰ_{f}. This means that one cannot consider a filler with zero conductivity ϰ_{f}=0, i.e. a suspension of ‘holes’ in a matrix of given heat conductivity ϰ_{m}. Such a case has to be considered separately and goes beyond the scope of the present paper. Thus, we introduce the dimensionless variables *t*′=*t*ϰ_{f}/*a*^{2}, *s*′=*sa*^{2}/ϰ_{f}, . Then
6.7
and the dimensionless memory kernel is defined as *K*′(*t*′)=ℒ^{−1}[*H*′(*s*′)].

Here and henceforth the primes will be omitted, and dimensionless variables will be presumed in all expressions.

## 7. Approximate expression for the memory kernel

We present *H*(*s*) in a form that is better suited for inverse Laplace transform
Denote . Then *H*(*s*) adopts the form

Now, by applying the inverse Laplace transform to *H*(*s*) and we find
7.1

For the in-depth validation of the analytical result shown in equation (7.1), we assess the error introduced when the series is truncated at certain finite *n*. First, we find a numerical solution of the inverse Laplace transform of function *H*(*s*). We tried a couple of different algorithms available in the literature. For our problem, the most reliable and dependable appears to be the algorithm (Valkó & Abate 2009); see, also the original works of Abate & Valkó (2004) and Valkó & Abate (2004). For a reasonable physical case with ϰ_{f}=2 m s^{−2} and ϰ_{m}=1 m s^{−2} when , we conducted a number of numerical experiments with the algorithm in order to find the optimal parameters. The results for the numerically inverted *H*(*s*) are presented by the dotted line in figure 1*a*. Then we computed *K*(*t*) from equation (7.1) with different number of terms. The results are presented by the lines of different colors (different lengths of the dashes) in figure 1. In figure 1*a*, the comparison is shown with the numerical solution for the range 0<*t*<1000. One can see that taking just seven terms in the series provides a high accuracy of the truncated series. This conclusion is further verified in figure 1*b*, where it is seen that the difference between an expansion with six and seven terms is very small up to *t*=30 000. This means that for all practical purposes, one can use the expression for the kernel *K*(*t*) of the memory integral as given by equation (7.1) with merely seven terms. Note that for a larger *δ*, one may need a larger number of terms.

## 8. Discussion

Equation (6.6) presents a very important result: the constitutive law for the heat conduction of a composite, whose both phases are governed by the Fourier law, is by necessity a non-Fourier material involving memory (‘retardation’) of the temperature gradient due to the convolution integral with kernel *K*(*t*) in equation (7.1). The conduction law for a composite of two Fourier materials is non-Fourier because of the enhancement/obstruction of the heat propagation due to the boundaries between the phases. We note that a general relation like the one in equation (6.6) can be rewritten in different forms, if one moves some terms to the left-hand side and then takes the inverse Laplace transform. A significant insight can be gained from such kind of reformulations, and they will be performed elsewhere. Here, we just focus on the most basic consequences of our results.

It is important to emphasize that in equation (7.1), apart from the memory integral, a term giving a linear dependence of the flux on the temperature gradient is also present. This means that the influence of the value of the gradient at the current moment of time is very strong. Our law can be incorporated into the general description from Gurtin & Pipkin (1968), provided a delta function is added to the kernel −*K*, which is another way to express the importance of the value of the gradient in the present moment. The most salient trait of the kernel *K*(*t*) is that it is singular at *t*=0 and it has the opposite sign of the linear term, acting to diminish the importance of the gradient at the current moment of time. Yet, the singularity being of type *t*^{−1/2} is not really strong, i.e. there is a significant contribution from the moments of time prior to the current moment. In this sense, the composite material has rather ‘long memory’.

The second important trait of the kernel *K*(*t*) is that around dimensionless time *t*=1, the decay of the kernel abruptly changes from *t*^{−1/2} to *t*^{−3/2}. This means that the predominant effect of the presence of a memory is felt for . For longer time the memory fades quicker. It is to be mentioned here that the memory effect in Coimbra & Rangel (2000) decays uniformly as *t*^{−1/2}. The difference stems from the fact that we are dealing here with a gradient field, which interacts with the uniform field and creates the faster decay of the kernel at infinity. Physically speaking, this is related to the mean size of the particles, more specifically to the ‘diffusion time’, *a*^{2}/ϰ_{m} needed for the heat changes to be felt considerable at a distance commensurate with the mean radius of the particles comprising the filler.

An interesting limiting case is encountered for *δ*≫1 when the medium contains spheres with infinite conductivity. In this case
8.1
where stands for the Riemann–Liouville fractional half integral (Podlubny 1998). Note that infers and 〈ϰ〉≃*c*[[ϰ]]. Then one can introduce a normalized flux, and to arrive at the following alternative expression for the constitutive law:
8.2
which resembles the Maxwell–Cattaneo (MC) law but now the relaxation is given by the half time derivative rather than the first derivative, and additional fractional retardation is present. The MC law is a phenomenological assumption that allows one to account for the finite speed of the propagation of heat disturbances (see the illuminating review Joseph & Preziosi 1989, 1990). The issues connected with the material invariance in the case of heat conduction in a moving continuum were discussed in Christov & Jordan (2005) and Christov (2009), where the role of the convective terms was identified. An important justification for the MC law has been recently provided in Ostoja-Starzewski (2009) from the point of view of molecular dynamics, albeit the effects in homogeneous fluids of relatively light molecules are expected to be quantitatively small. As a result of the present investigation, a heterogeneous mixture with particles of much larger size is expected to possess memory in the heat conduction even without investigating the molecular dynamics in detail. This is a very important result, because it provides a justification for using fractional hyperbolic heat conduction laws for suspensions of nanoparticles as discussed in Vadasz & Govender (2010).

## 9. Conclusion

In the present paper, the method of stochastic functional expansions with random-point basis function is applied to the case of the heat conduction of a particulate medium subjected to a time-dependent spatially constant mean temperature gradient. It is shown that within the first order of approximation with respect to the concentration, the equation for the kernel of the first-order functional integral is the equation of the disturbance introduced by a single spherical inclusion in a matrix subjected to a time varying spatially constant gradient at infinity. The resulting initial-boundary value problem is solved by means of Laplace transform. After ensemble averaging of this solution, the effective connection between the time varying gradient and the flux is established. It turns out that the effective law involves a memory integral (retardation) of the temperature gradient. A method based on infinite series expansion is used to invert the Laplace transform, and an approximate expression for the memory kernel is found.

For the interesting limiting case of filler material with infinite conductivity, the memory integral becomes merely the Riemann–Liouville half integral, and a simple inversion of the relationship gives a Maxwell–Cattaneo type law of heat conduction, in which the relaxation of the flux is given by the half fractional time derivative.

Thus, the key result of the present paper shows rigorously for dilute suspensions that the constitutive relation for the overall heat conduction of a suspension contains memory terms even if the two materials involved are obeying the Fourier law.

## Acknowledgements

The work of CIC was supported, in part, by an 2009 ASEF/ONR Summer Faculty Fellowship. CIC gratefully acknowledges the hospitality of Dr Stanley Chin-Bing and the members of his group at the U.S. Stennis Space Center, Mississippi.

## Footnotes

This paper is dedicated to the memory of Professor Konstantin Markov on the occasion of the 65th anniversary of his birth.

- Received March 4, 2010.
- Accepted April 16, 2010.

- © 2010 The Royal Society