Multi-scale computational homogenisation to predict the long-term durability of composite structures

A coupled hygro-thermo-mechanical computational model is proposed for fibre reinforced polymers, formulated within the framework of Computational Homogenisation (CH). At each macrostructure Gauss point, constitutive matrices for thermal, moisture transport and mechanical responses are calculated from CH of the underlying representative volume element (RVE). A degradation model, developed from experimental data relating evolution of mechanical properties over time for a given exposure temperature and moisture concentration is also developed and incorporated in the proposed computational model. A unified approach is used to impose the RVE boundary conditions, which allows convenient switching between linear Dirichlet, uniform Neumann and periodic boundary conditions. A plain weave textile composite RVE consisting of yarns embedded in a matrix is considered in this case. Matrix and yarns are considered as isotropic and transversely isotropic materials respectively. Furthermore, the computational framework utilises hierarchic basis functions and designed to take advantage of distributed memory high-performance computing.


Introduction
Fibre reinforced polymer (FRP) composites have exceptional mechanical and chemical properties, including light weight, high specific strength, fatigue and corrosion resistance, low thermal expansion and high dimensional stability. They are commonly used in engineering application including aerospace, ships, offshore platforms, automotive industry, prosthetics and civil structures [1,2]. Textile or woven composites is a class of FRP composites, in which interlaced fibres are used as reinforcement, which provides full flexibility of design and functionality due to the mature textile manufacturing industry [3]. A detailed review, explaining the design and fabrication of textile preforms including weaving, knitting, stitching and braiding with their potential advantages and limitations is given in [4]. As compared to the standard laminated composites, textile composites have better damage and impact resistance, better through-thickness properties and reduced manufacturing cost. However, waviness of the yarns in the textile composites reduces the tensile and compressive strengths [5].
Due to their complicated and heterogeneous microstructure, Computational Homogenisation (CH) provides an accurate modelling framework to simulate the behaviour of FRP composites and determine the macro-scale homogenised (or effective) properties, including mechanical stiffness, thermal conductivity and moisture diffusivity, based on the physics of an underlying, microscopically heterogeneous, representative volume element (RVE) [6][7][8][9][10].
The homogenised properties calculated from the multi-scale CH are subsequently used in the numerical analysis of the macro-level structure. A variety of analytical and numerical homogenisation schemes have been developed for textile based FRP composites, which are normally based on the existence of an RVE and focus attention on the mechanical behaviour. Analytical methods are quick and easy to use but generally give poor estimates of the homogenised properties and are normally based on oversimplified assumptions of the microstructure and states of stress and strain. In the literature, some of the analytical homogenisation schemes, with their potential applications and limitations highlighted, are given in [10][11][12][13][14]. Numerical techniques, on the other hand, can accurately estimate the homogenised properties by capturing accurately the intricate micro-structure exactly but are computationally expensive. Examples of numerical homogenisation schemes applied to FRP composites can be found in [15][16][17][18][19][20]. A review article, summarising some of the analytical and numerical homogenisation techniques for the mechanical properties of the textile composite is given in [21].
During their service lives, FRP structures can be exposed to harsh hygrothermal environmental conditions in addition to mechanical loading, which can lead to matrix plasticisation, hydrolysis and degradation of fibres/matrix interfaces [22][23][24]. In the long-term, these processes significantly reduce the mechanical performance of these structures. Therefore, understanding heat and moisture transport mechanisms and their effect on the mechanical performance are fundamental for assessing the long-term use of FRP structures. Effective diffusivities of FRP composites with impermeable fibres were studied in [25] within the context of FEM, considering the variation in fibre volume within square and hexagonal unit cells. In [26], effective moisture diffusivity as a function of temperature and fibre volume fraction were investigated for FRP composites with permeable fibres using a unit cell approach. For textile composites, moisture transport was studied in [22,23] as a function of variation in tow architectural parameters, e.g. tow waviness, tow crosssection shape and wave pattern. Transient moisture transport in multilayer textile composites was investigated in [27]. An analogy between thermal and moisture transport analysis was used in [28] to study the moisture diffusion and corresponding weight gain for carbon braided composites. In [24,29], a multi-scale CH framework based on the hygro-mechanical analysis was proposed while using a two-dimensional RVE with randomly distributed fibres in 0 o and 90 o directions. A two way coupling was considered between the mechanical and moisture transport analysis and a model reduction scheme was used to reduce the computational cost. For the composite material, deformation dependent diffusion at finite strains was considered in [30]. Masonry wall reinforced with FRP reinforcement was studied in [31] while considering hygro-thermo-mechanical analysis. A recent review article [32], explains different degradation mechanism for FRP composites in connection with different environmental conditions.
In this paper, a coupled hygro-thermo-mechanical computational framework based on the multiscale CH is proposed for FRP composites. At each integration point, an RVE consisting of single plain weave textile composite is considered, consisting of yarns embedded in the matrix. Elliptical cross sections and cubic spline paths are used to model the geometry of these yarns.
Separate RVEs are considered for the heat transfer, moisture transport and mechanical CH. One-way coupling is considered in this case, i.e. mechanical analysis is assumed to be dependent on both moisture transport and thermal analyses but any influence on the moisture or thermal behaviour due to the mechanical behaviour is ignored. A degradation model, developed from experimental data relating evolution of mechanical stiffness over time for given exposure temperatures and moisture concentration was also developed and incorporated in the proposed computational framework. A unified approach is used to impose the RVE boundary conditions, which allows convenient switching between the different RVEs boundary conditions (linear Dirichlet, uniform Neumann and periodic) [33]. For a given size of RVE the periodic boundary conditions gives a better estimation of the homogenised material properties as compared to linear Dirichlet and uniform Neumann boundary conditions, which give an upper and lower limit [33][34][35] respectively.
The developed computational framework utilises the flexibility of hierarchic basis functions [36], which permits the use of arbitrary order of approximation, thereby improving accuracy for relatively coarse meshes. For the thermal and moisture transport analyses both matrix and yarns are assumed as isotropic materials, while for the mechanical analysis, the yarns are considered as transversely isotropic materials. The required principal directions of the yarns for the transversely isotropic material model are calculated from potential flow analysis along these yarns. Furthermore, distributed memory high performance computing is used to reduce the computational cost associated with the current multi-scale and multi-physics computational framework.
This paper is organised as follows. The multi-scale CH framework and corresponding implementation of the RVE boundary conditions are described in §2. Transient heat and moisture transport analyses along with their FE implementation are discussed in §3. The derivation of the degradation model from the experimental data is next explained in §4. The overall multi-scale and multi-physics computational framework is described in detail in §5. Computation of yarns directions are explained in §6. A three-dimensional numerical example and concluding remarks are given in §7 and §8 respectively.

