© S.R. Rangasamy Mahendren et al., published by EDP Sciences 2020

1 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 [9–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 closed-form 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.

2 Theoretical framework

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

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): (1) (2)

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: (3)

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 conductivityand 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 is the scalar crack density ( is the number of cracks per unit area) as defined by Budiansky and O’Connell [24] and 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]: (4)

where A is the second-order gradient localization tensor. Similar to (3), the overall behaviour of the RVE can be given by: (5)

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: (6)

where = 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 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: (7)

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 thecracks are open, λ_c = 0, which supports the adiabatic conditions on the crack lips,

• when thecracks 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.

Elasticity		Heat conduction
stress σ		heat flux q
strain ε		temperature gradient g
stiffness		thermal conductivity λ
Hooke’s law		Fourier’s law

Fig. 1 (a) RVE with single family of parallel microcracks, (b) crack geometry.

3 Calculation of the effective thermal conductivity

We impose auniform 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 conductivitythrough 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: (8)

Substituting (8) in (6), we get the general expression: (9)

where tensor R is defined as: (10)

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 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 ,

• closed crack: one has λ_c = λ^*≠0, so ξ≠0 and .

Accordingly (9) can be simplified into: (11)

When we are to consider some interactions between cracks, the Mori-Tanaka scheme may provide an interesting solution [5]. The boundary condition here is givenby and the localization tensor reads: (12)

This leads to: (13)

As before, the specific behaviour of flat cracks according to their status gives the following: (14)

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: (15)

For simplicity, a circular spatial distribution is adopted, for which S_d = . Now (6) can be written as: (16)

Based on the state of the flat defects, one gets: (17)

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 . On the other hand, we observe a complete deactivation of microcracking when the defects are closed . Note that detailed developments and extension to 3D case can be found in [27].

4 Numerical simulations

In the following, numerical simulations are performed using finite element software Abaqus®. (t, v) denotes an orthonormal coordinate system. The simulated area 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 cracks in the area, the radius of the crack can be given by . 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 crackclosure. 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 designedwith 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 or closed . 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 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: (18)

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 . At the same time, temperatures T₁ and T₂ (ΔT_v = T₁ − T₂ > 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 temperature gradient , creates an overall heat flux inside the simulated area. 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: (19)

where Q_v is the average heat flux along the v direction. It can be calculated as with HFL2 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. RFL_i.

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 thatcase, 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 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 and shows continuity in heat transfer.

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

Fig. 3 Simulated area mesh .

Fig. 4 Influence of cracks position for (n, v) = 0°.

Fig. 5 Heat flux vectors at integration points inside the simulated area .

5 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: (20)

This can be estimated for different schemes 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, (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 also illustrates 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).

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 givenas: λ(v) = λ_m, ∀ v. Just like theopen 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 and 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.

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. 7 Heat flux vectors at integration points for various orientations of open cracks .

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 .

Fig. 9 Normal thermal conductivity λ(n) normalized by its initial value for various aspect ratios ; points represent simulation results (RC); log10 scale is usedfor abscissa.

Fig. 10 Influence of scalar conductivity λ* on the normal thermal conductivity λ(n) normalized by its initial value ; log10 scale is usedfor abscissa.

6 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 the matrix property. On the other hand, cracks closure leads to a complete deactivationof 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.

Nomenclature

A: 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⁻¹.K⁻¹

n: Unit vector normal to the crack’s plane

S_E: Depolarization tensor

: 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₁ − T₂: Difference between temperatures T₁ and T₂ 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⁻¹

Q_t, Q_v: Average heat flux along the direction t and v respectively, W.m⁻²

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

References

All Tables

All Figures

	Fig. 1 (a) RVE with single family of parallel microcracks, (b) crack geometry.
In the text

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

	Fig. 3 Simulated area mesh .
In the text

	Fig. 4 Influence of cracks position for (n, v) = 0°.
In the text

	Fig. 5 Heat flux vectors at integration points inside the simulated area .
In the text

	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.
In the text

	Fig. 7 Heat flux vectors at integration points for various orientations of open cracks .
In the text

	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 .
In the text

	Fig. 9 Normal thermal conductivity λ(n) normalized by its initial value for various aspect ratios ; points represent simulation results (RC); log10 scale is usedfor abscissa.
In the text

	Fig. 10 Influence of scalar conductivity λ* on the normal thermal conductivity λ(n) normalized by its initial value ; log10 scale is usedfor abscissa.
In the text