Hygroscopic phase field fracture modelling of composite materials

This paper investigates the effect of moisture content upon the degradation behaviour of composite materials. A coupled phase field framework considering moisture diffusion, hygroscopic expansion, and fracture behaviour is developed. This multi-physics framework is used to explore the damage evolution of composite materials, spanning the micro-, meso- and macro-scales. The micro-scale unit-cell model shows how the mismatch between the hygroscopic expansion of fibre and matrix leads to interface debonding. From the meso-scale ply-level model, we learn that the distribution of fibres has a minor influence on the material properties, while increasing moisture content facilitates interface debonding. The macro-scale laminate-level model shows that moisture induces a higher degree of damage on the longitudinal ply relative to the transverse ply. This work opens a new avenue to understand and predict environmentally-assisted degradation in composite materials.


Introduction
Lightweight composite materials have been widely used in aerospace, wind and marine structural applications due to their high strength/stiffness to weight ratios and corrosion resistance.However, composite materials are very sensitive to environmental conditions such as moisture and temperature [1].Moisture contents were found to degrade the mechanical properties of composite materials significantly via fibre-matrix interface debonding [2], hygroscopic expansion [3] and plasticisation of the matrix [4][5][6]; see Fig. 1 for examples of moisture-induced composite material degradation.It has been reported that the absorbed moisture in glass fibre epoxy composites results in up to 60% reduction in tensile strength [7].Hygroscopic expansion of composite materials is known to cause fibre-matrix interface debonding.This issue becomes more evident in natural fibres such as flax fibre, which are intrinsically hydrophilic and very sensitive to the degree of humidity in the environment [8].Moisture contents also decrease the glass transition temperature of composite materials, affecting their long-term durability [9].Therefore, fundamental understanding of moisture-related degradation mechanisms is essential to enable the development of moisture-resistant composite materials.
A large number of experimental studies have been devoted to the measurement of moisture intake and subsequent material degradation [4,5,[10][11][12][13][14]. Buehler and Seferis [4] found that water absorption led to a decrease in flexural modulus/strength and mode I/II interlaminar fracture toughness of glass fibre reinforced plastics (GFRP).Sateesh et al. [5] reported that the flexural modulus of GFRP decreased by 20-25% when subjected to a hygrothermal environment for six months.LeBlanc and LaPlante [14] investigated mixed-mode delamination growth, showing that the delamination toughness decreased by 25-62% across loading modes.Lu et al. [10] reported that the flexural properties of flax-epoxy composites were reduced by up to 20% when the relative humidity increased to 97%.However, there are inherent challenges associated with the experimental quantification of micro-scale damage evolution over long periods, and it is often difficult to gain fundamental insight into material degradation mechanisms under coupled mechanical and hygroscopic environmental conditions from laboratory tests.This has triggered an interest in the development of computational models to understand moisture-assisted degradation in composite materials.Some studies are focused on modelling moisture diffusion in composite materials with stationary pre-cracks or cavities.Sinchuk et al. [15] developed a meso-scale model with Fickian diffusion theory to characterise the moisture diffusion path and local moisture concentration.Gagani and Echtermeyer [16] combined a representative volume element (RVE) model with Fickian diffusion to understand transport in composite laminates containing pre-existing cracks and delamination regions.Bourennane et al. [17] used the same methodology and accounted for the role of cavities in the RVE model.Chilali et al. [2] simulated the hydro-elastic behaviour of flax fabric-reinforced epoxy composites and obtained the stress distribution at the fibre-matrix interface, although interface debonding was not explicitly modelled.Some authors have also considered the role of growing cracks.For example, Wong et al. [18] simulated delamination growth in the presence of moisture using moving-node elements.However, this brings in a degree of mesh sensitivity and requires a complex implementation, as remeshing and an explicit treatment of the boundary conditions at the moving crack faces are needed.Cheng et al. [19] developed a hygro-thermal continuum meso-scale model considering damage evolution.Their model solves the diffusion-stress-damage problem in a decoupled fashion and is inherently mesh-dependent.These limitations can be overcome by the use of variational phase field fracture models.
The phase field fracture model builds upon Griffith's thermodynamics principles [20] and has been widely used to predict the evolution of cracks in a wide range of materials, including rock-like materials [21,22], ceramics [23,24], ductile metals [25,26], functionally graded materials [27,28], porous media [29,30], shape memory materials [31,32], ice [33,34] and composites [35,36].In addition, it has been successfully employed to model fracture and bridging behaviour in fibre-reinforced composite materials at the microscale level [37,38].Moreover, phase field fracture modelling can readily be coupled to equations describing other physical phenomena and thus its use has been particularly popular in addressing multi-physics problems.Examples include hydrogen embrittlement [39,40], Li-Ion battery degradation [41,42], thermo-mechanical fracture [43,44], and stress corrosion cracking [45,46].However, the use of multi-physics phase field-based formulations to investigate environmentally-assisted degradation in composites materials is an area that remains relatively unexplored.Only very recently two studies have been published in this regard.Arash et al. [47] presented a phenomenological phase field-based model to investigate the effect of hydrothermal effects on the failure of the behaviour of nanocomposite materials.In this study, the transport problem was not explicitly resolved and, instead, a uniform moisture distribution was assumed.Ye and Zhang [48] formulated a coupled phase field model that accounted for the interaction of stresses, moisture and crack propagation in composite materials.Their framework incorporated the hygroscopic expansion phenomenon due to moisture diffusion and the viscoelastic behaviour of polymer matrix.They also presented a so-called "crack filter theory" to regularise the sharp fibre-matrix interface for phase field fracture and moisture diffusion.The crack filter controls the moisture fluxes that can diffuse through the crack.Nevertheless, finding the coefficients for the crack filter functions remains a challenge.
In this work, we develop a new phase field-based multi-physics framework for moisture-assisted fracture in composite materials.The framework combines: (i) moisture diffusion, (ii) hygroscopic expansion, (iii) phase field fracture, and (iv) a diffuse interface description of fibre-matrix interface debonding.The potential of the proposed model is demonstrated by addressing three representative benchmark studies spanning the fibre, ply and laminate scales.The results obtained show that the proposed framework can shed new light into the fundamental mechanisms governing environmentally-assisted degradation in composite materials, while also serving as a tool capable of delivering durability predictions over technologically-relevant scales.

Numerical model
In this section, we formulate a new multi-physics phase field model to capture the coupling of moisture diffusion, damage evolution and stress redistribution, as shown in Fig. 2. First, the phase field description of fracture is presented (Section 2.1).Then, in Section 2.2, we proceed to formulate the governing equations for the transport of moisture, including the definition of the hygroscopic strains.A diffused interface theory is subsequently proposed to simulate fibre-matrix debonding (Section 2.3).The coupling of these three main aspects enables a thermodynamically consistent framework that can capture the interplay between mass transport, hygroscopic expansion, crack evolution and stress redistribution.

Phase Field Fracture
We use the phase field model to predict the cracking of fibres, matrix and fibre-matrix interfaces (debonding).The phase field fracture model builds upon Fig. 2: A coupled multi-physics phase field model capturing the interaction between stress, moisture and cracking in the fibre, matrix and fibre-matrix interface.
Griffith's thermodynamics framework [20] -for a crack to propagate, the reduction in potential energy that occurs during crack growth must be balanced with the increase in surface energy resulting from the creation of new free surfaces.In an elastic solid undergoing prescribed displacements, this can be expressed mathematically as follows, where E is the total energy, dA is an incremental increase in the crack area, Ψ = ψ dV is the elastic strain energy, with ψ being the elastic strain energy density, ε e is the elastic strain tensor, and W c is the work required to create two new surfaces.The last term in Eq. ( 1) is typically referred to as the critical energy release rate G c = dW c /dA, a material property that characterises the fracture resistance of the solid.For a solid of domain Ω containing a crack of surface Γ, Griffith's energy balance can be formulated in a variational form as [49]: The evolution of cracks can then be predicted as an exchange between stored and fracture energies.However, minimisation of Eq. ( 2) is computationally challenging due to the unknown nature of Γ.To overcome this, Griffith's functional can be regularised using the phase field paradigm.The crack-solid material interface is no longer sharp and discontinuous but it is instead smeared over a finite domain and tracked using a phase field order parameter φ.The phase field resembles a damage variable, going from φ = 0 in intact material points to φ = 1 inside of cracks.Also, following damage mechanics arguments, the stiffness of the solid is degraded by means of a degradation function g(φ) = (1 − φ) 2 .Thus, the regularised functional reads where we have adopted constitutive choices for the degradation function and the crack density function intrinsic to the so-called AT2 phase field fracture model [50].As rigorously proven using Gamma-convergence, the regularised functional E (3) converges to E (2) for a choice of phase field length scale → 0 + .The presence of a phase field length scale makes the model nonlocal and thus ensures mesh objectivity [51].The strong form of the coupled fracture-deformation problem can be readily obtained by taking the variation of E with respect to ε e and φ, noting that the elastic strain tensor is a function of the displacement vector, ε e (u), and that the stress tensor is given by σ = ∂ψ/∂ε e , and applying Gauss' divergence theorem.Accordingly, the strong form equations read where σ 0 denotes the undamaged stress tensor; σ = g(φ)σ 0 .We adopt a strain energy decomposition to prevent crack growth in compressive stress states.Based on the volumetric-deviatoric split method [52], the volumetric and deviatoric strain energies can be subjected to damage but not the compressive volumetric strain energy.Thus, the elastic strain energy is decomposed as ψ = ψ + + ψ − , with where ε e is the deviatoric part of the elastic strain tensor and x ± = 1/2 (x± x ) are the Macaulay brackets.Also, damage must be irreversible, such that φ(x, t) ≥ 0. To enforce this, we introduce a history field H to store the elastic contribution to the damage driving force.I.e., for a current time t, over a total time t t , the history field can be defined as, and the phase field evolution equation is reformulated as,

Moisture Diffusion
Moisture diffusion is governed by Fick's first law.Thus, for a moisture diffusion coefficient D, the change in moisture concentration C with time t is given by, The driving force for the moisture diffusion is governed by the gradient of the chemical potential µ.The chemical potential of moisture is given by, where R is the gas constant, T is the absolute temperature and µ 0 denotes the chemical potential on the standard state.Then, the mass flux is defined based on a linear Onsager relationship, such that The mass conservation relates the rate of moisture concentration with moisture flux through the external surface.Using the divergence theorem, the strong form of the balance equation can be obtained, For an arbitrary scalar field, δC, the variational form of Eq. ( 11) reads, Using the divergence theorem and rearranging the equation, the weak form becomes, where q = J•n is the concentration flux leaving the solid through a surface ∂Ω q .By inserting the mass flux, Eq. ( 10), into the chemical potential of moisture, Eq. ( 13), the weak form of the moisture diffusion becomes, Changes in time of the moisture concentration lead to the presence of hygroscopic strains ε m .The hygroscopic strain scales linearly with the change of moisture concentration where α is the moisture expansion coefficient and C 0 is the initial moisture concentration in the bulk polymer.Accordingly, the constitutive behaviour of the material incorporates such strains, resulting in the following Cauchy stress definition, where C is the linear elastic stiffness matrix and ε is the total strain tensor.In this regard, it is worth noting that the fibre is assumed to be a transversely isotropic material.

Diffuse interface
The definition of sharp interfaces complicates the modelling of coupled multiphysics moisture-driven fracture problems.Instead, we choose to use a diffuse interface paradigm to interpolate relevant material properties across the fibre-matrix interface and thus capture fibre-matrix debonding and moisture diffusion across the fibre-matrix interface.To this end, we exploit once again the phase field paradigm to define, as an initial step, a transition zone between the interface and the bulk materials.Thus, a phase field indicator parameter is defined as: Similar to the phase field fracture description provided above (see Section 2.1), a length scale parameter d comes into play to govern the width of the interface region.For simplicity, throughout this work we choose to assume the same magnitude for the phase field fracture length scale and the interface indicator parameter length scale d .Also, the phase field order parameter d is used to interpolate between two limit states: d = 0, corresponding to the bulk, and d = 1, denoting the interface; see Fig. 3a.An interpolation function n is used to smoothly transition between the interface, denoted by the superscript I, and the bulk, denoted by the superscript (i).The latter is chosen to highlight that the bulk can correspond to two domains: the polymer matrix, (i) = m, or the fibre, (i) = f .With this in mind, we use this diffuse interface theory to interpolate the material toughness, the diffusion coefficient and the moisture expansion coefficient, Effectively, this interpolation of material properties prior to the analysis enables simulating fibre-matrix interface debonding, as characterised by the fracture energy of this process, G I c , and capturing the role that the different nature of the fibre-matrix interface has on moisture transport and hygroscopic swelling.In this regard, the different diffusion coefficients assigned to the interface and the bulk regions capture the different energy barriers inherent to the transport of water molecules through bi-material interfaces.This paradigm is new and provides a simple avenue for incorporating the characteristics of moisture diffusion.In Appendix A we provide details of the implementation of this diffuse interface approach into the commercial finite element package ABAQUS by means of user subroutines1 .

Finite element implementation
In the following, we describe the finite element (FE) discretisation of the deformation and fracture problems (Section 3.1), followed by the discretisation of the moisture transport equation (Section 3.2), and the description of the coupled problem and its implementation (Section 3.3).Common to these sub-sections, the nodal values of the primary fields are interpolated as follows: where Voigt notation has been adopted, n denotes the number of nodes and N i are the interpolation matrices -diagonal matrices with the nodal shape functions N i as components.Similarly, the corresponding gradient quantities are discretised as follows, where B u i denotes the standard strain-displacement matrices and B i are vectors with the spatial derivatives of the shape functions.

FE discretisation of the deformation-phase field problem
Recall the principle of virtual work and consider the constitutive relations outlined in Section 2. The weak form for the coupled deformation-phase field fracture problem can be formulated as Making use of the FE discretisation given in Eqs. ( 21) and ( 22), and considering that Eq. ( 23) must hold for arbitrary values of the primal field variables, the residuals can be derived as follows: with κ being a sufficiently small numerical parameter introduced to keep the system of equations well-conditioned when φ = 1.We choose to adopt a value of κ = 1×10 −7 .The components of the stiffness matrices can then be obtained by differentiating the residuals with respect to the incremental nodal variables as follows:

FE discretisation of the mass transport problem
The residual vector for the moisture transport problem can be readily obtained by discretising Eq. ( 14), given that δC is an arbitrary variation of the moisture concentration; A diffusivity matrix can then be defined, Together with a concentration capacity matrix, and a diffusion flux vector, Subsequently, the discretised moisture transport equation reads,

Coupled scheme
The deformation, diffusion and phase field fracture problems are weakly coupled.First, mass transport affects the stress field in the fracture process zone via hygroscopic expansion.Secondly, the resulting mechanical fields drive the evolution of the phase field and the reduction in load carrying capacity of the solid.The linearised finite element system,   is solved in an incremental manner, using the Newton-Raphson method.The solution scheme follows a so-called staggered approach [53], in that the solutions for the displacement, moisture concentration and phase field problems are obtained sequentially.The implementation is carried out in the finite element package ABAQUS by means of a user element (UEL) subroutine.As described in Appendix A, the overall framework is implemented in two steps.In the first step, a HETVAL user subroutine is used to introduce the diffuse interface, by generating a scalar field as the interface indicator.In the second step, the physical simulation is carried out, with a UEL user subroutine being used to define the residuals and stiffness matrices defined above and to assign properties based on the phase field indicator defined in the first stage, as per Eq. ( 18).

Results
We proceed to showcase the abilities of the computational framework developed to predict moisture-induced composite material degradation across the scales.Firstly, in Section 4.1, the role of moisture in generating mechanical deformation and damage is investigated using a single fibre model undergoing wetting and drying.Secondly, multiple fibres are considered in a micro-scale model aimed at investigating the role of fibre distribution in moisture diffusion and subsequent material degradation (Section 4.2).Thirdly, in Section 4.3, we turn our attention to the coupling between mechanical load and moisture transport through a micro-scale single-edge cracked model.Fourthly, a meso-scale model is used to understand the distribution of damage at the ply level and compare model predictions with experiments (Section 4.4).Finally, a macroscale model with 8 plies is studied to investigate the role of hygroscopic swelling on laminate-level structural integrity.Common to all these case studies is the choice of materials adopted.We focus on a composite of epoxy matrix and flax fibre, whose properties are given in Table 1, together with those of the fibrematrix interface.Although natural fibres are promising sustainable materials compared to synthetic materials, they are hydrophilic and extremely sensitive to humid environments.The single fibre model depicted in Fig. 4a is first investigated to gain fundamental understanding into interfacial damage due to hygroscopic swelling.The composite dimensions are H = W = 0.02 mm with a fibre diameter of d = 0.01 mm.The vertical displacement is fixed at the top and bottom boundaries, while the horizontal displacement component is fixed at the right boundary.In other words, no mechanical loading is applied and deformation is purely the result of hygroscopic swelling due to moisture gradients.At the left edge, a boundary condition for the moisture concentration is applied.Specifically, a two-stage process is followed, whereby an initial moisture absorption stage is followed by a drying stage.In the first stage the moisture content is prescribed to be equal to 7.45% [2] for 2,000 seconds.This is followed by a second stage where the moisture is reduced to 0% and the calculation runs for further 5,000 seconds.Initially, the moisture concentration is assumed to be C(t = 0) = 0 in the entire domain.As described in Section 2.3, the interface indicator model is used to capture fibre-matrix debonding, with different work of fractures required for rupturing the epoxy matrix, the flax fibre and the interface (see Table 1).The entire domain is discretised using 1,862 8-node quadratic elements, with the characteristic element size being equal to 0.0005 mm, two times smaller than the phase field length scale [53].Fig. 4b shows the evolution of the moisture content in the centre of the domain (centre of the fibre) as a function of time.It can be seen that initially the moisture concentration increases sharply and eventually a saturation or equilibrium regime is attained, which is followed by a drop in the moisture content after 2,000 seconds.The results obtained are shown in Fig. 5. Contours of moisture concentration are provided for t = 100 s in Fig. 5a.The C distribution shows a noticeable gradient, which results in hygroscopic strains.The moisture front is not uniform due to the role that the fibre and the matrix-fibre interface play, due to their lower diffusivities.These gradients in concentration eventually disappear as the plateau stage is reached.The impact of moisture diffusion on deformation and damage is given in Figs.5b-f for a time of t = 2, 000 s. Namely, Fig. 5b shows the displacement magnitude contours, revealing the impact of hygroscopic expansion.The resulting stresses are given in Figs.5d and e, which provides contours of σ xx and σ yy , respectively.It can be seen that the maximum tensile stresses are attained at the interface due to the material property mismatch.The damage distribution is shown in Fig. 5d; damage accumulates at the interface, indicative of fibre-matrix interface debonding.This failure mechanism has been experimentally observed in flax fibre composites [10] -see Fig. 1b.Finally, Fig. 5f shows the evolution of the vertical reaction force as a function of time.Three stages can be observed, that resemble the three moisture content stages shown in Fig. 4b.The load increases as hygroscopic swelling effects increase with moisture diffusion until a more homogeneous C distribution is attained, and the load finally drops during the drying process.Strains and stresses decrease significantly through the drying process but the phase field damage along the interface remains the same due to the irreversibility of damage.

Micro-scale: multi-fibre model
The modelling framework is subsequently used to shed light into the role of fibre distribution on moisture transport and ensuing material degradation, an area that remains unexplored in the literature.To this end, a micro-scale model of a composite material with multiple reinforcing fibres is adopted.Two scenarios are considered: (i) a Square Arrayed (SA) multi-fibre model, where 36 fibres are uniformly distributed over 6 columns and 6 rows (Fig. 6a), and (ii) a Randomly Distributed (RD) multi-fibre model, with the same number of fibres as the SA model but randomly distributed (Fig. 6b).The diffuse interface approach is adopted to introduce the desired properties in all fibre-matrix interfaces and capture multiple fibre-matrix debonding.In both SA and RD scenarios, the composite domain has dimensions of H = W = 0.1 mm and the fibre diameter equals d = 0.01 mm.The boundary conditions, sketched in Fig. 6, involve fixing the vertical displacement in the top and bottom boundaries, and prescribing the horizontal displacement on the right edge.The moisture boundary conditions follow a two-step wet-dry analysis.Thus, for an initially dry sample, a Dirichlet boundary condition C = 7.45% is first adopted for a total of 30,000 s, which is sufficient time to reach the equilibrium state.Then, a drying stage begins, where the boundary condition is switched to C = 0 and the calculation is run for further 60,000 seconds.Fig. 6b shows the evolution of the moisture content at the centre of the domain, revealing three distinct stages that resemble those observed for the single-fibre analysis (Section 4.1).Approximately 60,000 quadratic elements with reduced integration are used to discretise the composite domain, with the characteristic element length being equal to 0.00045 mm, two times smaller than and d .
The results obtained are shown in Fig. 7 in terms of contours of moisture content at t = 2, 000 s (Fig. 7a), phase field damage at t = 30, 000 s (Fig. 7b), and stress component σ xx at t = 30, 000 s (Fig. 7c); for both SA and RD fibre models.The results show similar qualitative trends for both scenarios.The trends are also qualitatively similar to those reported for the single-fibre model (Section 4.1) albeit some degree of fibre to fibre interaction is observed in terms of stress fields and phase field damage.The force versus time results confirm the impression extracted from the qualitative contours; as shown in Fig 8 the mechanical responses of the SA and RD fibre models are almost indistinguishable.However, it should be noted that this conclusion is relevant to the material parameters and conditions assumed here.That is, different conclusions might be drawn if the fibre, matrix and interface exhibited higher differences in diffusivity, or if the role of fibre distribution had been explored beyond (e.g., dissimilar volume fractions and fibre radii).
Following the moisture absorption, a drying stage was applied until t = 60,000 seconds.At the end of the drying stage, the moisture content level decreased from 7.45% to 0.829% in the SA model and 0.835% in the RD model, 0.006% lesser in the SA model.Fig. 8a) shows the result of the vertical reaction force measured from the bottom fixed boundary.During the equilibrium state of moisture, the reaction force of the SA and RD models are 26.04N and 25.78 N, respectively.The reaction force is in a linear relationship with the hygroscopic expansion, hence the SA model obtained a higher reaction force from a higher moisture content.Fig. 8b) shows very similar behaviour in both models during the processes of moisture absorption, equilibrium state, and drying.This is mainly attributed to the relatively uniform distribution of fibre in both models and the close diffusion coefficients of fibre and matrix.It is possible that a non-uniform fibre distribution, for instance, a cluster of fibres, will affect the results significantly.From the above analysis, we can conclude that the uniformly-distributed fibre formation in the matrix will not significantly affect the behaviour of the composite under the hygroscopic expansion effect.