Multi-scale computational homogenisation
In multi-scale CH, a heterogeneous RVE is associated with each Gauss point of the macro-homogeneous structure. Multi-scale CH gives us directly the macro-level constitutive relation, allows us to incorporate large deformation and rotation on both micro-and macro-level and both physical and geometric evolution can be included on both micro-and macro-level [34]. The multi-scale CH procedure and corresponding implementation of RVE boundary conditions is described initially for the mechanical case, which is subsequently extended to corresponding thermal and moisture transport cases.
The first order multi-scale CH is used in the paper, the basic principle of which is shown in Figure 1, where Ω ⊂ R 3 and Ω µ ⊂ R 3 are macro and micro domains respectively. Macro-strain ε = ε 11 ε 22 ε 33 2ε 12 2ε 23 2ε 31  the micro length scale is considered to be very small compared to the macro length scale and the macro-strain field attributed to each RVE is assumed to be uniform. Therefore, first-order CH is not suitable for problems with large strain gradient (but can be used for problems subjected to large strain) and cannot be used to take into account micro-level geometric size effects [7, 34].
On the micro-level at any point y = y 1 y 2 y 3 T the displacement field is written as [37][38][39] where εy is a linear displacement field and u µ is a displacement fluctuation.
The micro-strain associated with point y is written as whereε (y) = ∇ sũ µ is the strain fluctuation at the micro-level and ∇ s is the symmetric gradient operator. Furthermore, volume average of the microstrain is equivalent to the macro-strain: where V is the volume of the RVE. It is clear from Equation (3) that the volume average of the strain fluctuation is zero, i.e.
The micro-equilibrium state in the absence of body force is written as or in the variational form where t is an external applied traction and η µ is a virtual displacement.
No constitutive assumption is made for the macro-level problem, and the macro-stress σ is determined by volume averaging of the micro-stress σ µ , i.e.
To determine consistent RVE boundary conditions Hill-Mandel principle is used, which equates the work done at the micro and macro levels, i.e.
After combining Equations (8), (7), (6) and (1), we obtain ∂Ωµ t (y) · η∂Ω µ = 0, i.e. the virtual work associated with the traction t vanishes and the variational form of the equilibrium Equation (6) reduce to Thus, on the RVE level, the problem reduces to the calculation of the displacement fluctuation u µ for a given macro-strain ε.
In this work, we consider three types of RVE boundary conditions, which can be shown to satisfy the Hill-Mandel principle [33,35] 1. Linear boundary displacement (Dirichlet): In this case, the displacement fluctuation u µ are assumed as zero over the boundary of the RVE, which leads to fully prescribed displacements on the RVE boundary 2. Periodic boundary conditions: In this case it is assumed that displacement fluctuation u µ is periodic and traction is anti-periodic, i.e.
where y + and y − are a pair of points on the opposite boundary of the RVE, i.e. ∂Ω + µ and ∂Ω − µ respectively. 3. Uniform traction (Neumann) boundary conditions: In this case the traction on the boundary of the RVE is written in term of the macro- where n is the outward normal of the RVE boundary.
The imposition of these different type of RVE boundary conditions leads to different responses of the RVE. The linear displacement boundary condition gives the stiffest response followed by periodic, while the uniform traction boundary condition imposes the least kinematic constraint. In this paper, RVE boundary conditions are prescribed using the procedure given in [33] with extension to three-dimensions. The final discretised system of equations in the case of the RVE is written as where K and u are the standard FE stiffness matrix and displacement vector and λ is the unknown vector of Lagrange multipliers required to impose the RVE boundary conditions. Matrices C and D are given as In Equation (15), N is a matrix of shape functions and X is a matrix of spatial coordinates, evaluated at Gauss points during numerical integration of the surface integrals in Equation (15) and is given as Furthermore, H is a matrix that is specific to the type of boundary conditions used, each row of which represents an admissible distribution of nodal traction forces on the RVE boundary [33]. The specific choice of H in case of linear displacement, periodic and uniform traction boundary conditions can be found in [33,39] and is not repeated here. Method of static condensation and projection matrices are alternative procedures to solve system of equations (14) [33].
To determine the homogenised stress, the right hand side of Equation (8) is written in term of surface quantities, i.e.
Furthermore, the work done by the traction on displacements is equal to the work of the generalised tractions (Lagrange multiplers) on the generalised displacements (strains): With reference to Equations (17) and (18), the homogenised stress can be determined as: For the macro-level finite element analysis, material stiffness matrix C need to be determined at each Gauss point. The relationship between macro-stress and strain is written as To compute the homogenised stiffness matrix C, six linear system of equations for the RVE need to be solved for a given set of unit strain vectors.
This will give a set of homogenised stresses, i.e.
where for example: In each of the six cases, only the right-hand side of the system of Equations (14) changes, which is solved very efficiently as the left-hand side matrix is factorised only once.
Computational homogenisation for thermal and moisture transport cases are very similar and is presented as a single case in the following. The com-plete derivation is not repeated due to its similarity with the mechanical counterpart but only the key differences are highlighted. The gradient of T is used to formulate the boundary value problem on the micro-level, where ψ = T, c is a scalar field (temperature and moisture concentration respectively). The homogenised flux and corresponding conductivity matrix K ψ are next obtained after solution of RVE's boundary value problem. The complete homogenisation process for this case is also shown in Figure 1. Similar to Equation (14), the final discretised system of equations in this case is written where K ψ and ψ µ are the standard FE conductivity matrix and vector of field values (temperature or moisture concentration) respectively. λψ is a vector of Lagrange multipliers, requires to impose the RVE boundary conditions.
In contrast to the three Lagrange multipliers per node for the mechanical case, only one Lagrange multiplier per node is required for temperature and moisture transport RVEs. Equations for both C ψ and D remains the same as their mechanical counterpart in Equation (15), albeit with reduced dimensions due to there being only one degree of freedom per node. The matrix for special coordinates, Equation (16) is modified as Furthermore, with reference to Equation (19), the homogenised flux is written as: The homogenised flux q ψ and gradient of field values ∇ψ are related with the following equation:

Macro-level transient field problems
FRP structures on the macro-level are generally exposed to hygro-thermal environmental conditions and require a detailed heat and moisture transport analysis. In this paper, conduction and diffusion models are considered for the heat and moisture transport analysis respectively, both of which can be represented by the following governing equation (conservation of energy for heat transfer and conservation of mass for moisture transport) where t is time, ρ and c p are macro-level density and specific heat capacity ψ are thermal or moisture conductivities in x 1 , x 2 and x 3 directions respectively, and are assumed to be independent of temperature T or moisture concentration c. In Equation (27), it is assumed that heat and moisture transport are governed by Fourier's law and Fick's law respectively. Fourier's law for the heat conduction is expressed as where is a vector of heat fluxes and K T is matrix of thermal conductivities and is written as Similarly, Fick's law is written as where is a vector of moisture fluxes and D c is matrix of moisture diffusivities and is written as where D c are moisture diffusivities in x 1 , x 2 and x 3 directions. Using Equation (27), the relationship between diffusivity and conductivity is written as Equation (27) is discretised in the standard way: with N = N 1 N 2 · · · N n and ψ = ψ 1 ψ 2 · · · ψ n T the vectors of shape functions and nodal degrees of freedom. Thus, the discretised form Equation (27) is written as [40] C ∂ψ ∂t where and where q s = q ψ n with n = n x 1 n x 2 n x 3 T , the vector of normal direction cosines. Finally, B is the matrix of shape functions derivatives and ∂Ω q is the boundary with specified fluxes.

