Numerical Modelling of the Transient Hygroscopic Behavior of Flax-Epoxy Composite

This contribution deals with the development of a three-node triangular plane finite element to analyze the transient hygroscopic behavior of 2/2 twill flax fabric-reinforced epoxy composite. Several plates of this material were fabricated using the vacuum infusion process and composite specimens were then cut and aged in tap water at room temperature until saturation. To simplify, a plane modelling of water diffusion in the aged specimens is adopted and Fick’s model is used to describe the water diffusion kinetics. To highlight the heterogeneity of the flax-epoxy samples, the twill flax fabrics waviness is modelled with a sinusoidal undulation. In particular, we show that the proposed finite element formulation allows estimating the flax fiber radial diffusion coefficient by an inverse approach.


Introduction
Over the past two decades, the environmental requirements have encouraged several research and development works to find an ecological alternative to traditional oil-based composite materials, in particular those reinforced with glass fibers. Within this context, natural fiber-reinforced polymer composites and especially those reinforced with flax fibers have gained ground in several industrial areas such as automotive, aeronautic, sports, packaging and construction industries [1][2][3]. In fact, natural fibers combine multiple advantages such as their interesting specific mechanical properties and especially their availability and biodegradability which constitute a major guarantee with respect to environment constraints [4,5]. However, their hydrophilic nature remains the major drawback to their development in semi-structural and structural applications exposed to a wet environment. Indeed, several research works have already shown that the exposure of natural fiber-reinforced composites to humidity induces an important variation of their mechanical and dynamical properties (see [6][7][8][9][10][11][12][13][14][15] among others). These variations are principally related to the swelling of natural fibers due to moisture absorption, which causes different degradations particularly at the fiber-matrix interface and consequently reduces the sustainability of these composite materials. Therefore, the control of these variations passes by the understanding of diffusion kinetics within natural fiber-reinforced polymer composites. Within this context, many research works have focused on these aspects by analyzing the behavior of natural fiber-reinforced composites in wet environment using one-dimensional (1D) Fick's model [8,[16][17][18]. However, this approximation cannot be sufficient in some cases due to the anisotropic nature of the natural fiber, which makes the water diffusion process more complex. Therefore, it remains necessary to take into account the threedimensional (3D) effects in order to better understand water diffusion kinetics in natural fiber reinforced composite materials. To this end, 3D Fick's model has been considered in some works [19,20] to identify the 3D moisture diffusion coefficients. For example, Chilali et al. [20] showed that water diffusion kinetics predicted by 3D Fick's model is in good agreement with the experimental curves. Moreover, they found through the study of several geometric parameters that water diffusion in their studied flax composites is strongly influenced by the morphology of the flax fiber.
Although very practical, these analytical models permit only to study composite materials at their macroscopic scale and for simple geometrical shapes. It is for example difficult to find an analytical solution of 3D Fick's law for an anisotropic material as discussed in [20]. In addition, the analytical models do not allow to predict the diffusion behavior of the composite constituents like the diffusion properties of natural fibers which are relatively difficult to measure experimentally. In this case, the use of a numerical model could outweigh these limitations as it allows for instance to estimate the moisture content at any point of a composite structure with complex shape.
In the literature, several works have dealt with the finite element modelling of moisture diffusion particularly in synthetic fiber-reinforced polymer composites [21][22][23][24][25][26][27]. In these latter papers, finite element simulations have been conducted at the microscopic level of the composite materials to take into account their heterogeneity. For example, Joliff et al. [23] used a 2D finite element model to simulate water diffusion in unidirectional glass fiber-reinforced epoxy composites aged at 70°C in deionized water up to 16 weeks. In particular, they modelled the fiber-matrix interphase and studied the influence of glass fibers distribution with and without barrier effects on the diffusion kinetics of the glass-epoxy composite. Joliff et al. [23] showed that the matrix diffusion coefficient must be higher than that of the bulk epoxy resin to correctly fit the experimental water uptake curve which indicates that the epoxy diffusivity behavior is highly modified by the presence of glass fibers. Contrary to synthetic fiber-reinforced composites, the number of contributions dealing with numerical modelling of moisture diffusion in natural fiberreinforced composites remains limited. Among these contributions, we cite the work by Regazzi et al. [27] which concerns numerical modelling of moisture diffusion in short flax-PLA composites aged in water at different temperatures. In this latter contribution, the studied composite material was assumed homogeneous to simplify the finite element analysis.
In this work, we present the formulation of a three-node triangular membrane finite element to analyze the transient hygroscopic behavior of 2/2 twill flax fabric-reinforced epoxy composite aged in tap water at room temperature. For this purpose, Fick's model is considered to describe the diffusion kinetics within the flax-epoxy composite specimens and their heterogeneity is also taken into consideration by modelling the twill flax fabrics waviness with a sinusoidal undulation.
The present paper is structured as follows. After the introduction, we present in Section 2 the experimental water uptake curve of the flax-epoxy composite. Section 3 is devoted to the finite element approximation of the hygroscopic problem and presents the formulation of the hygroscopic three-node membrane element. Before the concluding remarks, we discuss in Section 4 the obtained numerical results.