Micro-scale: coupled hygroscopic-mechanical model
A micro-mechanical single-edge cracked plate model is used to explore the influence of moisture diffusion on damage tolerance under the action of both mechanical strains resulting from the application of mechanical load and hygroscopic swelling.The geometry and boundary conditions are given in Fig. 9a).The composite plate has a width of W = 0.1 mm, a height of H = 0.2 mm and a random distribution of fibres with diameter d = 0.01 mm.An initial crack is introduced geometrically with a length of a 0 = W/2 = 0.05 mm.The plate is mechanically loaded by prescribing a vertical displacement at the upper edge, while both vertical and horizontal displacements are prescribed at the bottom boundary.Regarding moisture, three scenarios are considered.In one case (referred to as "No moisture"), C = 0 in the entire domain at all times, and damage is driven only by mechanical loading.A second case study, referred to as "Moisture absorbed", refers to a sample that is initially dry but that is exposed to a C = 7.45% on the left edge for 30,000 s.And a third and final scenario, denoted "Moisture dried", assumes that the sample has a uniform, initial moisture content of C = 7.45% and then is exposed to a C = 0% environment on the left side for 60,000 s (Dirichlet boundary condition).Fig. 9b shows the moisture content at the centre of the domain for the 'Moisture absorbed" and "Moisture dried" case studies.The mechanicalmoisture analysis is done in a sequential fashion; that is, the environmental analysis is carried out first, and then a mechanical load is applied where the remote displacement is ramped at a rate of 0.04 mm/s.The entire domain is discretised using 71,390 8-node quadratic elements, with the characteristic element size being equal to 0.00045 mm, two times smaller than the phase field length scale.
The finite element results obtained are shown in Fig. 10.Consider first Fig.10a, showing the distribution of moisture within the domain.In agreement with expectations, the case of moisture absorbed shows a higher content of moisture close to the surface exposed to the environment, while in the moisture dried case there is a reduction in the magnitude of C near the left free surfaces of the plate.Since the crack is introduced geometrically, it has a distinct impact on the moisture contours.Fig. 10b shows the reaction force evolution for the moisture absorbed and moisture dried cases.The latter exhibits a drop in the load carrying capacity as a result of the loss of moisture while a monotonic increase is observed for the moisture absorbed case, as represented in Fig. 10b, which shows the force versus time responses before mechanical loading is considered.The mechanical and fracture responses obtained after applying a mechanical load are shown in the bottom row of Fig. 10.It can be observed that, for both the damage contours (Fig. 10c) and the force versus displacement curves (Fig. 10d), material behaviour is dominated by mechanical effects and not moisture.This is arguably related to the small degree of damage that is triggered by the moisture, for the material parameters adopted and the environmental conditions considered -see Figs.5c and 7b, the magnitude of φ remains small throughout the analysis for similar environmental conditions.In all three cases (no moisture, moisture absorbed, moisture dried), matrix cracking and multiple matrix-fibre interface debonding events are observed and these lead to a reduction in load carrying capacity.Fig. 7d shows multiple load drop events which are associated with the coalescence of interface debonding and matrix cracks.The peak load of the moisture absorbed model is slightly lower than the one attained for no moisture and moisture dried models but differences are nevertheless small.Considering a moisture concentrationdependent toughness could bring larger differences but this dependency is not clear from the experimental side [57].

