Three-dimensional nonlinear micro/meso-mechanical response of the fibre-reinforced polymer composites

A three-dimensional multi-scale computational homogenisation framework is developed for the prediction of nonlinear micro/meso-mechanical response of the fibre-reinforced polymer (FRP) composites. Two dominant damage mechanisms, i.e. matrix elasto-plastic response and fibre-matrix decohesion are considered and modelled using a non-associative pressure dependent paraboloidal yield criterion and cohesive interface elements respectively. A linear-elastic transversely isotropic material model is used to model yarns/fibres within the representative volume element (RVE). A unified approach is used to impose the RVE boundary conditions, which allows convenient switching between linear displacement, uniform traction and periodic boundary conditions. The computational model is implemented within the framework of the hierarchic finite element, which permits the use of arbitrary orders of approximation. Furthermore, the computational framework is designed to take advantage of distributed memory high-performance computing. The accuracy and performance of the computational framework are demonstrated with a variety of numerical examples, including unidirectional FRP composite, a composite comprising a multi-fibre and multi-layer RVE, with randomly generated fibres, and a single layered plain weave textile composite. Results are validated against the reference experimental/numerical results from the literature. The computational framework is also used to study the effect of matrix and fibre-matrix interfaces properties on the homogenised stress-strain responses.


Introduction
Compared to conventional materials, fibre-reinforced polymer (FRP) composites can offer exceptional physical and chemical properties (including high strength, low specific weight, fatigue and corrosion resistance, low thermal expansion and high dimension stability), making them ideal for a variety of engineering applications, including aerospace, marine, automotive industry, civil structures and prosthetics [1][2][3]. Phenomenological or macro-level models cannot accurately describe the complex behaviour of FRP composites due to their underlying complicated and heterogeneous microstructure. Furthermore, nonlinearities associated with the matrix elasto-plasticity and fibre-matrix decohesion make the computational modelling even more challenging. Multiscale computational homogenisation (CH) provides an accurate modelling framework to simulate the behaviour of FRP composites and determine the macro-scale homogenised (or effective) response, based on the physics of an underlying, microscopically heterogeneous, representative volume element (RVE) [4][5][6][7][8][9]3]. The homogenised properties calculated from the multi-scale CH are subsequently used in the numerical analysis of the macro-level structures.
A variety of numerical techniques have been developed to model the nonlinear micro-mechanical response of unidirectional (UD) FRP composites, mostly based on finite element analysis. For UD glass/carbon (G/C) FRP composites, a computational model was developed in [10,11] within the framework of finite deformation. Both in-plane shear and compressive loading scenarios were considered. The Mohr-Coulomb yield criterion and cohesive interface elements were used to model the response of epoxy matrix and fibre-matrix interfacial decohesion respectively. Fibres were generated randomly within the RVEs using the algorithm presented in [12] and were modelled as a linear-elastic and isotopic material. A parametric study, including the effect of matrix and interface properties on the stress-strain response, was also conducted. The idea of [10] was extended further in [13] by incorporating thermal residual stresses (due to cooling of FRP composites after curing process, caused by the mismatch in thermal expansion coefficients of matrix and fibres) in the simulation, in addition to transverse tensile and cyclic loading for the CFRP composites. The nearest neighbour algorithm (NNA) [14] developed by the same authors, was used to randomly generate the fibres within the RVEs. Using the same constitutive models for matrix, fibres and fibre-matrix decohesion as in [10], a multilayer multi-fibre (M 2 ) RVE was used in [15] for laminates. Each lamina was modelled as a cube with randomly distributed but axially aligned fibres, generated using a fibres randomisation algorithm in DIGIMAT FE [16]. Both cross [0/90] ns and angle [AE45] ns (where the subscript ns represents n layers with the same sequence and symmetric about the mid plane) GFRP composites were considered with in-plane shear loading and results of stress-strain behaviour were validated against the experimental results. A combined transverse compression and axial tension loading scenario was considered in [17] for UD GFRP composite. In addition to matrix plasticity and fibre-matrix decohesion, fibre breakage was also included in the FE simulation. The pressure dependent, Drucker-Prager yield criterion was used to model matrix plasticity and both fibre breakage and fibre-matrix interfacial decohesion were modelled with cohesive interface elements. A simple periodic, hexagonal fibre arrangement was assumed. In [18], a modified von Mises yield criterion was used to model the behaviour of the matrix material, while a maximum tensile stress criterion was used to model fibre breakage. Fibre-matrix decohesion was also included in the simulation and was modelled with cohesive interface elements. The random distribution of the fibres was also included within the RVE based on the optical microscopy of real composites. A variety of loading conditions was used subsequently to study the response of the UD FRP composite. The limitations of different plasticity models for modelling matrix materials, including Mohr-Coulomb and Drucker-Prager were argued in [19,20], especially in complex loading scenarios. Instead of the conventional plasticity models, a pressure dependent thermodynamically consistent plasticity model [21] was used. A statistically proven random distribution algorithm proposed by the same authors in [22] was used to randomly generate UD fibres within the RVEs. Similar to previous studies, fibres were modelled as linear-elastic and isotropic material and fibre-matrix decohesion was modelled with cohesive interface elements. A variety of RVE loading scenarios was considered including transverse tension and compression, transverse and longitudinal shear and combined transverse compression and transverse shear.
A number of numerical modelling approaches have been used to simulate the behaviour of textile composites subjected to different loading scenarios. A comprehensive review of these methods can be found in [23]. Continuum damage mechanics (CDM) was used in [24] to model both matrix and yarns for glass and carbon plain weave textile composites. Dissipated energy density was used as damage parameter and both material and geometric nonlinearities were included in the simulation. Further use of CDM in the simulation of textile composites can also be found in [25][26][27]. Moreover, a three-dimensional CDM based approach was used to simulate the progressive damage in laminated FRP composites in [28,29]. A variety of failure mode, including matrix tensile and compressive cracking, fibre tensile and compressive failure, fibre-matrix shearing and delimitation between the layers were included in the simulations. For a twill weave textile CFRP subjected to in-plane loading, a meso-mechanical analysis was performed in [30]. The matrix was modelled as elasto-plastic material with the same plasticity model as in [19][20][21], while yarns were modelled as linear-elastic and transversely isotropic material. Results of the RVE strain fields and homogenised stress-strain response were validated against the experimental results and found in a good agreement. These numerical simulations, described above, of FRP composite behaviour are limited to specific RVE type (2D or 3D, UD or woven/textile) or loading scenarios (normal or shear). In contrast, this paper develops a fully generalised three-dimensional micro/ meso-mechanical framework, which is subsequently implemented in the authors' open source FE software, MOFEM [31]. The dominant damage mechanisms (observed experimentally [10]), i.e. matrix elasto-plasticity and fibre-matrix decohesion, are included in the computational framework. Matrix material is modelled using a pressure dependent paraboloidal yield criterion [19,20,30,21] with an exponential hardening law. Fibre-matrix decohesion is modelled with zero thickness cohesive interface elements. Yarns are modelled as linear-elastic and transversely isotropic materials. Rather than simplified fibre arrangements for UD FRP composites, e.g. in [32,33,17,34], which are not the actual representation of the real FRP composites and can lead to erroneous results. this study adopts a statistically proven random distribution algorithm proposed in [22] to generate fibre arrangements within the RVE. The RVE boundary conditions are imposed in a unified manner which allows convenient switching between displacement, traction and periodic boundary conditions [35]. Hierarchic finite elements are adopted, which permits the use of arbitrary order of approximation, leading to accurate results for relatively coarse meshes. The computational framework is designed to take advantage of distributed memory high-performance computing. Moreover, CUBIT [36] and ParaView [37] are used as pre-and post-processor respectively. This paper is organised as follows. The computational framework is fully described in Section 2. The material models are given in Section 2.1, consisting of material model for matrix (Section 2.1.1), yarns/fibres (Section 2.1.2) and fibre-matrix interfaces (Section 2.1.3). The nonlinear multi-scale CH with corresponding RVE boundary conditions are explained in Section 2.2. Calibration and validation of the matrix plasticity model is given in Section 3.
Three numerical examples are given in Section 4, including UD GFRP composites (Section 4.1), M 2 RVE (Section 4.2) and plain weave textile composites (Section 4.3). Finally, the concluding remarks are given in Section 5.

