Steady-state heat transfer in microcracked media

. Material behaviour is often affected by the heterogeneities existing at the microscopic level. Especially the presence of cracks, voids, etc collectively known as defects, can play a major role in their overall response. Homogenization can be used to study the influence of these heterogeneities and also to estimate the effective properties of a given material. Several research works have been dedicated to determining the elastic behaviour of microcracked media. Yet, thermal properties are not investigated as much. Moreover, the question of unilateral effect (opening/closing of cracks) still remains an important issue. So, this paper aims to provide the effective thermal conductivity of 2D microcracked media with arbitrarily orientated cracks, either open or closed. With the help of Eshelby-like approach, homogenization schemes (dilute and Mori-Tanaka) and bounds (Ponte Casta˜neda-Willis) are developed to provide the closed-form expressions. In addition, these results are compared to numerical simulations performed based on finite element modelling.


Introduction
Defects have an influence on the macroscopic behaviour of a material, each on a different scale.The overall behaviour of the material can be characterized by its microstructure.This transition from micro-to-macro can be modelled using averaging techniques (homogenization) in order to derive the effective properties of a material.
Homogenization studies often concentrate on the elastic behaviour of a microcracked material.In so-called direct methodology, cracks are represented as material discontinuities with parallel faces.The displacement jumps induced by the cracks allow deriving their contribution to the overall response.For instance, Kachanov [1] and Nemat-Nasser and Hori [2] have provided effective stiffness expressions for arbitrarily oriented microcracks.Eshelby's equivalent inclusion method [3] also offers relevant solutions when considering cracks as flat ellipses (in 2D) or ellipsoids (in 3D).For instance, Mura [4] has studied various ellipsoidal limit cases and Mori and Tanaka [5] have enhanced the representation to the case of multiple interacting inhomogeneities.Note that energy-based bounds developed by Ponte Castañeda and Willis [6] allow accounting for different spatial cracks distribution.
Based on the physical analogy with elasticity (as in pioneering works of Bristow [7], see also [8]), some authors have extended these modelling approaches to thermal, electrical and permeability properties of cracked media * e-mail: helene.welemane@enit.fr[9][10][11].For steady-state heat conduction, Sevostianov [12] and others [13,14] apply the direct methodology based on temperature jump across insulating crack lips.For Hoenig [15], Hatta and Taya [16], Benveniste and Miloh [17] and more recently Shafiro and Kachanov [18], the equivalent inclusion method appears again as a key issue.While several studies account for the arbitrary value of matrix/inclusion conductivity and arbitrary crack's orientation or shape, most of the existing papers generally provide thermal conductivity of microcracked media in the non-interacting case.Nguyen et al. [19] give closedform expression for different schemes but consider only one orientation of the crack.Nevertheless, that is not the only challenge.Opening or closing of microcrack (also known as unilateral effect) can have a different influence on the material, in turn on the overall properties.Consequences of both induced anisotropy and unilateral effect on the elastic problem have been studied by few authors [1,20], but the same cannot be said for the heat conduction problem.
The modelling of the steady-state behaviour within microcracked media can also be achieved through numerical simulation.Carson et al. [21] apply Finite Element Method (FEM) to find the conductivity of non-insulated porous of various shapes and sizes, while Tang et al. [22] propose a similar modelling for concrete with conductive heterogeneities.Shen et al. [23] use a plastic damage model to create cracks under tensile load and then consider steady-state conduction to find the conductivity of the microcracked concrete with high aggregate volume.One can cite also works of Tran et al. [14] based on an adaptive scheme Boundary Element Method (BEM) to find conductivity of a domain containing several cracks.Once again, we note that the crack orientation and unilateral effect are not given enough attention.
The present work intends to propose an Eshelby-like modelling approach for the steady-state heat transfer in a 2D microcracked medium.The effective thermal conductivity is derived based on the geometry of cracks considered as thin aspect ratio inclusions, and on the relevant choice of cracks properties according to their status (open or closed).The theoretical basis of the 2D linear thermal problem is stated in Section 2. As a demonstration, for different estimations (dilute and Mori-Tanaka schemes and Ponte Castañeda-Willis bound), closed-form expressions for a single-family of parallel cracks are presented in Section 3. In addition to the analytical solution, we also propose a numerical analysis of the thermal problem by means of finite element simulations.Modelling and description of the simulated area are given in Section 4. The results obtained from micromechanics and numerical simulation are finally compared and discussed in Section 5.