Meso-scale: ply-level model
A ply-level model is investigated using the same geometry, and boundary conditions as Chilali et al. [2] (see Fig. 11a).The ply dimensions are L = 10 mm and H = 1.5 mm, with the fibre diameter being d = H × V f = 0.24 mm, where the volume fraction equals V f = 32%.The symmetry conditions are located at the bottom and right boundary.A moisture content of 7.45% is prescribed at the top and left boundaries of an initially dry sample for a total of 25×10 6 seconds, which is sufficient to reach the equilibrium state -see Fig. 11b.After that, the second step reduces the moisture content to 0% at the same boundary for 50×10 6 seconds, as shown in Fig. 11c.A characteristic element size of 0.013 mm was used, two times smaller than the phase field length scale, which results in a FE model of 107,333 8-node quadratic elements.The result obtained are shown in Fig. 12.The moisture concentration distribution at time t = 1 × 10 6 s (Fig. 12a) shows a sharper gradient along the vertical direction, relative to the horizontal one.More importantly, as shown in Fig. 12b, the damage resulting from hygroscopic expansion appears to concentrate mainly on the fibre while also spreading slightly into the interface.The internal horizontal stresses resulting from hygroscopic swelling achieve their maximum values at the interface, while notable compressive horizontal stresses are attained inside of the fibre (see Fig. 12c).In terms of the σ yy stress component, its maximum values are attained at the twist of the fibre, see Fig. 12d.In addition, this case study enables us to compare our numerical predictions with experiments.Specifically, as shown in Fig. 12e, moisture absorption predictions can be compared with the experimental results by Chilali et al. [2].In turn, the same figure illustrates the sensitivity of moisture transport to the flax fibre diffusion coefficients.The results reveal that, while the value of D = 1.19 × 10 −6 mm 2 /s adopted (see Table 1) delivers a result that is close to the experimental one, a perfect agreement requires taking a diffusion coefficient of D = 3.47 × 10 −4 mm 2 /s, which is close to the value of 2×10 −4 mm 2 /s experimentally measured by Célino et al. [58] at a relative humidity of 80%.Finally, Fig. 12f shows the horizontal expansion of the laminate due to hygroscopic expansion.During the absorption of moisture, the laminate undergoes a linear horizontal expansion depicted by the movement of the free vertical edge, which ends when the equilibrium state of moisture is reached.At this point, the expansion stops leading to a maximum increase of the laminate length of 0.63 mm.The composite starts to contract during the drying process and the original shape is effectively recovered after 5.0×10 7 s.