Computational framework
The computational framework developed for FRP composites consists of a set of constitutive models for individual components including matrix, yarns/fibres and fibre-matrix interface and implemented within the formulation of first-order multi-scale CH.

Material constitutive models
Typical RVEs in the case of UD FRP and textile composites are shown in Fig. 1(a) and 1(b) respectively, consisting of yarns/fibres embedded within a polymer matrix. The constitutive model for FRP composites is a combination of constitutive models for these individual components, together with fibre-matrix interface decohesion. In the following, each of these constitutive model is explained in detail.

Matrix
The polymer matrix is modelled as an elasto-plastic material using a non-associative pressure dependent paraboloidal yield criterion [21,19,20,30,9]. This plasticity model can incorporate different yield strengths in tension and compression and is shown Fig. 2 (a) in the principal stress space. The yield function is written as where r is Cauchy stress tensor, I 1 ¼ trðrÞ is the first invariant of Cauchy stress tensor, J 2 ¼ 1 2 g : g is the second invariant of deviatoric stress g ¼ r À 1 3 I 1 and r t and r c are yield strengths in tension and compression respectively. A non-associative flow rule is used, for which the plastic potential function is written as where m plas is a material parameter and is known as plastic Poisson's ratio. Furthermore, the Helmholtz free energy in the case of linear isotropic hardening is written as where k and l are the Lame parameters, e is the strain tensor, r t 0 and r c 0 are the initial yield strengths in tension and compression respectively, a 0 and a 1 are internal kinematic variables and H t and H c are hardening parameters in case of tension and compression respectively. Following Eq. (3), yield strengths in tension and compression are written as As compared to the linear, a more realistic, exponential hardening law is presented in this paper, due to which Eq. (3) is rewritten as where H t ; H c are the difference between the ultimate and yield strengths, n t and n c are material parameters and determine the rate of convergence between yield and ultimate strengths. Following Eq. (5), yield strengths in tension and compression are written as A similar hardening law as a function of equivalent plastic strain was also used in [30]. In addition to linear and exponential, the developed computational framework can accommodate any isotropic and kinematic hardening laws, for further details, the readers are referred to [31].