Theoretical framework
Current works on the effective thermal properties are influenced by the similarities between elasticity and steady-state heat conduction variables [8]:

Elasticity
Heat conduction stress σ heat flux q strain ε temperature gradient g stiffness C thermal conductivity λ Hooke's law Fourier's law Let denote Ω the area of the 2D Representative Volume Element (RVE) of the microcracked media, ∂Ω its outer boundary and u the outward unit normal to ∂Ω (Fig. 1a).The macroscopic temperature gradient G (respectively heat flux Q) can be defined as the mean temperature (resp.external heat flux) on the boundary ∂Ω.Under stationary thermal conditions, the macroscopic temperature gradient G (resp.heat flux Q) corresponds to the average of the corresponding microscopic quantity g (resp.q): where T (x), g(x) and q(x) respectively represent the local temperature, local temperature gradient and local heat flux at any point x of Ω.
The RVE studied here includes two phases (matrix and microcracks), where each phase is homogeneous and locally follow the Fourier's linear thermal law: where λ is the local symmetric second-order thermal conductivity tensor.The matrix is considered to be isotropic, continuous and has the thermal conductivity tensor λ m = λ m I (λ m is the scalar thermal conductivity and I is the second-order identity tensor).This matrix is weakened by randomly distributed single-family of parallel microcracks (Fig. 1a).These cracks are modelled as a flat oblate ellipse (mean semi-axes a and c with c a; Fig. 1b) with unit normal n and 2D volume fraction f c = π d ω.Here d = N a 2 is the scalar crack density (N is the number of cracks per unit area) as defined by Budiansky and O'Connell [24] and ω = c /a 1 is their mean aspect ratio.A uniform macroscopic temperature gradient G is imposed at the outer boundary δΩ of the RVE.Assuming an initial natural state, the microscopic and macroscopic quantities can be linked linearly as [25]: where A is the second-order gradient localization tensor.Similar to (3), the overall behaviour of the RVE can be given by: where λ hom is the overall thermal conductivity of the microcracked media.Assuming the condition A = I and (4), the effective thermal conductivity of the microcracked media comes to: where • r = 1 Ωr Ωr • dΩ denotes the mean value over the area Ω r of the phase r for r = {m, c}.
The single-inhomogeneity problem studied by Eshelby [3] considers a single ellipsoidal inclusion embedded inside an infinite matrix subjected to macroscopic stress or strain tensors at infinity.This elasticity problem can be extended to thermoelasticity [8,26].Here, matrix and cracks exhibit matrix-inclusion topology.In such case, the local temperature gradient in the crack can be approximated by the uniform local field obtained for an ellipsoid embedded in an infinite matrix subjected to uniform boundary conditions G ∞ .
Taking all this into account, the estimated solutions for localization tensor over the crack's phase A c can be determined.And they depend on the depolarization tensor S E (similar to the Eshelby tensor of the elastic problem) whose factors are given in [8].In our case, the said tensor for a flat oblate ellipse can be given as: Note that the configuration of flat cracks corresponds to the limit case where ω → 0, which has to be introduced only at the very end of the mathematical developments.As mentioned before, the unilateral effect is one of the main focus of this work.So, two different results can be obtained at the end, based on the state of the crack (open or closed).In either case, cracks are assumed to be isotropic λ c = λ c I (λ c is the crack's scalar thermal conductivity) with a different value of λ c depending on the state of the crack: -when the cracks are open, λ c = 0, which supports the adiabatic conditions on the crack lips, -when the cracks are closed, λ c = λ * , which accounts for some level of heat transfer continuity; this assumption is inspired by the works of Deudé et al. [20], where closed cracks are represented by a fictitious isotropic material with scalar conductivity λ * = 0.