Macro-scale: laminate-level model
The last case study involves a macro-scale, laminate-level model inspired by the work of Chilali et al. [2].Thus, a laminate model comprising 8 plies with a layup of [0 • /90 • ]2s and dimensions L = 10 mm and H = 6 mm is adopted -see Fig. 13a.The fibre diameter is assumed to be d = 0.24 mm.Due to symmetry conditions at the bottom and right boundaries, only a quarter of the model is simulated.A moisture concentration of 7.45% is enforced at the top and left boundaries for 1×10 8 s which, as shown in Fig. 13b, is sufficient to achieve the equilibrium state.This is followed by a drying stage where the moisture concentration is fixed at C = 0 in the same boundaries (top and left) for further 3×10 8 s.The entire domain is discretised using 113,565 8-node quadratic elements, with the characteristic element size being equal to 0.013 mm, two times smaller than the phase field length scale.The computational results obtained are shown in Fig. 14 in terms of moisture content contours (Fig. 14a), elongation as a function of time (Fig. 14b), phase field damage contours (Fig. 14c), and σ xx (Fig. 14d) σ yy (Fig. 14e) stress contours.The moisture contours show gradients in the horizontal and vertical directions, as it could be expected from the boundary conditions adopted.Fig. 14b reveals that the composite elongates by up to 0.61 mm before returning to its original shape as a result of the drying process.Damage appears to concentrate in the longitudinal fibres but a degree of damage is also observed in the transverse ones (see Fig. 14c).The location of the maximum σ xx tensile stresses is at the interface between the longitudinal fibres and the matrix (Fig. 14d), while vertical tensile stresses appear to localise at the interface between the matrix and the transverse fibres (see Fig. 14e).The distinctive behaviour observed for the longitudinal and transverse plies is the result of the assumed anisotropy in Young's modulus, Poisson's ratio and hygroscopic expansion coefficient.Overall, the results suggest that moisture will induce a higher degree of damage in the longitudinal plies of composite laminates.