Degradation model
is assumed, where T is exposure temperature in o C, G | T is shear modulus at temperature T , G o = 3.76 GPa is shear modulus for the dry FRP sample, where α (T ) is function of temperature and is written as where T g is the glass transition temperature and β is a model parameter. In leading to a value of β = −0.001682. Using this value, a comparison with the fitted exponential curves is shown in Figure 3 indicating very good agreement.
Furthermore, the least square fitted curve adheres to the aforementioned assumption for the rate of degradation at T = 0K and T = T g . Combining Equation (39) and (40), the degradation model is written as A comparison between the experimental data, fitted exponential curves and the proposed degradation model (Equation (42)) is shown in Figure 5.
The degradation model in Equation (39) In this paper, due to the lack of experimental data, a simple linear relationship for γ (c) is proposed (shown in Figure 6), i.e. for fully saturated sample We have assumed that the degradation process is very slow compared to the daily variation of temperature T and moisture concentration c. This  the assumption that it degraded in a similar way as shear modulus, and is expressed as where ω is the isotropic damage parameter, with ω = 0 equivalent to no degradation and ω = 1 is full degradation. This degradation model is discre-

Computational framework
The proposed computational framework, composing a series of micro and macro-level analyses, is implemented in our group's FE software MoFEM (Mesh Oriented Finite Element Method) [41]. The detailed flow chart for the computational framework is shown in Figure 7. In this paper, one-way coupling is assumed, i.e. mechanical analysis is dependent on thermal and moisture transport analyses but not vice-versa. Furthermore, it is also assumed that thermal and moisture transport analyses are independent of each other. The conductivity matrices K T and K c for these analyses are deter- The residual associated with Equation (46) is written as Discretisation of the degradation parameter follows the standard form: where (1 − ω) h is the approximation of (1 − ω) and W (1−ω) is vector of unknown nodal values of (1 − ω). In the discretised form Equation (47) is written as Furthermore, the Jacobean associated with Equation (47) is written as where (

Heuristic method for computation of yarn directions
In order to re-orient the known stiffness matrix for the transversely isotropic material from the local coordinates to global coordinates, the yarn directions at each Gauss point need to be determined. To do this it is possible to simply use the cubic splines that were used to construct the yarns. However, this can lead to inaccuracies in the case of yarns with non-uniform cross-sections along their length. An alternative approach is proposed here in which the yarn directions are determined by solving the potential flow field separately along each yarn. The governing equation for potential flow is given as where φ is the stream function, the gradient of which, i.e. ∇φ = v subsequently determines the yarn directions. A detailed description of how to transfer the stiffness matrix for transversely isotropic material between local and global coordinate axes is given in previous publications [39,42,43].

Numerical example
A macro-level three-dimensional plate structure with a hole is considered, the geometry for which is shown in Figure 8 The RVE used in this case, comprising a plane wave textile composite, is shown in Figure 9(a). The geometry parameters for the RVE are shown in Table 1. Both macro-structure and RVE are discretised with tetrahedral elements with 780 elements in the case of macro-structure and 10,285 elements in the case of the RVE. An example of the yarn directions for the currently used RVE calculated from the potential flow analysis is shown in Figure 9(b). For the matrix material, thermal conductivity, density, specific heat and moisture diffusivity are 190 kg mm/s 3 C, 1.2x10 −6 kg/mm 3 , shown in Figure 13(c). It is clear that this increase in vertical displacement at constant load is due to the degradation of the macro-structure stiffness.

Concluding remarks
In this paper, an ongoing work, consisting of a coupled hygro-thermomechanical computational framework based on multi-scale CH is described    [5] A. Miravete (Ed.), 3-D textile reinforcements in composite materials,

Fibres properties Matrix Properties
(mm)