Calculation of the effective thermal conductivity
We impose a uniform macroscopic thermal gradient G at the outer boundary δΩ of the RVE.This is similar to the classical strain-based formulation in elasticity.As a first, we will estimate the effective conductivity through different schemes and bounds.When there is a dilute concentration of cracks (small d), it is considered that there is no interaction between them.The remote condition in this case can be given by G ∞ = G.Hence, the localization tensor can be given by: Substituting ( 8) in ( 6), we get the general expression: where tensor R is defined as: The above equation is valid for all the mean aspect ratio ω 1 and all the ratio ξ of scalar conductivity between defects and matrix.The present study focuses on the case of flat ellipse-shaped microcracks (c a) for which aspect ratio tends to zero.Besides, we intend to account for different crack status: -open cracks: one has λ c = 0, so ξ = 0 and Accordingly ( 9) can be simplified into: When we are to consider some interactions between cracks, the Mori-Tanaka scheme may provide an interesting solution [5].The boundary condition here is given by G ∞ = g m and the localization tensor reads: This leads to: As before, the specific behaviour of flat cracks according to their status gives the following: Ponte Castañeda-Willis developed an energy-based upper bound to find effective stiffness [6].This bound takes into account the shape of the inclusion (through S E ) and also the spatial distribution of cracks through an additional tensor S d .The simplified localization tensor can be given by: For simplicity, a circular spatial distribution is adopted, for which S d = 1 2 I. Now (6) can be written as: Based on the state of the flat defects, one gets: Note that PCW bound will provide the same result as dilute and MT schemes when no spatial distribution is considered for the former and elliptical distribution for the latter [27].Some main comments can be made regarding these theoretical developments.First, the three modelling approaches show crack induced anisotropy for open cracks.Yet, equations ( 11), ( 14) and ( 17) provide different expressions of the effective conductivity tensor through the tensorial term n ⊗ n.We also observe that as d → 0, all estimations lead to the same result which corresponds to the matrix conductivity λ dil hom , λ M T hom , λ P CW hom → λ m .On the other hand, we observe a complete deactivation of microcracking when the defects are closed λ dil hom = λ M T hom = λ P CW hom = λ m .Note that detailed developments and extension to 3D case can be found in [27].

Numerical simulations
In the following, numerical simulations are performed using finite element software Abaqus®.(t, v) denotes an orthonormal coordinate system.The simulated area A is a square (size L = 1 m) that follows steady-state heat conduction.The matrix is designed as an unit 2D shell with its own scalar conductivity λ m .Assuming there are N = 10 cracks in the area, the radius of the crack can be given by a In what follows, the range of considered density is less than 0.1, related crack's radius has, therefore, a maximum value of 0.1 m.Cracks are usually represented as seams for the open state (duplicated nodes).Yet, this cannot account for the heat transfer during crack closure.So, the crack is modelled here as an elliptical inclusion (created as a partition on the 2D shell) with normal n and scalar conductivity λ c .Since creating a crack with zero aspect ratio is not possible (Ω c = 0), the cracks are designed with an aspect ratio 0 = ω 1 (so f c 1).For a given f c , the value of the scalar conductivity λ c determines if the cracks are open (λ c = 0) or closed (λ c = λ * = 0).Such a description of the crack geometry and the unilateral effect is in line with the theoretical framework used in Section 2.
The cracks are positioned inside the simulated area using circular spatial distribution (in agreement with spatial distribution assumed for PCW bound).Since we want to study the influence of the crack's orientation on conductivity, further simulations are done by rotating the whole group of cracks which maintains a constant distance between them for all orientations.To be precise, the so-called Reference Configuration (RC) corresponds to the distribution of cracks rather grouped near the centre of A to reduce edge effects (Fig. 2).While keeping the circular spatial distribution, other configurations are also studied in the following to show the influence of cracks' location.
The generalized scalar conductivity of a material λ(v) related to the direction of unit vector v is defined by: when the material is subjected to uniform temperature gradient G = G v v.For the numerical model, zero flux condition is imposed on the left and right edges (with outer normals ± t) of area A. At the same time, temperatures T 1 and T 2 (∆T v = T 1 − T 2 > 0) are applied respectively on the top and bottom edges (with outer normals ± v) of the cell; the temperature on each side is uniform.Such latter boundary condition, namely tem- On a global point of view, the two edges with zero flux act as adiabatic walls, allowing the heat flux Q to be mainly oriented along the v direction.From definition (18), the numerical effective conductivity in direction v is then provided by: where Q v is the average heat flux along the v direction.
It can be calculated as being the heat flux density in v direction along the path on the top/bottom edge.Alternatively, Q v can be found using Reaction flux RFL i calculated on each node i on the top/bottom edge when considering unit dimension in the transverse direction, i.e.
The finite element type used for both the matrix and crack is quadratic triangular DC2D6 (see Fig. 3 for RC; (n, v) = 45 • ).Fixing 100 elements inside each crack, the influence of the size of the matrix elements on the heat flux has been studied (Tab.1).We have used a very fine mesh to prevent improper scattering of the flux around the crack tips.In that case, the model has approximately 73500 elements and 148000 nodes including 201 nodes on the outer edges.Note that the computation time remains acceptable (less than a minute).Moreover, estimations of λ num (v) obtained for different cracks distributions have been compared for the most critical case, i.e. (n, v) = 0 • (see Fig. 4).From Table 2, it is observed that cracks   distribution has no major influence on the resulting conductivity.Accordingly, the RC will be considered for all further simulations done in the study.As first illustration, Figure 5 shows the heat flux vector at integration points for density d = 0.1 (RC; (n, v) = 45 • ). Figure 5a corresponds to the open case (λ c = 0) and shows that the cracks acts as a thermal barrier according to the adiabatic behaviour on their lips.Figure 5b corresponds to the closed case (λ c = 50% λ m ) and shows continuity in heat transfer.