Conclusions
We have presented a new phase field-based multi-physics framework to model moisture-induced degradation in composite materials.The framework couples the diffusion of moisture, the associated hygroscopic expansion and the mechanical and fracture behaviours of the matrix, fibres and matrix-fibre interfaces.Also, a novel diffuse interface approach is presented to interpolate relevant properties along the fibre-matrix interface.The framework is numerically implemented using the finite element method and numerical experiments are conducted to gain insight into moisture-material interactions at the micro-, meso-and macro-scales.The main findings are: • The micro-scale analysis of a single-fibre model shows how gradients in concentration lead to hygroscopic strains that result in damage accumulation at the fibre-matrix interface (debonding).• Uniformly distributed and randomised multi-fibre models show that the fibre distribution plays a secondary role for the material parameters and conditions assumed (equal volume fraction and fibre radius).• Micromechanical models involving multi-fibres and combined mechanical and hygroscopic straining show multiple matrix-fibre interface debonding events but a small influence of moisture on the failure process (for the properties and environments considered).• A meso-scale ply-level model reveals a good agreement with experimental measurements of moisture uptake and predicts that damage localisation will take place mainly along the flax fibre.• The macro-scale laminate-level model shows that moisture induces a higher degree of damage on the longitudinal ply, relative to the transverse due to the anisotropy of flax fibres.
The developed framework constitutes a comprehensive multiscale virtual tool to understand the environmentally-assisted degradation of composite materials.There are multiple avenues of future work.One is to further enhance the efficiency of the calculations by incorporating adaptive mesh refinement techniques, which have shown to be compelling in handling large-scale problems [59][60][61][62].From a physical viewpoint, other possibilities include accounting for the role of hydrostatic stresses on mass transport and the sensitivity of the mechanical properties to moisture content.model using initial conditions and predefined field variable.
The first step involves introducing a phase field distribution so as to interpolate the bulk and interface material properties accordingly.This is an initialisation step that is carried out in the absence of any mechanical or chemical load.To achieve this conveniently in ABAQUS, the analogy between the partial differential equations describing phase field evolution and heat transfer is exploited by means of a HETVAL subroutine [63,64].Thus, the equation describing the evolution of the indicator parameter d, Eq. ( 17), can be rearranged as where d is the length scale parameter.This equation has the same structure as the heat transfer equation under steady-state conditions when a source term r due to internal heat exists; Here, T is the temperature field and k is the thermal conductivity.Thus, the T ≡ d analogy can be readily attained by defining k = 1 and formulating the source term as, In addition, the HETVAL subroutine requires the definition of the rate of change of r relative to the primary field, which is given by ∂r ∂T ≡ ∂r ∂d = − 1 This HETVAL subroutine is used to run a pre-processing steady-state heat transfer analysis in which the Dirichlet boundary condition T ≡ d = 1 is enforced in the set of nodes that describes the (sharp) interface.Then, the physical simulation is conducted, recovering the diffusion interface through the use of a so-called predefined field variable.For this purpose, an initial step with initial conditions is defined, which reads the temperature distribution (output variable NT) from the pre-processing analysis.This field variable is used to interpolate the material properties within a user subroutine.The approach is general, in that it can be combined with any choice of user subroutine (UMAT, UELMAT, UEL) in the second step.A simple example of this implementation, accompanied by documentation, is made available to download at www. imperial.ac.uk/mechanics-materials/codes and www.empaneda.com/codes.