Yarns/fibres
In textile composites, as shown in Fig. 1(b), yarns are interwoven together to form a textile structure. On the micro-level, yarns are the same as UD FRP composites consisting of bundles of glass/carbon fibres within the polymer matrix. On the meso-level, yarns are modelled as homogenous and transversely isotropic material with homogenised or effective properties obtained from the multi-scale CH of the UD FRP composites. Five material parameters are required for transversely isotropic materials, which are E p ; m p ; E z ; m pz and G zp , where z and p are fibres and transverse directions respectively as shown in Fig. 2(b). E p and m p are Young's modulus and Poisson's ratio in the transverse direction respectively, while E z ; m pz and G zp are the Young's modulus, Poisson's ratio and shear modulus in the fibre directions respectively. In order to reorient the known stiffness matrix for the transversely isotropic material from the local coordinates to global coordinates, the yarns 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 used in this paper, in which the yarns directions are determined by solving the potential flow along these yarns. A detailed description of this approach and how to transfer the stiffness matrix from the local to global coordinate axes is given in [6,3,[7][8][9]].
On the micro-level, in the case of UD FRP composites, fibres are modelled as linear-elastic and isotropic material (a special case of transversely isotropic material model), for which only two material parameters, i.e. Young's modulus and Poisson's ratio (E f , and m f ) are required. In this case, the following material parameters are used:  In the literature, glass fibres have been modelled as linear elastic and isotropic material (e.g. [20,15]) while carbon fibres are modelled both as isotropic elastic (e.g. [38]) and transversely isotropic elastic material (e.g. [39]). Therefore, the developed computational framework can be used for both glass and carbon fibres. Moreover, the framework allows for nonlinear behaviour for the fibres. However, at this stage, we restricted ourselves to linear.