Results and discussion
This section intends to compare theoretical developments and FE numerical simulations.From ( 18) and ( 5), the  theoretical scalar conductivity λ(v) comes to: This can be estimated for different schemes (th = {dil, M T, P CW }) and compared to λ num (v).
Recalling previous results from Section 3, open cracks contribute to the degradation of the thermal conductivity, mainly along the direction n normal to the crack surface.This case is true for the simulations as well (see Fig. 6).Both the theoretical and simulated results show us damage-induced anisotropy irrespective of the scheme or crack density.As pointed out earlier, for the theoretical models, we see that as d → 0, λ dil hom ≈ λ M T hom ≈ λ P CW hom (d = 0.1 in Fig. 6a, d = 0.05 in Fig. 6b).This can be attributed to the fact that as d decreases, the size of the crack decreases (respectively from a = 0.1 m to a = 0.07 m), making the interaction between the cracks less influential and at one point there is no interaction between the cracks essentially leading to a dilute configuration.We also see that as the crack becomes smaller, so does its influence on the conductivity (λ(n) ≈ 0.73 λ m for a = 0.1 m whereas λ(n) ≈ 0.86 λ m for a = 0.07 m). Figure 6   the consistency between the theoretical and simulated results.It is interesting to observe that for lesser angles (n, v) < 45 • , simulated results tend towards PCW and for higher angles, they approach the dilute case.Indeed, interactions are greater when the cracks are mostly orthogonal to the heat flux.But, if cracks tend to be aligned with the direction of the temperature gradient then the influence of cracks decreases and heat flux is less disturbed, tending to the dilute case (see Fig. 7).

also illustrates
On the other hand, dilute, Mori-Tanaka and PCW approaches show that closed cracks do not contribute to the degradation of conductivity (see (11), ( 14) and ( 17)), i.e. the effective conductivity in any direction is recovered to its initial value at the cracks' closure.So the generalized scalar conductivity in unit direction v for closed cracks can be given as: λ(v) = λ m , ∀ v. Just like the open crack, simulated and theoretical results are consistent for the closed crack (see Fig. 8).We also see that the former has only a negligible amount of degradation of thermal conductivity (less than 0.035% for d = 0.1 and less than 0.02% for d = 0.05 when considering λ * = 50% λ m ).From ( 11), ( 14) and (17) we know that the theoretical results are not a function of the aspect ratio ω since they all correspond to the limit case ω → 0 λ dil hom , λ M T hom and λ P CW hom only depend on λ m , d and n .But as discussed earlier, it is not possible to simulate an ellipse with zero aspect ratio.So it seems natural to study the influence of the aspect ratio on the simulated results.Since the maximum degradation is along the direction n normal to the crack, we intend to focus only on λ(n).Figure 9a corresponds to open crack and Figure 9b corresponds to closed crack with fixed values of dilute, MT and PCW denoted as reference.In both cases, the simulated results are really sensitive to the aspect ratio ω and get closer to the PCW bound when ω → 0. Especially in the closed case, the simulations tend to the full recovery of λ(n), same as the theoretical models.Note that all the simulations linked to varying aspect ratio are performed by varying the crack thickness c and keeping the crack density d and radius a as constants.Also, the theoretical results for the closed case do not depend on the fictitious scalar conductivity λ * .This may not be true for the simulations.So, a series of simulations were performed with varying λ * and for different aspect ratios (d and a are still constants).The values for λ * are given as a proportion of λ m such that λ * = α λ m with α = {1, 5, 10, 25, 50, 80, 100}[%].Figure 10 shows that there is a clear influence of the scalar conductivity λ * on the numerical thermal conductivity.For α ≤ 10%, we see a drastic decrease in the conductivity, this is due to the fact that we are slowly approaching the open case (α = 0).We also observe that as ω → 0 the influence of λ * diminishes and representation of closed cracks by means of an ellipse with fictitious scalar conductivity λ * becomes independent of the λ * value, just like the theoretical results.