Experimental Water Absorption Curve
In this work, the hygroscopic behavior of 2/2 twill flax fabric-reinforced epoxy composite is investigated. The flax fabrics were provided by Depestele group and present an aerial weight of 330 g/m 2 and a fiber density of 1450 kg/m 3 . The epoxy resin is the SR 8100 and was provided with its SD 88225 hardener by Sicomin. Several flax-epoxy laminate plates of dimensions 200 × 200 × 3 mm 3 were prepared by the vacuum infusion process (see [20,28] for more details about the fabrication process). To obtain a thickness of 3 mm, four layers of flax fabrics were used. After that, the flax-epoxy plates were cut and shaped in square specimens of side 20 mm. Then, the obtained specimens were totally immersed into tap water at room temperature to accelerate their water absorption kinetics. During their ageing, the weight variation of each composite specimen was measured and the moisture absorption rate was determined by the following expression: where W0 is the dry initial weight of each specimen and Wt is its weight at time t. The experimental water uptake curve of the flax-epoxy specimens is shown in Fig. 1. Initially, the water absorbed by these samples increases quasi-linearly with the ageing time and then slows down before reaching the equilibrium after 60 days of ageing. This diffusivity behavior can be described by Fick's model [29].
The macroscopic diffusion coefficients D1, D2 and D3 of the flax-epoxy specimens are summarized in Tab. 1 where 1 is the warp direction, 2 is the weft direction and 3 is the thickness direction. These parameters were identified by minimizing the quadratic error between the analytical solution of 3D Fick's model and the experimental results of Fig. 1. In particular, we remark that water diffusion occurs principally across the thickness of the flax-epoxy samples as D3 is found approximately 3 times the plane diffusion coefficients D1 and D2. This result shows that water diffuses more easily across the thickness of the composite specimens than towards their plane directions.

Finite Element Approximation of the Hygroscopic Problem
We consider a 2D composite solid B0 as shown in Fig. 2. To describe its hygroscopic behavior, we introduce the variable ( , ) c t x which represents the moisture content at time t of one point P inside B0 characterized by its position vector

Figure 2: A 2D composite solid with moisture content boundary conditions
We recall that Fick's model is considered to describe water diffusion kinetics in the flax-epoxy specimens: where c is the moisture concentration and D is the diffusion tensor which is supposed symmetric.
For plane problems, D and the gradient of the moisture concentration c ∇ are given by: where x D and y D are the diffusion coefficients along the x and y directions, respectively, while xy D represents the diffusion rate in the x-direction due to a moisture concentration gradient in the y-direction.
The weak form of the diffusion Eq. (2) is obtained by introducing a test function * c that verifies * 0 c = on Sc: By remarking that . 0 where n φ is the normal diffusion flux and 0 B ∂ is the boundary of the 2D composite domain B0.
The finite element approximation of the weak form (6) necessitates spatial and time discretizations. To this end, the heterogeneous plane composite domain B0 is firstly divided into three-node plane elements for spatial discretization and secondly, the backward Euler integration scheme is adopted for time integration.
The proposed hygroscopic plane element is shown in Fig. 3. The following interpolation = associated with the classical three-node membrane element are used to approximate the moisture content within the hygroscopic plane element in terms of the nodal variables: where i c are the nodal moisture contents.
Using the finite element approximation (7), it is possible to show that c and c ∇ are related to the nodal degrees of freedom vector cn as At the element level, the finite element approximation of the weak form (6) reads This approximation could be rewritten as 0 .   (11) The elementary tangent stiffness matrix of the proposed three-node membrane element is obtained by differentiating the residual vector R with respect to n ∆c as: (12)

Numerical Results
In this section, we use the hygroscopic three-node triangular element to study the water absorption behavior of the aged flax-epoxy specimens. In particular, this study permits to obtain an approximation of the flax fiber radial diffusion coefficient and its evolution with thickness. To take account of the heterogeneity of the flax-epoxy samples, we need to describe the 2/2 twill weave of the flax fabrics as detailed in the next section.