Fibres/matrix interfaces
Fibre-matrix decohesion is modelled using standard cohesive interface elements with a straightforward material model, i.e. linear traction-separation law (shown in Fig. 2(c)). Only three material parameters are required for the material model, including cohesive strength f t , fracture energy G f (shaded area in Fig. 2)) and material parameter b, which assign different weight to opening and shear displacements. Mathematically the material model for cohesive interface elements is written as is the displacement jump with d n and d s as its normal and shear components and x is damage parameter. j is a history parameter and is equal to the highest value of displacement jump d. Furthermore, d 0 and d max are respectively the displacement jumps at the onset of damage (x ¼ 0) and when the interface is fully damaged (x ¼ 1). E 0 and d 0 are written as where E m is the Young's modulus of matrix material and h is the interface thickness. The initial value of E 0 and/or h must ensure displacement continuity in the case of no damage. These are related by Eq. (9). Theoretically, h tends to zero. However, to have a finite penalty term, h is chosen to be very small. Consistent with [10], From Fig. 2(c), the damage parameter x is written as The constitutive matrices for cohesive interface elements in the local and global coordinate systems are written as where I is the unit matrix and R is the transformation matrix [40]. Equations for stiffness matrix and corresponding internal forces are written as where U is a matrix of shape functions and t loc is the traction vector in local coordinates, details of which are given in [40].

Multi-scale computational homogenisation
In multi-scale CH, a heterogeneous RVE is associated with each Gauss point of the macro-homogeneous structure, the boundary conditions for which are implemented using the generalised procedure proposed in [3,35]. Small displacement and small strain formulations are used within the framework of first order multiscale CH, the basic concept of which is shown in Fig. 3 For a global step n þ 1, the discretised system of equations in case of an iteration i of the Newton-Raphson algorithm is written as where K and u are the standard FE tangent stiffness matrix and displacement vector respectively and k is the unknown vector of Lagrange multipliers required to impose the RVE boundary conditions. Matrix C in Eqs. (13) are calculated over the boundary C of the RVE and are constant throughout the calculations [3,35] and are given as In Eq. (14), N is a matrix of shape functions and 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 [35]. The specific choice of H in the case of linear displacement, periodic and uniform traction boundary conditions can be found in [35,3] and is not repeated here.
Matrix K comprises contributions of the matrix, yarns and yarnmatrix interface elements. F i nþ1 is a vector of residuals and is written as where X is a matrix of spatial coordinates, evaluated at Gauss points during numerical integration of the surface integrals and is given as

5: ð16Þ
At Newton-Raphson iteration i, variable n ¼ u; k is calculated using n i nþ1 ¼ n n þ P i m¼1 n m nþ1 . In Eq. (15), F int i nþ1 is a vector of internal forces. Furthermore, Cu i nþ1 and C T k i nþ1 are associated with the RVE boundary conditions and are written as where u h and k h are displacements and Lagrange multipliers calculated at a Gauss point, i.e. n h ¼ u h ; k h ¼ Nn e i nþ1 , where n e is a matrix of displacements or Lagrange multipliers associated with element e. Finally, the homogenised stress for global increment n þ 1 is written as: To compute the homogenised stiffness matrix C at the end of global increment n þ 1, the converged matrix K is subjected to six different macro-strain perturbations of unit vector leading to six linear system of equations. 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 Eqs. (13) changes, which is solved very efficiently as the left-hand side matrix is factorised only once.