Conclusion and perspectives
In this work, we have presented the closed-form expressions for the effective thermal conductivity of 2D microcracked media under the steady-state heat condition.The theoretical background is based on the equivalent inclusion method where cracks are represented as the limit case of thin elliptic inclusions.Special attention has been paid to the unilateral effect by considering specific properties of the cracks according to their state.Different estimations of the overall thermal behaviour that take into account (Mori-Tanaka scheme and PCW bound) or not (dilute) cracks interactions have been provided, both for open and closed cracks.The determination of effective thermal conductivity was also performed by means of finite element simulations.For this numerical part, geometry and properties of cracks were taken into account in the same way as for theory.The consistency of theoretical and numerical results have been demonstrated through following points.For open cracks, we observe that the microcracked medium exhibits an induced anisotropy with main degradation of conductivity in the direction normal to the cracks.Also, as crack density tends to zero, all models recover property.On the other hand, cracks closure to a complete deactivation of damaging effects.Finally the sensitivity of numerical results based on the aspect ratio of defects has been shown.
Further studies could now be conducted to extend such analytic procedure (homogenization and numerical simulations) for flux-based boundary condition from which effective thermal resistivity can be derived.Considering earlier work of the authors [27], it would also be relevant to compare theoretical and numerical results in the 3D case.

Fig. 2 .
Fig. 2. Simulated area showing spatial distribution of cracks for the Reference Configuration (RC).

Table 1 .
Influence of matrix element size on the average heat flux Q v (RC; (n, v) = 0 • ).Size of the matrix elements 0

Fig. 6 .
Fig. 6.Generalized thermal conductivity λ(v) normalized by its initial value for a material weakened by a single array of parallel open microcracks of unit normal n.

Fig. 8 .
Fig. 8. Generalized thermal conductivity λ(v) normalized by its initial value for a material weakened by a single array of parallel closed microcracks of unit normal n (RC; λ * = 50% λm).

A 2 N
Second-order temperature gradient localization tensor I Second-order identity tensor λ c , λ m Thermal conductivity tensor of the crack and the matrix respectively λ hom Effective thermal conductivity tensor of the microcracked media λ c , λ m Scalar thermal conductivity of the crack and the matrix respectively, W.m −1 .K −1 n Unit vector normal to the crack's plane S E Depolarization tensor A Simulated area with dimension L × L, m Number of cracks per unit area ω Crack's mean aspect ratio Ω, Ω r Area of the RVE (with boundary ∂Ω) and the phase r respectively (t, v) Orthonormal coordinate system of the simulated area a, c Crack's mean semi-axes, m d Scalar crack density f c Crack's volume fraction in 2D g, G Microscopic and macroscopic temperature gradient respectively q, Q Microscopic and macroscopic heat flux respectively ∆T v = T 1 − T 2 Difference between temperatures T 1 and T 2 along the direction v, K u Outward unit normal to ∂Ω ξ Ratio of scalar conductivities between crack and matrix G t , G v Temperature gradient along the direction t and v respectively, K.m −1 Q t , Q v Average heat flux along the direction t and v respectively, W.m −2 T (x), g(x), q(x) Temperature, temperature gradient and heat flux respectively at point x HFL Heat flux in a structure RC Reference Configuration RFL Reaction flux RVE Representative Volume Element