<figure><img src="image_1.png" /></figure>

HAL
open science

# Computational Model of Effective Thermal Conductivity of Green Insulating Fibrous Media

Hamidou Sankara, Dominique Baillis, Ousmane Coulibaly, Remi Coquard,

Naïm Naouar, Zahia Saghrouni

► To cite this version:

Hamidou Sankara, Dominique Baillis, Ousmane Coulibaly, Remi Coquard, Naïm Naouar, et al.. Computational Model of Effective Thermal Conductivity of Green Insulating Fibrous Media. Materials, 2024, 17 (1), 10.3390/ma17010252. [hal-04662315](https://hal.science/hal-04662315v1)

HAL Id: hal-04662315

https://hal.science/hal-04662315v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Computational Model of Effective Thermal Conductivity of Green Insulating Fibrous Media

Hamidou Sankara<sup>1,2,*</sup>, Dominique Baillis<sup>1</sup>, Ousmane Coulibaly<sup>2</sup>, Rémi Coquard<sup>3</sup>, Naïm Naouar<sup>1</sup>
and Zahia Saghrouni<sup>4</sup>

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Sankara, H.; Baillis, D.; Coulibaly, O.; Coquard, R.; Naouar, N.; Saghrouni, Z. Computational Model of Effective Thermal Conductivity of Green Insulating Fibrous Media. *Materials* **2024**, *17*, 252. https://doi.org/10.3390/ma17010252

Received: 29 October 2023
Revised: 20 December 2023
Accepted: 29 December 2023
Published: 3 January 2024

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2024 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** Modelling effective thermal properties is crucial for optimizing the thermal performance of materials such as new green insulating fibrous media. In this study, a numerical model is proposed to calculate the effective thermal conductivity of these materials. The fibers are considered to be non-overlapping and randomly oriented in space. The numerical model is based on the finite element method. Particular attention is paid to the accuracy of the results and the influence of the choice of the representative elementary volume (REV) for calculation (cubic or rectangular parallelepiped slice). The calculated effective thermal conductivity of fibrous media under different boundary conditions is also investigated. A set of usual mixed boundary conditions is considered, alongside the uniform temperature gradient conditions. The REV rectangular slice and uniform temperature gradient boundary conditions provide a more accurate estimate of the effective thermal conductivity and are therefore recommended for use in place of the usual cubic representative elementary volume and the usual mixed boundary conditions. This robust model represents a principal novelty of the work. The results are compared with experimental and analytical data previously obtained in the literature for juncus maritimus fibrous media, for different fiber volume fractions, with small relative deviations of 7%. Analytical laws are generally based on simplified assumptions such as infinitely long fibers, and may neglect heat transfer between different phases. Both short and long fiber cases are considered in numerical calculations.

**Keywords:** effective thermal conductivity; thermal insulation; numerical model; homogenization; fibrous media

## 1. Introduction

During the last decade, bio-sourced thermal insulating materials have experienced a strong development, especially concerning building insulating applications. Examples include works on cellulose-based aerogel [1,2] and composite materials using granular cork [3,4] or natural fibers such as hemp shives [5–8], flax and straw rape [8], palm date fibers [9–12], and juncus maritimus fibers [13]. One of the most important factors that make composites increasingly attractive is the ability to manage properties by constitutional design [14]. In order to optimize by the design thermal performances of materials, knowledge and modeling of the thermal properties is of primary importance. Most previous research in the literature has focused on the characterization of the thermal properties of these materials, via analytical and experimental methods. The application of the current paper is particularly concerned with green fibrous media. Concerning this material, the recent work of Saghrouni et al. [13] showed that the analytical Glicksman law was the most

<footnote n="1">¹ LaMCoS, INSA-Lyon, CNRS UMR 5259, Université de Lyon, 69621 Villeurbanne, France; dominique.baillis@insa-lyon.fr (D.B.); naim.naouar@insa-lyon.fr (N.N.)</footnote>

<footnote n="2">² Laboratoire de Physique et de Chimie de l'Environnement (ED-ST/LPCE), Université Joseph KI-ZERBO, Ouagadougou 03 BP 7021, Burkina Faso; coulous2005@yahoo.fr</footnote>

<footnote n="3">³ EC2 Modélisation Campus Lyon Tech, 69603 Villeurbanne, France; remi.coquard@ec2-modelisation.fr</footnote>

<footnote n="4">⁴ Laboratory of Thermal and Energetic Systems Studies (LESTE), National Engineering School of Monastir, University of Monastir, 5019 Monastir, Tunisia; zahiasagrouni@yahoo.fr</footnote>

<footnote>* Correspondence: sankarah461@gmail.com; Tel.: +226-74-53-46-97</footnote>
appropriate by comparison with experimental data. Analytical models are often limited because they are based on approximations. The Glicksman's model is based on a resistor approach and assumes that the fibers are infinitely long and randomly oriented in the surroundings. Therefore, accurate and robust numerical models are required in order to compare with simplified analytical models. This is one of the purposes of the current paper. We propose a numerical model to calculate the effective thermal conductivity of fibrous media and to compare it with Glicksman's law. Numerical models present an advantage because they are more generalist. With the advancement of computer computational capabilities, numerical analyses have been employed for the computation of heat conduction, for example [15–19]. Li et al. [20] determined both the in- and out-of-plane thermal conductivities of composites using the representative volume element (REV) technique with two-unit cells created at varied scales and periodic boundary conditions. However, as heterogeneous media can have prohibitively large representative elementary volume (REV) sizes due to the randomness of the microstructure (non periodic), it is often necessary to estimate the true effective thermal conductivity (ETC) based on the apparent thermal conductivity (ATC) of computational domains (elementary volume) smaller than RVE, for which different boundary conditions may provide different results [21]. Thus, in this current paper, particular attention is placed on the precision of the numerical results and the influence of the choice and size of the computational representative elementary volume (REV). We propose considering two types of elementary volume: usual cubic or slice. Until now, elementary cubic volumes [21] have generally been considered, and there has been virtually no previous work on the influence of REV volume shape. This is an innovative aspect of this article, which aims to show the influence and importance of the choice of volume shape. The calculated effective thermal conductivity of fibers media under different boundary conditions is also investigated. A set of usual mixed boundary conditions is considered within the computational homogenization framework, alongside the uniform temperature gradient conditions. In the literature, temperatures are usually imposed on opposite sides, as in the hot-guarded plate experiment [22]. Recently, another type of boundary condition was proposed, that of gradients (UTGC) [21], but applied to foams. These boundary conditions have not been applied to fibrous materials. We propose studying these two boundary conditions in this article, which is also an innovative point in this paper. The numerical model is based on the finite element method. The fibers were studied in a vacuum, air, and finally in a mortar. The cases of short and long fibers are considered. Numerical results are compared with experimental and Glicksman analytical laws results previously used in literature for juncus maritimus fibers [13]. Deviations in numerical results are quantified for different fiber volume fractions.

## 2. Theoretical Models for Effective Thermal Conductivity

Heat conduction through heterogeneous media depends on the structure of the materials and the thermal conductivity of each phase. Different methods, including analytical, approximate, and numerical approached, allow for an evaluation of the effective thermal conductivity ($k_{eff}$) based on the assumption of the equivalent homogeneous medium.

### 2.1. Analytical Approach

Recent works [13] dealing with the thermal conductivity of insulating juncus maritimus fibrous mortar composites model tend to show that the analytical law proposed by Schuetz and Glicksman [23] is the most appropriate in comparison with the experimental data. The analytical model proposed by Schuetz and Glicksman [23,24] is based on a serial parallel approach applicable to fibers in fluid phase and open-cell foams (Equation (1)). They assumed infinitely long, nearly randomly oriented fibers. They considered the total length of all fibers per unit volume, which were oriented in some direction within an infinitesimal solid angle, to be constant. This is their definition of “random fibers”. Isotherms are assumed to be planes perpendicular to the heat flux. This means planes parallel to the face of elementary volume at imposed temperatures. The Glicksman model assumes that
there is no exchange between the fibers and the surrounding medium. Thus, when the surrounding gas is in a vacuum, the Glicksman model can be used as reference. However, there are approximations when the surrounding medium is not in a vacuum (such as in air or mortar). The effective thermal conductivity of the medium is given by the following:

$$k_{eff} = (1 - f_v)k_{fluid} + f_v \frac{1}{3}k_{fib} \quad (1)$$

where $k_{fluid}$ is the thermal conductivity of the fluid phase, $k_{fib}$ the thermal conductivity of the solid phase, and $f_v$ is the volume fraction of fibers.

Bruggeman's model [25] used Maxwell's model for cylindrical particles (Equation (2)). He obtained an expression for the effective thermal conductivity in the following form:

$$k_{eff} = \frac{k_{mat} \left[ 1 - \left( 1 - \frac{k_{fib}}{k_{mat}} \right) \frac{2}{3} f_v \delta \right]}{[1 + (\delta - 1) f_v]} \quad (2)$$

where $k_{mat}$ is the thermal conductivity of the continuous phase, $k_{fib}$ is the thermal conductivity of the dispersed phase and $f_v$ is the volume fraction of the dispersed phase, respectively. $\delta$ is a parameter that depends on the shape of the particles: spherical, cylindrical, or flat.

For cylindrical particles (fibers): $\delta = \frac{5k_{mat}+k_{fib}}{3(k_{fib}+k_{mat})}$.

## 2.2. Numerical Calculation of the Thermal Conductivity of Fibrous Composites by FEM

The basic idea of computational homogenization techniques is to obtain the effective properties of heterogeneous materials by solving the Fourier heat equation in a representative elementary volume (REV) of the composite materials with appropriate boundary conditions. The simplest way to determine the thermal conductivity of fibrous from numerical approaches consists of solving the steady state energy equation (Fourier law) in a representative sample subjected to a thermal gradient. As a reminder, the Fourier's law (Equation (3)) is written as follows:

$$\vec{\phi} = -\lambda \overrightarrow{grad} T \quad (3)$$

$T$ represents the temperature and $\lambda$ the thermal conductivity of the material. In our approach, a temperature difference is imposed between two sides of a parallelepipedal sample. The numerical simulation allows for the computation of the average flux crossing the sample, i.e., both sides with imposed temperatures, in a direction normal to these faces and in stationary regime. The effective conductivity of the heterogeneous material is then obtained using the Fourier relation:

$$k_{eff} = \frac{L Q}{S |\Delta T|} \quad (4)$$

where $\Delta T/L$ is the applied thermal gradient in which $L$ is the thickness of the sample and $Q$ refers to the heat flux through the sample over an area $S$, which is normal to $\Delta T/L$. When using (Equation (4)), the following points should be highlighted:

* First, one-dimensional (1D) heat transfer is considered; therefore, the area $S$ must be much larger than the thickness $L$ of the sample. To circumvent this constraint, it is common practice to apply adiabatic conditions on lateral boundaries where temperatures are not prescribed (i.e., boundaries parallel to the direction of heat flow).

* Secondly, this formula assumes that a linear variation in temperature prevails through the thickness of the sample. To achieve this, the absolute value of $\Delta T$ must be chosen to be sufficiently small.

* (Equation (4)) allows for estimating the effective conductivity in the direction of the imposed temperature gradient. Therefore, in order to evaluate the possible anisotropy
of the material thermal behavior, the same calculation must be performed in each of the three Cartesian directions.

* Finally, we also consider that the thermal contacts between the matrix material and the fibers are perfect: there is no thermal contact resistance at the fiber/matrix interfaces.

The main unknown in (Equation (4)) is the heat flux Q. Q is closely related to the conductivity and to the spatial arrangement of the different phases of the porous material, which govern the heat diffusion paths inside the heterogeneous material. Q also depends on the dimensions of REV and on the boundary conditions applied. The solution of the heat equation by means of the finite-element method is a powerful tool whose formulation is particularly well adapted to the estimation of this flux Q. This numerical method is implemented in numerous common software. In the present work, we used the ABAQUS FEM software.

## 2.3. FEM Numerical Modeling with ABAQUS Software

The ABAQUS finite-element software have a thermal solver, allowing for solving the heat equation (Fourier law equation (3)) in transient or stationary regimes, notably inside composite materials. The software is also equipped with an integrated interface for model development and visualization (ABAQUS CAE), which allows for generating and meshing physical models with the prescribed boundary conditions. The model generation could be automatized by means of Python scripts. Such automatization scripts were used in the present work. They allowed us to generate parallelepipedal representative elementary volumes (REVs) composed of arrangements of fibers immersed in a continuous matrix. The fibers were modeled as circular cylinders whose length-to-diameter ratio ($\chi = L_f/D_f$) could be varied. Moreover, the distribution of the fiber orientation was also monitored. Two different kinds of fiber orientation distribution were considered in the present work:

* Random orientation in the space ("random 3D");

* Random orientation in a plane ("random 2D")—this plane may be, for example, parallel to some sides of the parallelepipedal REV.

Additionally, the fiber volume fraction $f_v$ was also controlled. Finally, the dimensions of the REV (length, width, height of $L_x-L_y-L_z$, respectively) could be modified.

The algorithm that generated and meshed the REV was based on the following process:

* A first step, conducted before the recourse to the ABAQUS utilization, consisted of generating the arrangement of fibers with the aimed characteristics (fiber volume fraction $f_v$, 2D (Figure 1) or 3D random orientations, and diameter distribution $DIS(D_f)$). The algorithm developed in Fortran language, unfolded as follows:

1. Random location, orientation, diameter, and length of the new fiber are randomly chosen, i.e., by means of successive random numbers $\zeta$ comprised between 0 and 1; $1 < \zeta < 0$:

* $x_n = \zeta \cdot L_x; y_n = \zeta \cdot L_y; z_n = \zeta \cdot L_z$

* $(\theta_n, \varphi_n, \gamma_n)$ for a random 3D distribution, we have $\theta_n = \cos^{-1}(\zeta)$; $\varphi_n = 2 \cdot \pi \cdot \zeta$ and $\gamma_n = 0$; for a random 2D distribution of fibers in a plane normal to the Z-axis, we have $\theta_n = \pi/2$; $\varphi_n = 2 \cdot \pi \cdot \zeta$ and $\gamma_n = 0$

* $D_n$ is obtained from a random number $\zeta$ in order to satisfy:

$$ \zeta = \frac{\int_{0}^{D_n} DIS(D_f) \cdot dD_f}{\int_{0}^{D_n} DIS(D_f) \cdot dD_f} \quad (5) $$

* $L_n = \chi \cdot D_n$

2. Then, we check if this newly generated fiber intersects a previously generated fiber of the current arrangement. In this case, the process returns to Step 1.
3. Thereafter, the current fiber volume fraction $f_{v,c}$ (fraction of the REV $L_x \cdot L_y \cdot L_z$ occupied by the fibers) is updated by adding the volume of the new fibers comprised in the REV.

<figure><img src="image_5.png" /><figcaption>Figure 1. 2D visualization of fibers randomly oriented in space.</figcaption></figure>

Steps 1–3 are repeated until the volume fraction of the arrangement reaches the desired value ($f_{v,c} \ge f_v$).

Figures 2 and 3 show a cubic REV and a slice REV with randomly oriented fibers, respectively. Figure 4 illustrates the orientation of the fibers in the matrix, where $e_x, e_y$, and $e_z$ are the local unit vectors of the Cartesian coordinate system. They can be expressed using the spherical coordinates, as follows:

$$ \begin{cases} x = r\sin\theta\cos\varphi \\ y = r\sin\theta\sin\varphi \\ z = r\cos\theta \end{cases} \qquad (6) $$

<figure><img src="image_6.png" /><figcaption>Figure 2. infinite randomly oriented fibers in a parallelepipedal REV ($L_x = L_y \approx L_z$).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** infinite randomly oriented fibers in a slice REV ($L_x = L_y \gg L_z$).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Fiber orientation coordinate system.</figcaption></figure>

* Once the locations and orientations of the fibers in the arrangement with the provided characteristics have been determined, the generation process utilizes the ABAQUS integrated interface (ABAQUS CAE) in order to construct the physical model on which the thermal simulation will be conducted. Thus, several automated steps are successively performed in ABAQUS CAE using Python scripts:

* Step 1: Geometry generation

Firstly, a volume containing the N cylindrical objects (each of them with a diameter and length $D_n$, $L_n$ for $n = 1, N$) representing the N fibers of the arrangement is generated. These volumes are created from a homogeneous medium noted A. Thereafter, a second parallelepipedal volume (B) with the dimensions of the REV (length–width–height of $L_x-L_y-L_z$, respectively) containing a homogeneous medium B representing the matrix is generated. A boolean operation (implemented in ABAQUS CAE) is then applied to remove the N cylinders (Volume A) from the parallelepipedal volume (B) to create Volume C. Then, Volumes A and C are gathered together inside a physical model composed of two phases (Volume
A: fibers; Volume C: matrix). Finally, still in ABAQUS CAE, this physical model is cut according to the planes $x = 0; x = L_x; y = 0; y = L_y; z = 0$ and $z = L_z$ in order to comply with the dimensions of the parallelepipedal REV.

### Step 2: Meshing of the model

The meshing tool integrated in the CAE interface is used to generate finite elements of Volumes C (Matrix) and A (N fibers). Quadratic tetrahedron elements («DC3D10» in ABAQUS) illustrated in Figure 5c are used in the present work. They allow for an accurate estimation of the thermal fields with a reasonable mesh refinement. We checked that the mesh densities used were sufficient in order for the results to be independent of the refinement.

### Step 3: Thermal contact conditions assignment

As explained previously, the thermal contact at the interface between the matrix and fibers is considered perfect. To account for this hypothesis, in ABAQUS CAE, a thermal continuity condition is specified at the interface between the two volumes (Volume A/Volume C) with a very large thermal conductance leading to a negligible thermal contact resistance. Note that a non-perfect thermal contact between the two phases might be modeled by specifying a lower value of the thermal conductance. The fibers can come into contact with the rest of fibers.

### Step 4: Material Properties Assignment

In ABAQUS CAE, the thermal conductivities are affected at each phase ($k_{mat}$ for Volume C representing the Matrix; $k_{fib}$ for Volumes A representing the N fibers) and might depend on the temperature.

The ABAQUS model is then complete and ready for the realization of the thermal simulations according to prescribed boundary conditions. The different successive steps for the model generation are illustrated in the following Figure 6.

<figure><img src="image_9.png" /><figcaption>Figure 5. (a) Mixed boundary condition (MBC); (b) uniform thermal gradient condition (UTGC); (c) quadratic mesh.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6. Cont.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Illustration of the successive steps for the model generation in ABAQUS CAE.</figcaption></figure>

## 2.4. Boundary Conditions Used for Modeling

As indicated previously, the boundary conditions applied to the FEM computations could influence the value of the computed effective conductivity. The two different types of boundary conditions considered in this work are presented in this section. Thus, we use the so-called “mixed” boundary conditions (MBCs) and the “uniform thermal gradient” conditions (UTGC). These two types of boundary conditions have been commonly encountered in previous works on porous materials [26–29] and sporadically used in computer homogenization [30,31] in order to analyze the macro-homogeneity condition. For MBC (Figure 5a), a hot temperature on one side and a cold temperature on the opposite side are imposed (Equation (8)), while the other four sides are considered adiabatic by the ABAQUS software (Equation (9)):

$$ \frac{\partial}{\partial x} \left( k_{x,y,z} \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y} \left( k_{x,y,z} \frac{\partial T}{\partial y} \right) + \frac{\partial}{\partial z} \left( k_{x,y,z} \frac{\partial T}{\partial z} \right) = 0 \quad (7) $$

$$ T|_{z=0} = T_{\text{cold}}; \quad T|_{z=Z} = T_{\text{hot}} \quad (8) $$

$$ \left. \frac{\partial T}{\partial x} \right|_{x=0} = \left. \frac{\partial T}{\partial x} \right|_{x=X} = 0; \quad \left. \frac{\partial T}{\partial y} \right|_{y=0} = \left. \frac{\partial T}{\partial y} \right|_{y=Y} = 0 \quad (9) $$

For the UTGC (Figure 5b), the hot temperature on one side and the cold temperature on the opposite side are still prescribed, but on the other lateral sides, we impose a linear temperature variation, which is expressed by the following equations:

$$ T|_{z=0} = T_{\text{cold}}; \quad T|_{z=Z} = T_{\text{hot}} \quad (10) $$
$$ T|_{x=0} = T|_{x=X} = T_{\text{hot}} + \left( \frac{T_{\text{hot}} - T_{\text{cold}}}{e} \right) z \quad (11) $$

$$ T|_{y=0} = T|_{y=Y} = T_{\text{hot}} + \left( \frac{T_{\text{hot}} - T_{\text{cold}}}{e} \right) z \quad (12) $$

with Z being the thickness between hot and cold faces.

## 2.5. Estimation of the Effective Conductivity from the FEM Computations

The FEM simulations conducted with both types of boundaryconditions (UTGC or MBC) allow for evaluating the homogeneous effective conductivity of the composite material. Thus, the conduction heat flux density crossing REV can be evaluated numerically from the temperature field in a steady-state regime. In the ABAQUS software, this resulting flux density is directly available for each node at which the temperature is prescribed through the RFL11 variable (so-called «reaction Flux»). The effective thermal conductivity ($k_{eff}$) is then evaluated using (Equation (13)):

$$ k_{eff} = \frac{\left| \sum RFL11 \right| \cdot e}{S \cdot \Delta T} \quad (13) $$

where $\sum RFL11$ is the total heat flux density crossing the REV from the hot to the cold faces, S is the section of the faces, e is the thickness of REV according to the temperature gradient, and $\Delta T$ is the temperature difference between the two faces. This method was notably successfully used by Coquard et al. [32].

# 3. Results and Discussion

For the present simulations, let us consider a hot temperature $T_{hot} = 20$ °C on the one side and a cold temperature $T_{cold} = 0$ °C on the opposite side. The choice of hot and cold temperatures was actually quite arbitrary, and had no influence on the results in terms of thermal conductivity, as the conductivities of the constituent materials (fiber and matrix) used in our calculations were always independent of temperature. As the size of REVs must be large enough to statistically represent the microstructure, and the computation time increases with REV size, it is important to study the influence of REV size to ensure good convergence towards a fair and accurate result, and to determine an appropriate REV size. We considered that the fibers were randomly distributed in the matrix according to a uniform distribution (i.e., there was no preferential orientation). Finally, we approximated the shape of the fibers using cylinders of variable diameter and length, as well as observing a uniform distribution law.

## 3.1. Infinitely Long Fibers in Vacuum

### 3.1.1. Influence of Mesh Size

In order to select the appropriate mesh size for further simulations, we performed a mesh sensitivity study. We considered a cubic REV with edge 6, a target volume fraction of 4% and 15%, with an infinite fiber length (i.e., $1/a = 5000$) and a constant fiber diameter of 0.2. We then compared the computed results with the Glicksman relation. The parametric mesh study in a vacuum was conducted considering a solid thermal conductivity with a value of $(1 \text{ W.K}^{-1}.\text{m}^{-1})$ and a fluid thermal conductivity $(10^{-6} \text{ W.K}^{-1}.\text{m}^{-1})$. In the rest of the article, the value of conductivity $(10^{-6} \text{ W.K}^{-1}.\text{m}^{-1})$ is taken as the vacuum matrix conductivity. So, the finite element mesh construction from our binary image was as follows: (i) a geometric definition of the fluid and solid phases by triangular facets, and (ii) a tetrahedral mesh of the volume. We used four-node linear elements, the ABAQUS DC3D10 heat transfer elements [31]. The meshing step was performed using the open-source code Iso2mesh [33].

Figures 7 and 8 show the weak influence of the mesh size. For each of the two boundary conditions (MBC or UTGC), when the mesh size increased, the value of thermal conductivity was practically unchanged. The mesh size was sufficiently small to calculate
stable and precise results. We chose a mesh size of 0.1 for the rest of our simulations, which constituted a good compromise between the simulation accuracy and calculation time. We observed that the results of the gradient (UTGC) and mixed (MBC) boundary condition computations did not converge to the same value, probably because the REV was not large enough.

<figure><img src="image_12.png" /><figcaption>Figure 7. Influence of mesh size for a volume fraction of 4%.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8. Influence of mesh size for a volume fraction of 15%.</figcaption></figure>

### 3.1.2. Influence of the Cubic REV Size

To obtain more accurate results, we performed calculations in all three directions (X, Y, and Z) through averaging—running the same simulation 10 times to obtain an average value for each parameter.

As the fluid is in a vacuum, Glicksman's thermal conductivity value will tend towards 1/3 of the fibers' volume fraction. It can be noted that the Glicksman model is appropriate for fibers in a vacuum. Thus, it could be considered as a reference in this case. We can see that the MBC computation tended to converge with a deviation of 75% from the Glicksman model (Figure 9). The volume did not seem to be large enough to have a good convergence, similarly to the conclusions of Z.K Low et al. [21] in the case of foams. UGTC seemed
to converge faster but with oscillations towards the Glicksman model as the REV volume increased. These results underline the need to generate larger REVs to converge for more stable and accurate results, as well as the limitations of considering a cubic REV. In the remainder of our study, non-cubic but parallelepiped REVs (small thickness and large surface area, i.e., in slices) will be examined.

<figure><img src="image_14.png" /><figcaption>Figure 9. Influence of cubic REV size at 4% volume fraction.</figcaption></figure>

### 3.1.3. Influence of the Slice REV Size

In this section, we consider a thin slice $L_z$ with characteristic surface dimensions ($L_x$, $L_y$ and $L_z$) with $L_x \approx L_y \gg L_z$. $L_z$ is assumed constant, $L_z = 1$ (Figures 10 and 11).

<figure><img src="image_15.png" /><figcaption>Figure 10. Geometry and meshing of a 3D slice ($L_x = L_y = 10$ and $L_z = 1$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11. 3D fiber temperature fields ($L_x=L_y = 10$ and $L_z = 1$).</figcaption></figure>
Figures 12 and 13 show the influence of slice REV size ($L_x = L_y$) for a volume fraction of 4% and 15%. For a fiber volume fraction of 4%, when the lengths $L_x \approx L_y$ are very close to $L_z$, significant fluctuations in effective thermal conductivity ($k_{eff}$) are observed, as was the case for a cubic REV. The curve tends to stabilize when $L_x \approx L_y$ is greater than 20. For UGTC, for a REV size greater than 20, the deviations are very small, varying between 0% and 5%. UGTC converges faster. Figure 13 shows the influence of the REV size for a fiber volume fraction of 15%, the UGTC model converges better than the MBC model. We can retain that in the case of a slice, the Glicksman model and the numerical model are in very good agreement for $L_x \approx L_y$ greater than 20. Next, we consider a slice REV with $L_x \approx L_y = 20$.

<figure><img src="image_17.png" /><figcaption>Figure 12. Influence of REV size at 4% volume fraction.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 13. Influence of REV size at 15% volume fraction.</figcaption></figure>

### 3.1.4. Influence of Volume Fraction

Fiber volume fraction is an important parameter, as this value can influence the material's thermal and mechanical properties.

Figure 14 shows the variation in ETC as a function of the volume fraction ranging from 4% to 15%. Firstly, the effective thermal conductivity (ETC) of the material increases with the volume fraction. As expected, numerical and analytical results show that thermal conductivity depends linearly on the volume fraction. The effective thermal conductivity (ETC) of the material increases with the volume fraction, because the conductivity of fibers is higher than the air conductivity. It can be observed that there are very small deviations between the values obtained numerically and those calculated with the Glicksman model (deviations of the order of 0.39%). This confirms that the choice of an element size of 0.1 and slice REV with $L_x = L_y = 20$ are appropriate, inducing precise results.
<figure><img src="image_19.png" /><figcaption>Figure 14. Influence of volume fraction.</figcaption></figure>

## 3.2. *Random Fibers in Air: Infinitely Long and Short Fibers*

Glicksman's model is based on the assumption of infinitely long fibers, so it is important to quantify the difference of calculated conductivity between infinitely long fibers and real fibers (not infinitely long). In the case of infinitely long fibers, they extend from one face to another face of parallelepiped REV. In the case of short fibers, there are more embedded fibers in the volume, so there's more exchange with the surrounding environment. In this section, we considered long and short fibers with lengths of 1000 and 1, respectively, and diameters of 0.2. The fibers were considered in air with volume fractions of 4% and 15%. The conductivity of the air matrix was 0.026 W.K⁻¹.m⁻¹. Simulations were performed with a thermal conductivity of fiber bulk material $k_{fib} = 1$ W.K⁻¹.m⁻¹ in Figures 15 and 16, as well as considering the thermal conductivity of juntus maritimus fibers ($k_{fib} = 0.472$ W.K⁻¹.m⁻¹) on Figures 17 and 18. Figures 15–18 show the evolution of the effective thermal conductivity as a function of REV size. MBC and UTGC seemed to converge to the same value for increasing REV sizes. The UTGC converged better. In Figures 17 and 18, discrepancies exist with the Glicksman model that remain small, less than 10%. These discrepancies can be explained by exchanges between the air and fibers (especially for $k_{fib} = 0.472$ W.K⁻¹.m⁻¹) which are not taken into account in the Glicksman model. For the rest of our study, we chose REVs with dimensions of $L_x \approx L_y = 20$ and $L_z = 1$. We note that the deviations between short (SF) and long fibers (LF) results remained low.
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Comparison of numerical and analytical models for increasing the size of REV fibers in air, $k_{fib} = 1 \text{ W.K}^{-1}.\text{m}^{-1}$ with a volume fraction of 4%.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 16.** Comparison of numerical and analytical models for increasing the size of REV fibers in air, $k_{fib} = 1 \text{ W.K}^{-1}.\text{m}^{-1}$ with a volume fraction of 15%.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 17.** Comparison of numerical and analytical models for increasing the size of REV fibers in air, $k_{fib} = 0.472 \text{ W.K}^{-1}.\text{m}^{-1}$ with a volume fraction of 4%.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Comparison of numerical and analytical models for increasing the size of REV fibers in air, $k_{fib} = 0.472 \text{ W.K}^{-1}.\text{m}^{-1}$ with a volume fraction of 15%.</figcaption></figure>

### 3.3. Fibers in Mortar: Case of Infinitely Long and Short Fibers

This part of the study consists of investigating the numerical model for composite materials and comparing it with previous works on juncus maritimus fibers randomly arranged in the mortar [13]. As the matrix porosity increases with the fiber volume fraction, the thermal conductivity of the matrix needs to be adjusted for each of the fiber volume fractions. For this purpose, we used the experimental values presented in Table 1 for the matrix thermal conductivity [13]. A slice REV of dimension $L_x \approx L_y = 20$ was considered.
<table><caption>**Table 1.** Physical properties of composite samples; JMC1, JMC2, JMC3, JMC4, are the juncus maritimus/mortar samples corresponding to the juncus maritimus mass fractions (2%, 5%, 7%, and 10%), respectively [13].</caption><thead><tr><th>Sample</th><th>JMC1</th><th>JMC2</th><th>JMC3</th><th>JMC4</th></tr></thead><tbody><tr><td>Mass ratio of juncus maritimus to mortar mixture (%)</td><td>2</td><td>5</td><td>7</td><td>10</td></tr><tr><td>Juncus maritimus volume fraction (%)</td><td>3.25</td><td>6.62</td><td>8.86</td><td>11.06</td></tr><tr><td>Mortar thermal conductivity (W.K<sup>-1</sup>.m<sup>-1</sup>)</td><td>0.454</td><td>0.239</td><td>0.211</td><td>0.182</td></tr><tr><td>Composite Mortar/JM Exp (W.K<sup>-1</sup>.m<sup>-1</sup>)</td><td>0.44</td><td>0.236</td><td>0.197</td><td>0.171</td></tr><tr><td>Juncus maritimus Thermal Conductivity (W.K<sup>-1</sup>.m<sup>-1</sup>)</td><td colspan="4">0.472</td></tr></tbody></table>

Table 1 shows the thermophysical properties of the various composite samples obtained from the work of Saghrouni et al. [13]. These properties were useful as input data in the numerical simulations. Note that the local thermal conductivity of juncus maritimus was particularly difficult to determine. It was identified in Saghrouni et al. [13] by assuming that Glicksman's law could be used and by fitting experimental effective conductivities of samples using analytical methods. Thus, it explains the very good agreement between experimental and Glicksman results (Figure 19). Note that Glicksman's relation was used by replacing the thermal conductivity of the fluid with that of the mortar matrix in Equation (1).

<figure><img src="image_24.png" /><figcaption>**Figure 19.** Comparison of the effective thermal conductivities of JM/mortar composites obtained experimentally and numerically for short and long fibers with the Glicksman model calculated using the estimate $k_{JM,fib} = 0.472 \text{ W.K}^{-1}\text{.m}^{-1}$.</figcaption></figure>

The obtained numerical thermal conductivity results were compared with the experimentally [13] and analytically obtained values in Figure 19. First of all, we can see that the numerical, experimental, and analytical results followed a non-linear decreasing curve. This was due to the thermal conductivity of the mortar varying in a decreasing manner with the increased fiber volume fraction (Table 1). Indeed, the matrix porosity increased with the increased fiber volume fraction. A small discrepancy between the numerical results (MBC, UGTC, and short and long fibers) was observed. The relative differences observed between the numerical models and the Glicksman model were weak, ranging from 0.5% to
7%. These results confirmed that Glicksman's law can be used as a good approximation, as proposed in [13].

## 4. Conclusions

A new robust numerical approach to calculate the effective conductivity of insulating fibrous media is proposed and compared analytically. Fibers are considered as being randomly oriented. This study provides a useful tool for evaluating the effective thermal conductivity of a desired volume fraction of fibers, and quantifies the deviation from the simplified analytical models. This paper highlights the importance of using two types of boundary conditions, MBC and UGTC, to verify the convergence of the numerical model. It shows that the UGTC boundary conditions model converged faster than the MBC model. The numerical results also show that a slice-type REV was more appropriate for evaluating the effective conductivity than a cubic REV, which can generate errors due to the very large REV size required. The proposed numerical model can be compared with analytical models such as the Glicksman model. The results are compared with experimental and analytical data previously obtained in the literature for juncus maritimus fibrous media for different fiber volume fractions. There are no numerical results on JM's thermal conductivity in the literature. The analytical model is assumed to be valid, so deviations with an accurate numerical model have to be quantified. The deviations between numerical results for short and long fibers are weak. This paper tends to confirm that the Glicksman analytical law could be used as it has a good approximation with a deviation from numerical results lower than 10%.

**Author Contributions:** Methodology, H.S., D.B., O.C., R.C., N.N. and Z.S.; Writing—original draft, H.S.; Writing—review & editing, H.S., O.C. and D.B.; Investigation, D.B. and N.N.; Conceptualization, D.B. and O.C.; Supervision, D.B., O.C. and N.N.; Validation, O.C., R.C. and Z.S.; Data curation, O.C.; Software, R.C. and N.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** We declared that there is no conflict of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

<table><tr><td>ATC</td><td>Apparent thermal conductivity</td></tr><tr><td>BC</td><td>Boundary condition</td></tr><tr><td>ETC</td><td>Effective thermal conductivity</td></tr><tr><td>Exp</td><td>Experimental</td></tr><tr><td>FEM</td><td>Finite element method</td></tr><tr><td>JMC</td><td>Juncus maritimus composite</td></tr><tr><td>LF</td><td>Long fiber</td></tr><tr><td>MBC</td><td>Mixed boundary conditions</td></tr><tr><td>RD</td><td>Relative deviation</td></tr><tr><td>REV</td><td>Representative element volume</td></tr><tr><td>SF</td><td>Short fiber</td></tr><tr><td>UTGC</td><td>Uniform temperature gradient conditions</td></tr></table>

### Acronyms

* *e* Thickness (m)

* *k* Thermal conductivity (W.K⁻¹.m⁻¹)

* *Q* Heat flow (W)

* *S* Surface (m²)
<table><tr><td>T</td><td>Temperature (°C)</td></tr><tr><td colspan="2"><strong>Subscripts</strong></td></tr><tr><td>eff</td><td>Effective</td></tr><tr><td>fib</td><td>Fiber</td></tr><tr><td>fluid</td><td>Fluid phase</td></tr><tr><td>fv</td><td>Volume fraction</td></tr><tr><td>JM</td><td>Juncus maritimus</td></tr><tr><td>mat</td><td>Matrix</td></tr><tr><td colspan="2"><strong>Greek symbols</strong></td></tr><tr><td>χ</td><td>fiber length to diameter ratio</td></tr><tr><td>Δ</td><td>Difference</td></tr><tr><td>δ</td><td>cylindrical particles (fibers)</td></tr></table>

## References

1. Baillis, D.; Coquard, R.; Moura, L. Heat transfer in cellulose-based aerogels: Analytical modelling and measurements. *Energy* **2015**, *84*, 732–744. [CrossRef]

2. Coquard, R.; Baillis, D. Thermal conductivity of Kelvin cell cellulosic aerogels: Analytical and Monte Carlo approaches. *J. Mater. Sci.* **2017**, *52*, 11135–11145. [CrossRef]

3. Cherki, A.B.; Remy, B.; Khabbazi, A.; Jannot, Y.; Baillis, D. Experimental thermal properties characterization of insulating cork-gypsum composite. *Constr. Build. Mater.* **2014**, *54*, 202–209. [CrossRef]

4. Barreca, F.; Fichera, C.R. Thermal insulation performance assessment of agglomerated cork boards. *Wood Fiber Sci.* **2016**, *48*, 96–103.

5. Collet, F.; Pretot, S. Thermal conductivity of hemp concretes: Variation with formulation, density and water content. *Constr. Build. Mater.* **2014**, *65*, 612–619. [CrossRef]

6. Evrard, A.; De Herde, A. Hygrothermal performance of lime-hemp wall assemblies. *J. Build. Phys.* **2010**, *34*, 5–25. [CrossRef]

7. Le, A.T.; Maalouf, C.; Mai, T.H.; Wurtz, E.; Collet, F. Transient hygrothermal behaviour of a hemp concrete building envelope. *Energy Build.* **2010**, *42*, 1797–1806.

8. Rahim, M.; Douzane, O.; Le, A.T.; Langlet, T. Effect of moisture and temperature on thermal properties of three bio-based materials. *Constr. Build. Mater.* **2016**, *111*, 119–127. [CrossRef]

9. Benmansour, N.; Agoudjil, B.; Gherabli, A.; Kareche, A.; Boudenne, A. Thermal and mechanical performance of natural mortar reinforced with date palm fibers for use as insulating materials in building. *Energy Build.* **2014**, *81*, 98–104. [CrossRef]

10. Chikhi, M.; Agoudjil, B.; Boudenne, A.; Gherabli, A. Experimental investigation of new biocomposite with low cost for thermal insulation. *Energy Build.* **2013**, *66*, 267–273. [CrossRef]

11. Agoudjil, B.; Benchabane, A.; Boudenne, A.; Ibos, L.; Fois, M. Renewable materials to reduce building heat loss: Characterization of date palm wood. *Energy Build.* **2011**, *43*, 491–497. [CrossRef]

12. Boumhaout, M.; Boukhattem, L.; Hamdi, H.; Benhamou, B.; Nouh, F.A. Thermomechanical characterization of a bio-composite building material: Mortar reinforced with date palm fibers mesh. *Constr. Build. Mater.* **2017**, *135*, 241–250. [CrossRef]

13. Saghrouni, Z.; Baillis, D.; Naouar, N.; Blal, N.; Jemni, A. Thermal properties of new insulating juncus maritimus fibrous mortar composites/experimental results and analytical laws. *Appl. Sci.* **2019**, *9*, 981. [CrossRef]

14. Xu, Y.; Yagi, K. Automatic FEM model generation for evaluating thermal conductivity of composite with random materials arrangement. *Comput. Mater. Sci.* **2004**, *30*, 242–250. [CrossRef]

15. Penazzi, L. Construction d'une Fonction de Transfert par la Méthode Monte Carlo Symbolique: Application à la Thermique Couplée en Géométries Complexes. Ph.D. Thesis, Ecole des Mines d'Albi-Carmaux, Albi, France, 2020.

16. Khalid, M.Y.; Al Rashid, A.; Arif, Z.U.; Ahmed, W.; Arshad, H.; Zaidi, A.A. Natural fiber reinforced composites: Sustainable materials for emerging applications. *Results Eng.* **2021**, *11*, 100263. [CrossRef]

17. Nurazzi, N.; Asyraf, M.; Khalina, A.; Abdullah, N.; Aisyah, H.; Rafiqah, S.A.; Sabaruddin, F.; Kamarudin, S.; Norrrahim, M.; Ilyas, R.; et al. A review on natural fiber reinforced polymer composite for bullet proof and ballistic applications. *Polymers* **2021**, *13*, 646. [CrossRef]

18. Vijayan, R.; Krishnamoorthy, A. Review on natural fiber reinforced composites. *Mater. Today Proc.* **2019**, *16*, 897–906. [CrossRef]

19. Tan, L.; Shi, X.; Cheng, T.; Zeng, X.; Zheng, H. Numerical study on effective thermal conductivity of transparent electrospun polymer composite. *Appl. Therm. Eng.* **2019**, *160*, 113949. [CrossRef]

20. Li, H.; Li, S.; Wang, Y. Prediction of effective thermal conductivities of woven fabric composites using unit cells at multiple length scales. *J. Mater. Res.* **2011**, *26*, 384–394. [CrossRef]

21. Low, Z.K.; Blal, N.; Naouar, N.; Baillis, D. Influence of boundary conditions on computation of the effective thermal conductivity of foams. *Int. J. Heat Mass Transf.* **2020**, *155*, 119781. [CrossRef]

22. Randrianalsoa, J.; Baillis, D. Thermal conductive and radiative properties of solid foams: Traditional and recent advanced modelling approaches. *Comptes Rendus Phys.* **2014**, *15*, 683–695. [CrossRef]

23. Schuetz, M.; Glicksman, L. A Basic Study of Heat Transfer Through Foam Insulation. *J. Cell. Plast.* **1984**, *20*, 114–121. [CrossRef]

24. Schuetz, M.A. Heat Transfer in Foam Insulation. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 1983.
25. Öchsner, A.; Murch, G.E.; de Lemos, M.J. *Cellular and Porous Materials: Thermal Properties Simulation and Prediction*; John Wiley & Sons: Hoboken, NJ, USA, 2008.

26. Suter, S.; Steinfeld, A.; Haussener, S. Pore-level engineering of macroporous media for increased performance of solar-driven thermochemical fuel processing. *Int. J. Heat Mass Transf.* **2014**, *78*, 688–698. [CrossRef]

27. Petrasch, J.; Schrader, B.; Wyss, P.; Steinfeld, A. Tomography-based determination of the effective thermal conductivity of fluid-saturated reticulate porous ceramics. *J. Heat Transf.* **2008**, *130*, 032602. [CrossRef]

28. Mendes, M.A.; Ray, S.; Trimis, D. A simple and efficient method for the evaluation of effective thermal conductivity of open-cell foam-like structures. *Int. J. Heat Mass Transf.* **2013**, *66*, 412–422. [CrossRef]

29. Baillis, D.; Coquard, R. Radiative and conductive thermal properties of foams. In *Cellular and Porous Materials: Thermal Properties Simulation and Prediction*; John Wiley & Sons: Hoboken, NJ, USA, 2008; pp. 343–384.

30. Dirrenberger, J.; Forest, S.; Jeulin, D. Towards gigantic RVE sizes for 3D stochastic fibrous networks. *Int. J. Solids Struct.* **2014**, *51*, 359–376. [CrossRef]

31. Jiang, M.; Jasiuk, I.; Ostoja-Starzewski, M. Apparent thermal conductivity of periodic two-dimensional composites. *Comput. Mater. Sci.* **2002**, *25*, 329–338. [CrossRef]

32. Coquard, R.; Loretz, M.; Baillis, D. Conductive Heat Transfer in Metallic/Ceramic Open-Cell Foams. *Adv. Eng. Mater.* **2008**, *10*, 323–337. [CrossRef]

33. Zhang, Y.; Fang, Q. BlenderPhotonics: An integrated open-source software environment for three-dimensional meshing and photon simulations in complex tissues. *J. Biomed. Opt.* **2022**, *27*, 083014. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.