Declarations
• Conflict of interest: The authors declare that they have no known competing interests or personal relationships that could have appeared to influence the work reported in this paper.
• Data availability: All data and codes available from the author upon reasonable request.

Fig. 3 :
Fig. 3: Illustration of the diffused interface approach in a single fibre model: (a) phase field interface indicator d, and (b) interpolation of the material toughness across the matrix-interface-fibre domains.Subfigure (a) shows how the phase field indicator d varies between the interface (d = 1) and the bulk (d = 0) regions, while subfigure (b) illustrates how a material parameter is varied across the fibre (high toughness), interface (low toughness) and matrix regions (intermediate toughness).

Fig. 4 :
Fig. 4: Interfacial damage in a single-fibre model due to hygroscopic swelling: (a) geometry and boundary conditions, and (b) evolution of the moisture content in the centre of the domain as a function of time, showing a moisture absorption stage, followed by an equilibrium plateau and a final drying phase.

Fig. 5 :
Fig. 5: Interfacial damage in a single-fibre model due to hygroscopic swelling: contours of (a) moisture concentration, (b) displacement magnitude, (c) phase field damage, (d) σ xx stress component, and (e) σ yy stress component.The results are shown for times of t = 100 s (a) and t = 2, 000 s (b)-(e).Finally, subfigure (f) shows the evolution of the vertical reaction force in time.