Calibration and validation of plasticity model
Following [30], the plasticity model is first calibrated against the experimental results from [41,19,30] for epoxy resin subjected to tensile and compressive loading. A list of material parameters used in this case is shown in Table 1, where E; m; r to ; r to and m plas are given in [30,19]. Moreover, the hardening parameters, i.e. H t ; H c ; n t and n c are determined from the numerical simulation based on the experimental stress-strain curves. The estimated parameters H t ; H c are the same as given in [30] but in contrast, due to the use of hardening law as a function of internal kinematic variables leads to different n t and n c in our case.
The geometry considered in this case is a cube of dimension 1 mm, which is discretised with 1191 tetrahedral elements and 299 nodes. The cube is fixed at the bottom face and subjected to tension, compression and shear loading on the top face as shown in Figs. 4(a), 4(b) and 4(c) respectively. A comparison between numerical and experimental stress-strain responses for all the three loading scenarios are shown in Fig. 4(d). As expected the numerical and experimental responses in tension and compression are in good agreement, as a result of parameter fitting. The response in shear is not fitted and also shows fairly good   agreement. All of the three responses are highly non-linear and it is clear that the plasticity model can capture them well. It must be noted that the plasticity model requires input data for tension and compression and can be used subsequently to simulate more generalised loading scenarios. The aim here is to reproduce cases with homogeneous stress state ( [41,19,30]) with only matrix material. Moreover, softening in not included in the matrix model. Therefore, no shear band is observed, while modelling these cases. Furthermore, the plasticity model with exponential hardening law is unable to reproduce the hardening behaviour observed in [41] in the case of compression loading for the applied strain greater than 20%

Numerical examples
Three numerical examples are now given to demonstrate the correct implementation and performance of the developed computational framework.

Unidirectional GFRP composites
The first numerical example consists of polymer composites reinforced with unidirectional glass fibres subjected to transverse tension. A similar numerical example is also considered in [20]. Two RVE sizes, consisting of transverse side-measure of 10Â and 20Â the fibre radius were used in [20] for the evaluation of independence of the results from the RVE size. For the smaller RVE, five different fibre distributions were used. A similar study was also performed in [42]. It was concluded from the homogenised stress-strain responses that the pre-peak response and ultimate load is independent of the size of the RVE. In this paper, four different RVE sizes are considered, consisting of periodic, randomly distributed but axially aligned fibres and are shown in Fig. 5. The algorithm proposed in [22] is used to randomly generate the fibres within the RVEs with diameter of 5 lm and volume fraction of 60%.
The four RVEs are discretised with 1499, 3239, 13,023 and 21,140 tetrahedral elements. The detailed geometry for RVE-4, showing  the individual matrix, fibres and cohesive interface elements, is presented in Fig. 6(a).
A list of material properties used for the elasto-plastic matrix materials are the same as given in Table 1, while for the linearelastic and isotropic glass fibres Young's modulus and Poisson's ratio used are 74 GPa and 0.2 respectively. For the cohesive interface elements, interface strength and fracture energy are 50 MPa and 2 J/m 2 respectively. The material properties used here are consistent with [20]. The macro-strain (applied to the RVEs using periodic boundary conditions) versus homogenised stress response for all four RVEs are compared to reference numerical results from [20] and are shown in Fig. 6(b). The RVEs are subjected to a transverse strain e xx of 1 percent. It is clear from Fig. 6(b) that the developed computational framework accurately predicts the stress-strain behaviour in the pre-peak region (up to e xx ¼ 0:65%) for all of RVEs. The size effect can be clearly seen in the post-peak region (beyond e xx ¼ 0:65%); increasing the size of the RVE leads to a more brittle response. A similar behaviour was also reported in [20,43]. Issues related with the existence and size of the RVE and pre-and post-peak region behaviour are described in detail in [43][44][45][46], where the ill-posedness of the macro-level BVP and its non-objectivity with respect to the size of the RVE is discussed. A detailed description of dealing with these limitation of the classical CH schemes is given in [46]. These specialised treatments are not considered in this paper. The final damaged RVEs with clear localisation zones/debonding are also shown in Fig. 6 (b). The damaged zones consist of fully damaged cohesive interface elements that are perpendicular to the direction of the applied strain. It is clear that fibre-matrix decohesion interface leads to a reduction in load transfer from the matrix to the fibres, which results in the overall stiffness. Furthermore, strain localisation associated with the damaged zones subsequently leads to severe plastic deformation of matrix material.
A parametric study is also conducted to investigate the effect of different parameters on the macro-strain versus homogenised stress response. The effect of fracture energy on the stress-strain response is shown in Fig. 7(a), where fracture energies of 2, 3, 4 and 100 J/m 2 are used but all other parameters are kept constant. It is clear from Fig. 7(a) that the stress-strain response are the same for all the four cases up to e xx ¼ 0:6%. Furthermore, lower fracture energies leads to clear damaged zone with high strain localisation (severe plastic deformation). The effect of interface strength on the stress-strain response is shown in Fig. 7 parameters are kept constant. Furthermore, the effect of considering unlimited interface strength, i.e. f t ¼ 1 in shown in Fig. 7(b). A clear localisation can be seen for both f t ¼ 35 and f t ¼ 50 cases as compared to the f t ¼ 20 for the applied strain of 1%.
The effect of a linear elastic material, as opposed to an elastoplastic material, on the strain-stress response is shown in Fig. 7 (c). In addition, three cases with different interface strengths of cohesive interface elements, i.e. 20, 35 and 50 MPa are considered. It is clear from Fig. 7(c) that in the pre-peak regions, the use of either linear-elastic or elasto-plastic matrix material leads to almost similar stress-strain response while in the post-peak region the use of linear-elastic matrix material leads to relatively stiff response. The final damaged RVEs for both f t ¼ 20 and f t ¼ 35 with both linear-elastic and elasto-plastic matrix materials are also shown in Fig. 7(c). The high strain localisation in the damaged zones leads to severe plastic deformations leading to a more brittle stress-strain response.