Description of the Flax Yarn Undulation
As previously announced, the three-node hygroscopic triangle element is used in the following sections to reproduce the experimental water absorption curve of the flax-epoxy specimens. To simplify, a plane modelling of the diffusion problem is considered as depicted in Fig. 4(b).
To highlight the heterogeneity of the flax-epoxy samples, a sinusoidal model is considered to describe the 2/2 twill weave of the flax fabrics. The 3 mm thick flax-epoxy specimens are constituted of four layers of flax fabrics. Owing to symmetry, only one-fourth of the cross section of each composite specimen is modelled as shown in Fig. 4(b) This corresponds to a unit cell of the 2/2 twill flax fabric. Moreover, the twill flax fabrics used in this study are balanced and accordingly the warp and weft directions (directions x and y in Fig. 4(a) are considered identical.
The geometrical model used to describe the flax yarn undulation is based on sinusoidal functions. In the warp direction (direction x), the flax strand undulation within the modelled 10 1.5 × mm 2 plane section can be approximated by the following function [ where a and h are defined in Fig. 4(c) and correspond to the width and thickness of the flax fiber strand, respectively ( 2.5 a = mm in our case). The geometry of the flax fiber strand is then approximated by the following sinusoidal function: This function is used to calculate the cross section of the flax fiber strand as: The cross Section S and the perimeter of the warp flax strand allow then to calculate the total surface of the flax fiber reinforcement that respects the fiber volume fraction of the flax-epoxy specimens (32%).

Hygroscopic Behavior of the Flax-Epoxy Composite
To model water diffusion within the heterogeneous flax-epoxy composite specimens ( Fig. 4(b)), we need the diffusion coefficients of the epoxy resin and the flax fiber. The diffusion coefficient of the epoxy resin can be deduced from its experimental water absorption curve. It presents a Fickian behavior and the minimization of the quadratic error between the analytical solution of Fick's model and the experimental results gives an isotropic diffusion behavior characterized by we recall that it is difficult to measure its radial and longitudinal diffusion coefficients experimentally. As water diffusion occurs principally across the thickness of the flax-epoxy specimens as shown in Tab. 1, we neglect in this study the longitudinal diffusion coefficient of the flax fiber in order to facilitate the estimation of its radial coefficient designated by Df in the following. Consequently, we use 0 and for the weft direction as depicted in Fig. 4(b).
To obtain a first approximation of the flax fiber radial diffusion coefficient, we use the homogenization model of Halpin-Tsai [31] largely considered in the case of synthetic fiber-reinforced composites: epoxy epoxy epoxy 1 with and 1 1 where Dt is the transverse (across the thickness) diffusion parameter of the flax-epoxy specimens, Df is the flax fiber radial diffusion coefficient and Vf is the flax fiber volume fraction. It is also important to note that ζ is a measure of the reinforcement geometry and is taken equal to 1 in this work to simplify. ). Fig. 6 shows a comparison between the experimental and numerical absorption curves of the flaxepoxy specimens. It is important to note that the numerical water absorption curve of Fig. 6 is the arithmetic average of all nodal moisture concentrations of the finite element model. We remark that the numerical water absorption curve does not accurately describe the experimental results. The difference between these two curves reflects the difficulty to predict the flax-epoxy diffusion behavior using a heterogeneous model.   × mm 2 /h also allows to accurately describe the experimental curve. These results clearly show that the water diffusion coefficient of the epoxy resin (and even that of the flax fiber) should be adjusted from that of the pure epoxy resin to correctly predict the diffusivity behavior of the flaxepoxy specimens. This tendency has been already reported in the literature in the case of synthetic fiberreinforced composites [24,32]. In Joliff et al. [24] for example, the diffusion coefficient of the epoxy matrix was found at minimum 30% larger than that of the bulk epoxy resin in order to fit the experimental water uptake curve of unidirectional glass fiber-reinforced epoxy specimens.
This difference in the value of Depoxy (+ 22% in our case) could be explained by the fact that much more defects appear in the flax-epoxy specimens when compared with the pure epoxy resin samples and they are related firstly to the presence of porosities in the flax-epoxy specimens (the porosity content is approximately equal to 5%) and secondly to the microcracks principally located at the fiber-matrix interface (Fig. 8(a)). In addition to that, the flax fiber is known to be much more hydrophilic than the epoxy resin which induces a differential swelling between the reinforcement and the matrix and leads to more microcracks at the fiber-matrix interface (Fig. 8(a)). Moreover, the flax fiber morphology and microstructure offer many interfaces which increase water diffusion in the flax-epoxy specimens (Figs.  8(b) and 8(c)). All these defects make the epoxy matrix less resistant to water diffusion than the pure epoxy resin and explains in part the difference in Depoxy. For the flax fiber, the results of Fig. 7 permit to obtain a more accurate approximation of its radial diffusion coefficient than that deduced from the Halpin-Tsai model (Eq. (16)).

Effect of Thickness
In this last section, we study the effect of thickness on the diffusion behavior of the flax-epoxy specimens. The thickness effect has been already experimentally investigated by Chilali et al. [20] [20]. Fig. 9 shows the experimental water uptake curves of the sealed flax-epoxy specimens and Tab. 2 summarizes their diffusion coefficients across the thickness and their maximum water uptakes. These transverse diffusion coefficients were identified by minimizing the quadratic error between the analytical solution of 1D Fick's model and the experimental results.  The results of Fig. 9 and Tab. 2 show a drop of water diffusion kinetics when increasing thickness characterized by a decrease of the maximum water uptake and a drop of the water absorption rate defined as the slope of the linear part of the water absorption curve. This tendency has been already reported in the literature in the case of synthetic fiber-reinforced polymer composites [33][34][35]. Bunsell [33] has related this phenomenon to a molecular rearrangement of the polymeric network when specimens become thicker which slows down their diffusion kinetics. Moreover, this variation of diffusion kinetics could also be related to the increasing number of epoxy resin layers enclosing the flax fabrics which may decrease water absorption across the thickness.
Contrary to the maximum water uptake and the water absorption rate, the transverse water diffusion coefficient of the flax-epoxy specimens is found to increase with thickness which shows that the flax fiber and the epoxy resin diffusion coefficients are also sensitive to thickness change. Accordingly, we propose in the following to quantify this evolution by considering the same inverse approach as Section 4.2. Fig. 10 shows a plane modelling of water diffusion in the sealed flax-epoxy specimens. The sinusoidal model introduced in Section 4.1 is always considered to describe the 2/2 twill weave of the flax fabrics. Owing to symmetry, only one-fourth of the cross section of each flax-epoxy specimen is modelled as shown in Fig. 10 Concerning the finite element calculations, we used for each flax-epoxy specimen a refined mesh like that of Fig. 5 after a mesh convergence study. The moisture content boundary conditions were applied only to the top edge of the 2D model as depicted in Fig. 10(a) ( c M ∞ = ). To carry out finite element simulations, we need the flax fiber and the epoxy resin diffusion coefficients. The application of the Halpin-Tsai model (Eq. (16)) allows to obtain a first approximation of these coefficients based on the transverse diffusion coefficients of Tab. 2. After that, the three combinations presented in Section 4.2 were applied to fit the experimental water absorption curves of the flax-epoxy specimens. We depict in Fig. 11 the evolution of the flax fiber and the epoxy resin diffusion coefficients that correctly fit the experimental water uptake curves of Fig. 9. We remark that these diffusion coefficients vary exponentially with respect to thickness with a more pronouncing evolution of the flax fiber radial coefficient. In fact, the epoxy resin diffusion coefficient of the 10 mm thick specimen is approximately 6 times greater than that of the 3 mm thick specimen while the flax fiber radial diffusion coefficient is found to increase by 13-fold from 3 mm to 10 mm thick specimens. This increase of Df and especially that of Depoxy seems to be in accordance with that of the flax-epoxy specimens (D3 increases by 548% between the 3 mm and the 10 mm thick specimens as shown in Tab. 2). The high increase of Df when compared with Depoxy could be related to the hydrophilicity of the flax fabrics when compared with the epoxy matrix which could justify in part the more pronouncing evolution of Df.

Conclusion
In this paper, a three-node membrane finite element was developed to analyze the transient hygroscopic behavior of 2/2 twill flax fabric-reinforced epoxy composite aged in tap water at room temperature until saturation. To take into account the heterogeneity of the flax-epoxy composite in the finite element model, a sinusoidal description of the flax fabrics undulation was adopted and Fick's model was considered to describe the flax-epoxy water diffusion kinetics.
This study firstly allowed to estimate the flax fiber radial diffusion coefficient using an inverse procedure. It is worthy to recall that the flax fiber and more generally natural fibers diffusivity properties are very difficult to obtain experimentally. Secondly, we showed that the diffusion behavior of the epoxy resin is highly modified by the presence of flax fabrics. Finally, the numerical investigation of the thickness effect allowed to follow the exponential evolution of the flax fiber and the epoxy resin diffusion coefficients across the thickness.
As a future work, it would be interesting to estimate the flax fiber longitudinal diffusion coefficient using the same inverse approach. This coefficient was neglected in this study to facilitate the identification of the flax fiber radial diffusion coefficient as water diffusion was found to occur principally across the thickness of the aged flax-epoxy specimens. Besides, it would also be interesting to take account of the dependency of the flax fiber diffusion coefficients on the local moisture content.