Fig. 6 :
Fig. 6: Interfacial damage in a multi-fibre model due to hygroscopic swelling: (a) geometry and boundary conditions for both the Square Arrayed (SA) and Randomly Distributed (RD) cases, and (b) evolution of the moisture content in the centre of the domain as a function of time, showing a moisture absorption stage, followed by an equilibrium plateau and a final drying phase.

Fig. 7 :
Fig. 7: Interfacial damage in a multi-fibre model due to hygroscopic swelling: contours of (a) moisture content after 2,000 s, (b) phase field damage after 30,000 s, and (c) σ xx stress component after 30,000 s.The figures in the top correspond to the Square Arrayed (SA) multi-fibre model while those in the bottom have been obtained with the Randomly Distributed (RD) multi-fibre model.

Fig. 8 :
Fig.8: Interfacial damage in a multi-fibre model due to hygroscopic swelling: evolution of the vertical reaction force in time for both the Square Arrayed (SA) and Randomly Distributed (RD) multi-fibre models.

Fig. 9 :
Fig. 9: Micro-scale chemo-mechanical analysis of a single-edge cracked plate: (a) geometry and boundary conditions, and (b) evolution of the moisture content in the centre of the domain as a function of time for the cases of moisture absorbed (MA) and moisture dried (MD).

Fig. 10 :
Fig. 10: Micro-scale chemo-mechanical analysis of a single-edge cracked plate: (a) moisture concentration contours at relevant times, (b) vertical reaction force versus time before the mechanical load is applied, (c) phase field damage contours showing matrix cracking and interface debonding, and (d) force versus displacement responses.Results are obtained for three case studies: no moisture (NM), moisture absorbed (MA) and moisture dried (MD).

Fig. 11 :
Fig. 11: Ply-level model subjected to moisture absorption: (a) geometry and the boundary conditions of the composite, and (b) evolution of the moisture content in the centre of the domain as a function of time, showing a moisture absorption stage, followed by an equilibrium plateau and a final drying phase.

Fig. 13 :Fig. 12 :
Fig. 13: Laminate macro-scale model subjected to moisture absorption: (a) geometry and boundary conditions of the composite, (b) evolution of the moisture content in the centre of the domain as a function of time, showing a moisture absorption stage, followed by an equilibrium plateau and a final drying phase.

Table 1 :
Material properties of the flax fibre, epoxy matrix and epoxy-flax interface.