Multi-fibres multi-layer RVE
A multi-fibre multi-layer RVE subjected to in-plane shear is considered in the second example. A similar example is also analysed experimentally and numerically in [15], the stress-strain response from which is used here as a reference. The UD FRP composite used in this case, consist of E-glass (ER-459L) and epoxy matrix (EPOFINE-556) with FINEHARD-951 hardeners. Parametric study was performed in [15] to determine the size of each cube within the RVE. Different RVE sizes, ranges from 0.1 mm (with 6 fibres) to 0.5 mm (with 155 fibres) were considered and the homogenised stress-strain responses were compared. It was concluded that the effect of the size of the RVE was not appreciable on the homogenised stress-strain response and RVE with 0.1 mm were used for all the subsequent analyses due to its lower computational cost. The RVE geometry is shown in Fig. 8(a), consisting of two cubes of dimension 1 mm with randomly distributed fibres (generated using the algorithm in [22]) of 24 lm and volume fraction of 28% and are placed on the top of each other with an angle of 90°. In Fig. 8(b) and 8(c) individual matrix and fibres are shown respectively. The RVE is discretised with 32,818 tetrahedral elements and is shown in Fig. 8(d), while fibre-matrix interfaces are discretised with 3056 cohesive interface elements and are shown in Fig. 8(e). Moreover, a perfect bond is assumed between laminae.
For the linear-elastic and isotropic glass fibres, Young's modulus and Poisson's ratio are 73 GPa and 0.23 respectively. For the matrix, most of the material parameters used are the same as given in Table 1 with the only exception of Young's modulus and Poisson's ratio, which are 4.7 GPa and 0.3 respectively. For cohesive interface elements, interface strength and fracture energy used are 30 MPa and 100 J/m 2 respectively. These material properties are consistent with [15]. The RVE in this example is subjected to shear strain c zx = 4%, as shown in Fig. 9. The shear stress versus shear strain (c zx versus s zx ) response is compared with the experimental and numerical results from [15] and is shown in Fig. 9, which are in a very good agreement. Stress-strain response is almost linear up to c zx = 1.5%, beyond which the response is nonlinear due to the decohesion between fibres and matrix. The difference between the numerical and experimental results (especially between c zx = 1% and 2.5%), might be due to the assumption of perfect bonding between the 0°and 90°laminae. At the end of the simulation, contours of c zx over the deformed RVE are also shown  in Fig. 9. Strain is very small in the glass fibres as compared to matrix material due to the associated high stiffness. Decohesion between matrix and fibres can also be seen Fig. 9.

Plain weave textile composites
Finally, a plain weave textile composite subjected to a variety of normal and shear loading conditions is considered, consisting of E-glass fibres and epoxy matrix. A similar numerical example is also considered in [24]. An RVE, consisting of similar yarns in warp and weft directions is used in this example, for which the geometry with all of the required dimensions are shown in Fig. 10(a). Elliptical cross-sections and cubic splines are used respectively to model the cross sections and paths of the yarns. The volume fraction of fibres within the yarns is 65% while the total volume fraction of fibres within the RVE is 35%. The RVE is discretised with 11,516 tetrahedral elements and is shown in Fig. 10(b). For the elastoplastic matrix material, the same properties are used as given in Table 1, while for the linear-elastic and transversely isotropic yarns, material properties are given in Table 2 [24]. Furthermore, a perfect bond is assumed between yarns and matrix.
The yarns direction, calculated from the potential flow analysis, are shown in Fig. 11(a). Four loading conditions, including two normal (e xx and e yy ) and two shear (c yz and c zx ) are considered, where the RVE is subjected to a macro-strain of 3%. For the shear case, the stress-strain responses are shown in Fig. 11(a). The nonlinear response beyond a strain of 1.5% is due to matrix failure. The stress-strain response in the case of c yz is also compared with numerical results from [24], which are in a very good agreement, especially in the linear region (up to c yz =1.5%). Beyond c yz =1.5%, our simulation result is relatively stiffer, which might be due to the use of linear-elastic material for the yarns. A high strain gradient can also be seen in the matrix, especially in regions of the thin matrix layer. The response in the case of c zx involves shearing of yarns leading to stiffer behaviour as compared to c yz . Furthermore, response in the case of e xx and e yy are shown in Fig. 11(b). For the given range of applied strains, stress-strain responses for both e xx and e yy are linear. Moreover, e xx involves direct tensile load on the yarns and behave stiffer as compared to e yy case, where strain is applied directly on the matrix.

Concluding remarks
A three-dimensional, nonlinear micro/meso-mechanical multiscale CH framework is developed for FRP composites. The matrix material is modelled as elasto-plastic, using a paraboloidal yield surface. Decohesion of the fibre-matrix interface are modelled using cohesive interface elements. The yarns/fibres are modelled as linear-elastic and transversely isotropic material. It is shown that the two dominant damage mechanisms, i.e. matrix plasticity   Fig. 11. Stress-strain responses for the plain weave textile composites subjected to different loading conditions. and fibre-matrix interfacial decohesion control the strength of FRP composites. Experimental stress-strain results for epoxy resin for both tension and compression load cases are used to calibrate the plasticity model and is validated subsequently for the shear loading. Three numerical examples with a variety of RVEs and loading conditions are considered to demonstrate the performance of the proposed computational framework. In both UD FRP composite and M 2 RVE examples, fibres are randomly generated within the RVEs using a statistically proven random distribution algorithm. In the UD FRP numerical example, the developed computational framework can accurately predict the stress-strain behaviour in the pre-peak region, while in the post-peak region size dependent response is observed, which is natural in the case of first-order computational homogenisation. A parametric study is also conducted for the UD FRP numerical example, i.e. the effect of different matrix and fibre-matrix interface parameters on the stress-strain behaviour and it is shown that failure starts at fibre-matrix interface followed by localised deformation and matrix plasticity. Furthermore, from the M 2 RVE and plain weave textile composite numerical examples, it is shown that the computational framework can accurately predict the stress-strain behaviour of RVEs with complicated geometries subjected to different loading scenarios. The developed computational framework is implemented in the authors' open-source FE software MOFEM; this has additional capabilities, including generalised RVE boundary conditions, hierarchic finite elements and optimisation for highperformance computing. The developed computational framework provides the nonlinear micro/meso-mechanical response at lamina level, which will be used subsequently to simulate FRP composites at both laminate and structure level. The proposed computational framework is valid only in the pre-peak region. Furthermore, no localization or damage is included in the matrix plasticity model. These limitation will be addressed in the future work.