A failure model for the analysis of cross-ply Non-Crimp Fabric (NCF) composites under in-plane loading: Experimental & numerical study

In the present work, a set of physically-based failure criteria are applied to a cross-ply Non- Crimp Fabric (NCF) composite. The proposed criteria account for the stitch reinforcement and the induced in- and out-of-plane fibre misalignment. A distinction between matrix and fibre failures in tension and compression is introduced, while the in-plane shear non-linear response is directly incorporated to the model as well. The constitutive model is associated with an energy-based damage mechanics method, which consists of five in-plane damage variables per layer level. Within this study, a novel approach to obtain the effective longitudinal stiffness of the composite, as a function of the in- and out-of-plane misalignment, is described. An experimental study characterises the internal structure and the mechanical response of the composite in tension, compression, in-plane shear and Mode I & II interlaminar fractures. In this paper, validation examples for the failure criteria, implemented in Abaqus/Explicit as a VUMAT subroutine, are presented. Tensile, in-plane shear and compressive responses are modelled at a coupon level with continuum shell elements and correlate well with the experimental data. Fibre misalignment, mainly out-of-plane, had a strong influence on the compressive modulus and strength of the composite. kinking and inter-laminar delamination was presented, which appeared to have a discrepancy of approximately 25% against test results. This work aims to present a holistic research on cross-ply NCF composites starting from an extensive experimental study to the numerical modelling. The experimental findings and observations are utilised to form a set of physically-based failure criteria, capable of predicting the mechanical response of cross-ply NCF composites under the main in-plane loadings and relate the applied loads to adequate damage mechanisms. To the knowledge of the authors, none of the existing failure models explicitly incorporates fibre undulation in both planes, but instead they are assigned on the elements as local coordinate systems. Hence, the novelty of this study includes the introduction of the in- and out-of-plane fibre waviness to the constitutive model as angular misalignments, without being represented as geometric feature. The proposed set of failure criteria is implemented in an existing energy-based damage mechanics model of modelling unidirectional thin laminates and woven composite materials which is convenient for the analysis of hybrid structures. to the layerwise approach of the model, the failure of the NCF composite can be simulated at each layer of an NCF blanket for different failure modes, while the option to model a homogenised blanket is maintained. Also, because of the plane-stress formulation introduced into the damage model, continuum shell elements are applied, which offer significant computational time savings when compared against solid element analyses. The limitations of the plane-stress approach, namely neglection of through-the-thickness stresses and constant transverse shear stresses, are acknowledged and accounted when post-processing the numerical results, although they did not significantly influence the examined case studies.


Introduction
Non-Crimp Fabrics (NCF) have been under development for some years; they consist of fibre tows, usually made of carbon or glass, which assumed to be straight and aligned, stitched together with a fibre either from polyester or aramid. A NCF blanket can be consisted of two (cross-ply), three (triaxial) or four (quadaxial) layers. Among their advantages when compared with other types of composites, NCF are characterised by reduced manufacturing cost [1] and a high level of drapeability [2]. The combination of their low production cost with non-autoclave liquid moulding methods, e.g. Resin Transfer Mould and Vacuum Infusion, results in an increased application of NCF composites to the aerospace, automotive and wind energy sectors [3].
In terms of the internal structure, NCF composites could be classified between the UD pre-preg taped (UDPT) and woven composites, exhibiting fibre waviness both in-and out-of-plane. The stitching parameters, such as stitch length & pitch, are responsible for the in-plane waviness while the out-of-plane fibre waviness is attributed to the imperfections occurring during the manufacturing process. The influence of the stitching pattern to the mechanical performance of NCF composites has been investigated by Asp et al. [4] & Vallons et al. [5] in quasi-static & fatigue conditions and interestingly, did not reveal any difference between NCF composites with different stitching properties.
The influence of the stacking sequence on the tensile stiffness of cross-ply NCF composites was studied by Edgren et al. [6] and Mattsson et al. [7], which identified a significant reduction in stiffness for [0 • /90 • /0 • /90 • ] S lay-up when compared with D. Gouskos and L. Iannucci
In-plane linear shear moduli of + 45 • & −45 • layers 13 , 23 Transverse shear moduli 12 In-plane shear modulus of fibre bundle 12  , Out-of-plane fibre misalignment initial, mean, maximum & additional rotation 12 , 13 In-plane & out-of-plane Poisson's ratios at layer level 12 , , Poisson's ratio of fibre, resin & polyester , , Density of fibre, resin & polyester 11 , 22 Normal stresses in the fibre and matrix directions at layer level 12 In-plane shear stress at layer level , , , In-plane fibre misalignment initial, mean, maximum & additional rotation  [8] highlighted the fibre waviness of the 0 • layers as a determinant factor for the longitudinal stiffness degradation and developed a master curve approach to predict it. Earlier, Cox and Dadkah [9] introduced an analytical knock-down factor methodology to model the axial stiffness degradation. However, none of these approaches explicitly DCB Double cantilever beam FI Failure index OM Optical microscope UDPT Unidirectional pre-preg tapes distinguished the influence of the in-and/or out-of-plane waviness on the longitudinal stiffness reduction, therefore in this work, a novel global knock-down factor is presented, composed by two individual sub-terms accounting the explicit influence of the inand out-of-plane waviness in the form of angular misalignment.
The tensile strength of carbon fibre NCF composites is generally lower than of UDPT in a range of 3% [10] to 40% [1]. Recently, Arteiro et al. [11] attempted to improve it by forming blocks of 0 • layers, however no significant improvement was noted. Similarly, the compressive strength of NCF composites has been reported by several researchers [12][13][14][15] to be 10%-50% lower when compared to UDPT and highly influenced by the out-of-plane fibre misalignment [16]. Interestingly, Wilhelmsson et al. [17] noted that the compressive strength is less susceptible to off-axis angles between 0 • and 10 • while it gradually drops at angles over 10 • . Failure mechanisms under compression were initially studied by Argon [18] and identified that the initial fibre misalignment with respect to the laminate and loading direction results to a matrix-shearing failure and thus kink band formation. Recently, Shipsha et al. [19] investigated the compressive strength and the damage mechanisms of cross-ply NCF composites for different off-axis angles and identified that the governing failure mechanism for all specimens was fibre kinking, with only exception the specimen with 45 • off-axis angle where the governing failure mechanisms were shearing cracks within fibre bundles and cracks at bundle-matrix interface.
Equally important and challenging is the numerical modelling of NCF composites. Two main approaches are used: a mesoscopic in which the constituents of the NCF composite are modelled by FE methods as a Representative Volume Element, incorporating the heterogeneity of the NCF composite. Such models have been presented in several studies [12,[20][21][22][23][24][25], each study with different assumptions, but able to predict the stiffness and strength of the NCF composite accurately, while the failure mechanisms were also successfully simulated. Nevertheless, this approach is computationally demanding and impractical for implementation in large scale structural components. In the second method, the modelling of the NCF composite is employed at the macroscopic level where classical laminate theory (CLT) is used for the analysis. In this approach the composite consisted of homogenised layers of different orientations assuming evenly distributed fibre tows, whilst the mechanical properties at layer level are obtained either from the rule of mixtures (ROM) or analytical multiscale methodologies [26,27]. The internal structure, namely fibre undulation, is introduced usually as a local coordinate system at an element level or as a sinusoidal function assuming a periodicity in fibre waviness.
Failure prediction, however, is far more complicated due to the different damage mechanisms occurring until failure onset, particularly in fibre kinking [19] and matrix failure. In comparison to the variety of failure criteria, developed for UDPT, limited research has been performed until now for the generation of failure criteria specifically for NCF composites. Edgren [28] assumed a linear interaction between the longitudinal compressive stress and shear stress to predict the strength of a layer. Although, the constitutive model achieved a good agreement with experimental results for quasi-isotropic NCF composites, it required a high number of input data for every layer of the configuration. A defect severity model to assess the longitudinal compressive strength of UD NCF composite, was proposed by Wilhelmsson et al. [29] while emphasis was given to the stochastic nature of fibre composites. The 2D FE implementation of the model predicted accurately the compressive strength and the formation of kink-bands. Molker with Wilhelmsson et al. [30] proposed physically-based orthotropic criteria for the transverse failure of UD NCF composites following a meso-micromechanical approach and predicted the orthotropic strength behaviour, distinguishing intrabundle and interbundle failures. While the present work was being developed, Shipsha et al. [31] studied the effect of the stacking sequence on strength and failure mechanisms for various cross-ply NCF lay-ups under compression. Also, a failure criterion for intra-laminar fibre kinking and inter-laminar delamination was presented, which appeared to have a discrepancy of approximately 25% against test results.
This work aims to present a holistic research on cross-ply NCF composites starting from an extensive experimental study to the numerical modelling. The experimental findings and observations are utilised to form a set of physically-based failure criteria, capable of predicting the mechanical response of cross-ply NCF composites under the main in-plane loadings and relate the applied loads to adequate damage mechanisms. To the knowledge of the authors, none of the existing failure models explicitly incorporates fibre undulation in both planes, but instead they are assigned on the elements as local coordinate systems. Hence, the novelty of this study includes the introduction of the in-and out-of-plane fibre waviness to the constitutive model as angular misalignments, without being represented as geometric feature. The proposed set of failure criteria is implemented in an existing energy-based damage mechanics model developed by Iannucci et al. [32,33] already capable of modelling unidirectional thin laminates and woven composite materials which is convenient for the analysis of hybrid structures. Thanks to the layerwise approach of the model, the failure of the NCF composite can be simulated at each layer of an NCF blanket for different failure modes, while the option to model a homogenised blanket is maintained. Also, because of the plane-stress formulation introduced into the damage model, continuum shell elements are applied, which offer significant computational time savings when compared against solid element analyses. The limitations of the plane-stress approach, namely neglection of through-the-thickness stresses and constant transverse shear stresses, are acknowledged and accounted when post-processing the numerical results, although they did not significantly influence the examined case studies.

NCF material system
The selected fabric for this project is a biaxial (+45 • /−45 • ) pillar stitched preform of Toray T700S fibres -12k filament yarn, weighing 300 g/m 2 [34], supplied by Cristex Composite Materials. The fabric is reinforced by polyester stitches of 8.0 g/m 2 areal mass and 8.0 tex linear density. The resin system is the MTM57 medium temperature, multi-purpose epoxy in the form of film, supplied by Solvay [35]. Table 1 summarises the mechanical properties for the fibres and the resin system as obtained from the suppliers. The mechanical properties of the polyester stitches were obtained by performing tensile tests, according to the ASTM D2256 [36] method after removing them carefully from the preform, in an Instron 5960 test machine. The yarns were constrained at a gauge length of 150 mm and tested at displacement rates varying from 100 to 300 mm/min. The stress-displacement response is shown in Fig. 1 and reveals rate dependence only within a very small range. The shear modulus, G Pes was estimated from the conversion formula of the homogeneous isotropic linear elastic materials [G = E/(2 ⋅(1 + ))], while the Poisson ratio, Pes was received from literature sources [37] equal to 0.25. The fracture toughness, G T Pes is characterised by the area under the stress-relative displacement curve, while it can also be derived via the integration, = ∫ ∞ 0 as suggested by Hertzberg et al. [38] and Miracle et al. [39] for unnotched specimens or samples loaded in tension until fracture. The mechanical properties of the polyester stitches as obtained from the tests are included in Table 1. With the assistance of a pycknometre the densities of the fibres, resin and stitches were found equal to = 1.7651, = 1.2028 & = 1.3925 g∕cm 3 respectively.

Laminate fabrication
Four laminates were manufactured with the Resin Film Infusion method, vacuum bagged to 135 • C at a ramp rate of 100 • C/h, under autoclave pressure of 7 bar and dwelled for 1 h before a non-controlled cool down to room temperature. The ultrasound C-Scanning method was implemented to evaluate the quality of the panels. Laminate properties, namely panel name, stacking sequence and panel thickness are summarised in Table 2.

Constituents volume fractions
A 60% fibre volume fraction for each laminate was the target, however, the actual volume fractions of the three constituents (T700S, MTM57 and polyester stitches) were determined with Thermogravimetric Analysis (ASTM D 3171-99 [40]) and optical microscope (OM). The calculation of the average fibre, resin & polyester stitches volume fractions of the laminates was based on [41] & [42] and found V F = 0.62 ± 1.3%, V M = 0.36 ± 1.0% & V Pes = 0.02 ± 0.5% respectively. Using the OM on the cross-sections of the laminates, the average fibre volume fraction within the fibre bundles, (V Fb ) was determined to be 0.85 with the average number of fibres per tow to be 11.7k according to the experimental observations.

Fibre misalignment characterisation
In previous studies [13,21] the in-& out-of-plane fibre waviness was assumed as periodic wave, however, in the present study, the vast majority of OM observations revealed fibre misalignment with only limited periodicity. Therefore, it was opted to be defined by recording the inclination angles instead.

Measurements of waviness in cross-ply laminate
In recent studies [15,43], advanced image analysis techniques were developed for measuring accurately the in-or out-of-plane fibre waviness, and are implemented in the present work as well. In Fig. 2, the procedure that was followed is presented. Every specimen was partitioned to sections in the area of the gauge length at a size varying from 20 to 50 mm depending on specimen's dimensions. The OM greyscale images were converted to binary images with the assistance of Matlab. Each section was divided to smaller cells, in order to discretise the curved fibres to straight lines and thus measure the angular misalignment. A spatial resolution of 50 pixels x 50 pixels is determined per cell, while three cells per layer thickness (0.135 mm) are assigned. Then with the use of Fast Fourier Transformation (fft2) the slopes of the fibres are obtained [15]. For each laminate a mean initial fibre misalignment angle, for out-of-plane and for in-plane, is calculated as: , where, n is the number of measured angles per laminate. The misalignment angles were determined for all individual layers per laminate. In total, 6000 measurements were obtained to define the in-and out-of plane fibre misalignment for all specimens. Table 2 summarises the mean and maximum in-and out-of-plane initial misalignment for each laminate. It can be identified that the in-plane misalignment is significantly lower than the out-of-plane due to the geometric constraints that the stitching reinforcement applies to the in-plane direction. Furthermore, in Fig. 3 the distribution of the recorded misalignment angles -& -for all laminates is best described by two Gaussian fitting curves with mean value ( ) 2.02 • and standard deviation ( ) 1.02 • for and = 3.88 • , = 2.28 • for respectively. The bar chart reveals that the scatter of is much narrower compared to . This can be attributed to the essential dependence of the in-plane misalignment to the stitching parameters such as the pitch length and spacing, which define some periodicity in the misalignment.

Mechanical testing
A series of both standard and non-standard mechanical testing at coupon level was carried out to characterise the in-plane behaviour of the NCF composites in tension, compression and shear. The interlaminar fracture properties were obtained through Mode I (DCB) & Mode II (4pt-ENF) tests. Ten specimens were tested for each loading case with the test results used later to validate the newly proposed failure criteria.

Specimen preparation
All specimens were waterjet-cut from the laminates and C-scanned to identify potential damage. The specimen end-tabs for tension, compression and in-plane shear were produced from woven glass fibre/epoxy composite, bonded with Araldite 2011 [44] resin-epoxy, while flatness and parallelism tolerances of the specimens were compliant to the test standards. The compression D. Gouskos and L. Iannucci   order to identify the effect of the in-plane misalignment to the compressive behaviour. For the DCB specimens, two aluminium end blocks with centre holes of 7.5 mm diameter were bonded at the end parts of each half-specimen with Araldite 2011. Both DCB & 4pt-ENF specimens were painted in white through-the-thickness at one side and a scale of 1 to 5 mm incrementation was marked to facilitate the view of the crack growth using the Immetrum optical strain camera (Fig. 4). Moisture removal techniques were followed for all specimens according to ASTM D 5229 [45]. Table 2 summarises the geometric aspects of the specimens.

Testing setup
The tensile tests were performed in accordance with ASTM 3039 [46], at an Instron 5980 test machine at a cross-head speed of 1.5 mm/min. For tension as well as shear test the GOM DIC camera with maximum resolution 4096 x 3000 pixels 2 was used to measure the axial and transverse, strain and displacement fields and detect any in-and out-of plane bending effects of the specimens, while GOM Correlate software was utilised to postprocess the test results. For both tests the Region Of Interest (ROI) was the gauge length of the specimen (145 x 25 mm 2 -2900 x 500 pixels 2 ), therefore it was finely speckle patterned. The pixel resolution was 0.05 mm whereas the face size (subset size) was set to 19 x 19 pixels 2 and point distance of 16 pixels was used meaning a 50% overlap, as suggested by Aramis for 6 Megapixel CCD cameras. To filter undesired data, a simple mean filter was used in the DIC commercial software Aramis GOM setting value 17. The GOM DIC camera was synchronised with the testing machine at 15 Hz. All previously described DIC parameters were verified after performing sample tests and compared against already validated data.    Table 3 summarises the average tensile properties of the specimens, namely longitudinal Young modulus (E xx T ), Poisson ratio ( xy ) and tensile strength (X T L ) and the respective coefficients of variation (CV ). The bending effects were significantly less than 10%, which is the premature bending threshold as defined by ASTM D 3039 standards.

Fractography
All tensile specimens exhibited a fibre-dominated pull-out failure (Fig. 5 left), with failure section perpendicular to the loading direction. However, the earliest event during tensile loading, has been identified as transverse cracking within the bundles of the 90 • layers, which has also been reported in [1,4,6]. The transverse cracks were developed at the free edges and spread inwards to the centre of the specimen [47] (white arrows in Fig. 5 right) while the stitches either debonded from the surrounding matrix resulting in matrix failures or acted as initiation sites for damage evolution.

Testing setup
The in-plane shear characterisation utilised the ASTM D 3518 [48] standard, where the specimens were loaded axially at a displacement rate of 2 mm/min. The same apparatus as for the tensile tests was used, while the GOM DIC camera was set at a 5 Hz frequency.

Results
The shear response of the axially loaded ±45 • lay-ups was highly non-linear beyond the initial linear stage ( = 0.005 in Fig. 6) due to the degradation of matrix. Both shear modulus, G xy and strength, results show high consistency as can be seen in Table 3 from the low values of CV. Bending effects on the specimens, in-and out-of-plane, were not identified from the DIC analysis. Table 3 Experimental results for tension and in-plane shear tests.

Fractography
The main failure mode identified in the specimens is scissoring, initiated at shear strain, = 0.04 prior to the ultimate failure ( = 0.055). As scissoring is defined the fibre reorientations due to pure shear deformation, where the intern fibre tow angle changes due to matrix driven failure (either as fibre breakage and/or fibre/matrix debonding) in the pure resin pockets while it is assumed that fibre tows do not slip at their pivoting points [49]. This behaviour significantly reduces laminate shear stiffness, prolongs shear strain failure [50], and is illustrated through the DIC observations and measurements in Fig. 7. From the DIC camera, it was detected that the relative angle between the + 45 • fibre tows and the −45 • for the above mentioned shear strains is reduced to 78 • from 90 • (angle between white dashed and solid lines). Moreover, shear failure is generated by inducing intralaminar shear within the NCF blankets and delamination between them ( Fig. 8-left). In addition to these failure modes, blanket/layer splitting occurs in the form of broom/brush-like failure modes ( Fig. 8-left). This is a common shear failure mode in materials with poor fibre/matrix interface strength like NCF composites which then results to stepped fracture of fibres ( Fig. 8-right).  Table 4 Compression test results.

Testing setup
The compression tests were performed in an Instron 5960 test machine at a cross-head speed of 1 mm/min, in accordance with the IC compression method [51]. Linear strain gauges (FLAB-2-11-1LJC-F) were attached at each face of the specimen to identify any bending effects during compression tests.

Results
The stress-strain curves for the different in-plane fibre misalignments are presented in Fig. 9 which exhibit a non-linear behaviour before the ultimate failure. This suggests the development of permanent damage which was validated by unloading a number of specimens prior to the ultimate failure. Comparing the specimens with the increased misaligned lay-up, it is noticed that the longitudinal stiffness is reduced by 4% and 7% for = 2 • & 4 • respectively whilst for 6 • , 8 • & 10 • a steeper reduction up to 35% is exhibited. This behaviour is not only attributed to the increased in-plane misalignment but also to the recorded higher values of for those specimens. As expected, the ultimate compressive failure stresses showed a similar trend, however it appears to be more susceptible to the fibre misalignment in both planes, since the reduction of strength reached higher values in terms of percentage, namely from 6% for = 2 • to as low as 40% for = 10 • . The test results are summarised in Table 4.

Fractography
SEM analysis revealed that the dominant failure mechanism for all specimens is fibre kinking as shown in Fig. 10. Kinkband formation was observed in the areas of the specimen with considerable crimp, whereas in regions with less crimp, intralaminar and interlaminar delamination was developed ( Fig. 10 -left). The representative macroscopic failures are clearly noticed as an angled crack through-the-thickness of the specimens accompanied by severe delaminations as a result of the shear driven brittle compressive failure of the 0 • & 90 • layers. As the shear component increases, failure starts to be influenced by intralaminar split initially of the 90 • layers and later of the 0 • layers. The 90 • layer splitting and the subsequent delamination between the 0 • /90 • layers are visible in Fig. 10 (right). By performing interrupted compression tests, it is apparent that the load-bearing fibres have kinked into these localised delaminations.

Testing setup
ASTM D 5528 [52] standard was used for Mode I DCB interlaminar fracture toughness test. The tests were performed in an Instron 5960 test machine at a tensile cross-head speed of 1 mm/min. The crack growth was monitored by the Immetrum optical strain camera with pixel resolution of 0.05 mm which was synchronised at 3 Hz with the test machine. The load from the cross-head and the optical displacement curves were plotted purely from the pin locations.

Results
A load-displacement curve is shown in Fig. 11(a), and reveal an initial increase in the load when the crack initiates, which is followed by successive load drops and crack arrests. The interlaminar fracture toughness, G IC , was calculated both by the Modified Beam Theory (MBT) and the Modified Compliance Calibration Method (MCCM) ( Fig. 11(b)). The results are interpreted as fracture toughness of initiation and propagation (Table 5) as required for the numerical modelling.

Fractography
The fracture surfaces of a DCB specimen are presented in Fig. 12. Fracture extends along the fibres (as indicated by the white dashed arrows) and then spreads into the surrounding matrix. Stitches do advance crimp development on the tows generating an undulating fracture surface which increases the fracture toughness and contributes to a slip-stick fracture process, as has also been observed in [53]. At microscopic level, fibre splitting was also detected ( Fig. 12 -left) through SEM in areas with high degree of crimp, which also enhances fracture toughness. Furthermore, a stitch bridging event at the crack was also observed, during which stitches are elastically stretched providing resistant forces in crack opening resulting in higher propagation load, until their failure or pull-out occurs, where load drop is observed. This mechanism enlarges the area surrounded by the load-displacement curve and thus enhances the fracture toughness (Fig. 12 right).

Testing setup
To obtain the Mode II interlaminar properties tests, the non-standard 'Martin and Davidson' [54] method was selected instead of the standard 3pt-ENF, since in the first a method it has been reported that a stable crack growth under displacement control can D. Gouskos and L. Iannucci  Mode II results of a tested specimen. Table 5 Interlaminar fracture toughness of initiation and propagation for DCB and 4pt-ENF tests. be achieved [55]. Furthermore, during trial 3pt-ENF tests which were performed by the authors, excessive matrix damage on the top surface of the specimens at the interface with the loading pin was observed. The test was also performed in an Instron 5960 test machine, with the set-up as shown in Fig. 4 (right) whereas the loading pins were free to rotate about a central pin across the specimen width. The test was loaded in compression at a rate of 1 mm/min while the Immetrum optical strain camera was synchronised at 3 Hz with the test machine and set-up as for the DCB test.

Results
A load-displacement curve obtained from the tests is shown in Fig. 13(a), exhibiting some unstable crack propagation, leading to a sudden drop and the crack arrest at the next stitch line. The initial increase of the load with crack growth is mainly due to the development of bridging stitch thread zone, resulting in increased inter-laminar shear toughness, G IIC . The data reduction method implemented in 4pt-ENF test, is the Martin and Davidson Compliance Calibration method (CCM) [54] (Fig. 13(b)). Table 5 summarises the initiation and propagation mean values respectively for all specimens.

Fractography
As for 4pt-ENF tests, stitches significantly contribute in fracture toughness enhancement, due to the development of bridging stitch zone. At the resin rich areas, large cusps are often apparent, associated with the stitches. Determining the crack growth direction of the interlaminar shear fracture surfaces from visual observation, was challenging, nevertheless the SEM indicated the presence of surface debris and the crack growth direction as pointed by the white dashed arrows (Fig. 14). Matrix fracture initiates at the fibre/matrix debonding and then extends as scarps into the matrix. Cusps are adjacent to the fibre plane, retained to the surface of the one half specimen (Fig. 14 right), while the other half contains only scallops (Fig. 14 left).

Material properties -mesomechanical approach
In the present work, the material properties calculations were performed at layer level and were inspired by the methodology proposed by Olsson, Marklund & Asp. [26,27] and ROM. The stitches were included at a 45 • ( ) orientation to the fibre direction ( Fig. 15) due to the architecture of the current NCF, whereas in the previous studies stitching yarns were assumed parallel & perpendicular to the fibres. The longitudinal modulus E 11 is obtained by a modified ROM (Eq. (1)) while the transverse moduli from the Halpin-Tsai relations [56], enhanced to account the angular component of the stitches (Eq. (2)): The effect of waviness in fibre bundles is not accounted at this stage in Eq. (1), because it is later introduced in Section 4.2 of this document in the form of a knock-down factor. Moreover, matrix stiffness (Eq. (2)) has not been reported in previous studies to be influenced by fibre waviness, therefore deliberately fibre misalignment is not accounted.
where, E Pes 45  where the stitches are aligned perpendicular to the shear load. Eq. (3) estimates the shear modulus for both orientations, however for the −45 • layer the term of enhanced stiffness in the angled bracket is omitted.
where, j = + 45 • or −45 • for the respective layer orientations, , with = 26 as advised by Asp et al. [27]. 12 is the shear modulus of the fibre bundle. The calculation of Poisson ratio of the layer can be found in Marklund et al. [26]. The strength prediction of an NCF layer in the fibre direction for tension and compression is estimated by Eq. (4), and it is assumed to be controlled by the fibre failure during loading.  where, i = T or C for tensile or compressive loading, X i b is the strength of the bundle and E b 11 the elastic modulus of the fibre bundle, the elastic and strength properties of the bundles were obtained from the composite cylinder assemblage [57,58] and together with the predicted material properties of T700S/MTM57 system at layer level are presented in Table 6.
With regards to the transverse and shear strength estimation, it was decided to use the in-situ transverse and shear strengths as developed by Camanho et al. [59] & Dvorak et al. [60]. Although, they were originally developed for UD composites, the fundamental parameters that define the in-situ strengths are identical for NCF too, namely the thickness of the ply or layer, the lay-up of the laminate and the value of the fracture toughness. In recent studies for NCF [61], the in-situ strength formula for UD composites was successfully applied, therefore is followed in the present work too. The analytical formulations in transverse tension, , compression, , longitudinal, and transverse shear, can be found in [59,60,62]. At this point, it has to be stated that Dvorak & Camanho's formulation originally requires the intralaminar fracture toughness values, however the authors consider that it is a fair assumption to use the interlaminar fracture values instead, since as suggested by Pinho et al. [63] and Catalanotti [62] in the absence of experimental data and due to the similarity of the fracture modes for the same material systems, the intralaminar fracture toughness can be approximated with the interlaminar for Mode I and Mode II. The in-situ strengths were calculated for layer thickness 0.13 mm and , which defines the non-linearity of the shear stress-strain relation, has been obtained according to [59] equal to 2.8e−8 MPa -3 . Table 7 presents the calculated in-situ strengths.

Effective longitudinal modulus
The knock-down factor (KD) approach is introduced to estimate the effective longitudinal modulus, 11 . of a fibre bundle within a NCF layer. Inspired by Cox and Dadkhak's [9] approach which assumes a periodic fibre waviness, the present study suggests the KD, (k), to be employed at layer level, composed by two individual KD subfactors due to the in-plane ( ) and out-of-plane ( ) fibre misalignment. The initial misalignment angles, -and -are introduced instead of the wave properties (wavelength, and amplitude, A). The ratio A/ can be replaced by tan or , tan( , ) = 2 ⋅ ⋅ ∕ as illustrated in Fig. 16. Hence, Eqs. (5) can be developed, where the denominator equals to 1 for straight tows, while the magnitude within the square brackets corresponds to the anisotropy factor.
Thus, the effective Young modulus in the longitudinal direction equals: E 11 Eff = k ⋅ E 11 . Fig. 17 depicts the evolution of k as and increase.

Failure criteria for NCF composites
The set of failure criteria for NCF composites proposed in this work is based on a physical model for each failure mode, namely tension & compression in fibre & matrix directions. Furthermore, the structural properties of the NCF composites, i.e. induced fibre misalignment, are accounted, since they have a major effect in strength in specific loading directions at layer level. The failure modes are developed for in-plane loading conditions but they can be extended for a 3D stress state too. Failure is predicted when the failure indices, equal to 1.

Fibre tensile failure
The maximum stress criterion is used to predict damage onset: where FI FT is the failure index associated with fibre tensile failure and 11 is the tensile stress in the longitudinal direction.

Matrix tensile failure
The matrix tensile failure index, FI MT of a NCF layer is described by the quadratic interaction criterion of LaRC03 [64], which couples the shear stress, 12 with the transverse tensile stress, 22 .

Fibre compression failure
The proposed failure model for fibre compression is based on Argon's hypothesis [18] and on the subsequent developments of Davila et al. [64] and Pinho et al. [65]. Both the in-and out-of-plane fibre misalignments are included in the failure criterion. A coordinate system 1-2-3 is aligned with the material axes of the composite: 1 = fibre direction, 2 & 3 transverse directions. The fibre kink band, is described in Fig. 18 and it is obtained numerically at a range 0 • ≤ ≤ 180 • . The coordinate system associated with the kink plane is 1 ( ) -2 ( ) -3 ( ) and the resulting new stress tensor is:   where, = or . Davila et al. [64] originally developed the additional rotation component at one plane. In this work, the additional rotation parameter is extended to include the fibre misalignment components of both planes 1-2 & 1-3 (Eqs. (9)). Note, that since two different planes are referred, the corresponding material properties need to be used i.e. G 12 & G 13 .
The stress tensor of the misaligned coordinate system is then: The resultant rotation R m , in the where, & are the internal material frictions ( = −1∕ tan 2 0 , = ⋅ ∕ ) and is the transverse tensile strength in direction 3. The relevant tractions for the misaligned stresses are summarised as:

Matrix compression failure
The failure in transverse compression, FI MC is predicted by the Mohr-Coulomb failure criterion. However, the in-situ strengths were introduced to the criterion as was recently proposed by Catalanotti et al. [62] and thus the criterion is described by Eq. (12).
with the effective stresses in transverse and longitudinal directions to be: The fracture angle, 0 usually ranges between 0 • and 53 • with the extreme values referring to pure shear loading and pure compression respectively.

Continuum damage mechanics method
Extensive details about the analytical formulation of the energy-based damage model can be found in [32,33], however, in the following subsections a brief overview of the model is presented.

Effective stress-strain relationship
This damage mechanics approach introduces five damage variables for in-plane damage per NCF layer (one cross-ply NCF blanket is made of a 0 • and a 90 • layer). Damage modes, are subdivided according to the stresses acting in respect to the local material orientation: i = 1: fibre tension, 2: matrix tension, 3: in-plane shear, 4: fibre compression & 5: matrix compression. The degradation of 11 . , 22 & 12 . is defined as shown in Eqs. (13): where d = 0 represents a virgin material and d = 1 the fully damaged material. Also, the Poisson's ratios are degraded in a similar manner to the Young moduli, after satisfying the following relationship:

Damage evolution
Damage propagation is described in terms of the directly related true strains whereas the direct stress components, in tension and compression, are expressed by a bi-linear stress-strain behaviour. Hence, the stress in the damaged material at each current strain level > 0, is defined in Fig. 20 as, OABC for tension and OABDE for compression. However, in compression, BD relates to the residual strength (X CR ) of the kink and matrix crush zone until failure occurs at a strain , . Additionally, E C is defined as the energy dissipated due to 'plastic' deformation of the crushed zone. The bi-linear section of the damage evolution curve is defined as: wherėis the strain-rate enhancement. The shear damage variable, 3 is expressed through Eq. (15):

In-plane shear component
The non-linear in-plane shear response was described in Section 3 and it is modelled in incremental form where the degree of non-linearity and G Ult the ultimate shear modulus, both obtained from the in-plane shear tests. Hence the shear modulus is defined by Eq. (16) as the greater of 12 . and the non-linear transition modulus: where is the shear strain (-) and 0 is the strain at the linear to non-linear transition (-).

Smearing methodology
A smearing cracking methodology is implemented to model fracture, in which the stored energy of a unit volume is related to the fracture toughness, to alleviate mesh dependency and overcome strain localisation. This method requires an additional length parameter in order to achieve a constant energy release per unit area of crack generated [66,67] regardless the element dimensions. The limitations of this approach are acknowledged by the authors as described by other studies [68,69]. For instance, the smeared formulation is effectively applied in structured meshes which are not recommended for complex geometries with discontinuities, viscous regularisation is used to obtain mesh objectivity [68] while crack growth direction can only be parallel to one of the element edges, which imposes limitations to predict damage evolution in arbitrary directions [70,71]. Though, in the present study the simulations are performed on rectangular specimens, which can be modelled effectively with structured meshes, and the loading scenarios are not expected to form complex failure patterns. The total specific energy, E f can be expressed as the sum of the elastic energy, E e and the propagation energy, E p . The specific energy dissipated is the area under the stress-strain curve and is defined as: = 1∕2 ⋅ 0, ⋅ , , the maximum element size that can be used in the mesh is defined by inequality < 2 ⋅ 0, ⋅ 0, .

Interface modelling
The cohesive surface approach is used to model the cohesive zones (CZ), applying a bilinear traction-separation law. The interface stiffness for Mode I, K n is defined by Turon et al. [72], ≥ ⋅ 33 ∕ℎ where, h is the thickness of the arm of Mode I specimen and is a bulk coefficient, which ensures that the compliance of the bulk material is much larger than the initial compliance of the cohesive surface. After performing a series of numerical simulations of DCB and 4pt-ENF tests a value equal to 50 was selected. The interface stiffness for Mode II follows the proposal of Iannucci & Willows's [32], = 2 ⋅ 23 ∕ , where G 23 is the shear modulus, which for numerical applications is assumed to be equivalent to the values attributed to the matrix. The magnitude e corresponds to the resin interface thickness between the NCF blankets. The initiation of delamination is modelled by a stress-based quadratic power law (Eq. (18)): where, , and are the normal and shear traction components respectively, while is the normal & , are the shear interface strengths. Because both Mode II and III are matrix dominated interlaminar shear failures in two orthogonal directions, where tangential separation of the blankets occurs and since matrix is treated in this study as a transversely isotropic material, it is assumed that = . The interlaminar shear strengths and are obtained by Turon's proposal [72]: ∕ =  element length size, should also satisfy the requirement of . In the analysis presented herein, the damage evolution and mixed-mode behaviour is described by the quadratic power law, as expressed by Turon [72]: where, , and are the energy release rates, assuming = . The power exponent, normally has values between 1 and 2. A value of 1.5 for the analyses carried out in this work was considered acceptable.

Finite element implementation
The analytical formulations and the constitutive failure criteria presented in Section 4 were introduced to Simulia Abaqus/Explicit as a VUMAT material model and validated against the test results presented in Section 3. Due to the non-linearities (e.g. in-plane shear, matrix failure modes) of the material constitutive formulation as well as due to the modelling of contacts (i.e. cohesive surfaces), where any change in contact is also extremely non-linear, the explicit analysis is opted so that a convergence solution is ensured. Continuum shell elements with reduced integration (SC8R) were utilised to model the NCF blankets of the specimens, while each layer of the NCF blanket was modelled with 3 integration points through-the-thickness. The fibres within the layers were homogenised and assumed to be evenly distributed.
Moreover, fibre misalignment is not represented as a geometric feature, but instead, it is introduced at layer level directly to the material card based on the OM observations for the tensile, compression & in-plane shear specimens. In particular, the regions on the specimens in which similar fibre misalignment was observed, are grouped together to specific partitions of the FE model. There is no limitation regarding the number of partitions, as it solely depends to the heterogeneity of fibre misalignment within the specimen. An advantage of this approach is that the fibre misalignment does not depend on the size of mesh since it is assigned to the geometric partitions and not to each element individually (as local coordinate system), a labour intensive procedure particularly in larger structural components. With regards to the in-plane fibre misalignment, , the values assigned to the top and bottom layers of the specimen are obtained directly from the OM, whereas for the interior layers, the average value of between the top and bottom layers was estimated. On the other hand, the out-of-plane fibre misalignment, is assigned at each layer individually. Due to the applied high loading rates (i.e. 3 mm/s) in the FE analyses, Rayleigh damping was introduced to minimise the generation of unreal modes and to ensure that the kinetic energy/internal energy ratio of the system will not exceed the order of 0.1. Concerning the correlation between the numerical and test results, an up to 5% discrepancy was considered acceptable. Also, no artifice was used to initiate failure at any specific location or along any certain direction.

Material properties for FEA application
The constitutive failure model requires the input of certain parameters for proceeding with FE analysis. Fracture plane angle for pure compression, 0 is found to be 50 • based on fractography. After one element in-plane shear simulations, the following magnitudes have been verified:  Table 8. Also, Table 8 summarises the interlaminar properties, where it can be observed that the selected , values are closer to the propagation values than the initiation. The authors performed trial numerical simulations in compression, using the initiation values of interlaminar fracture toughness ( , ). It was noted that premature decohesion occurred within the blankets at the midsection of the specimen which was not validated by the experimental findings. On the other hand, when introducing values close to propagation fracture toughness -for crack growth = 8 mm for DCB ( Fig. 11(b)) and 4pt-ENF ( Fig. 13(b)) tests, equivalent to compression specimen gauge length-, the authors identified that interlaminar failure was in agreement with the test observations. Therefore, since the numerical model requires a single fracture toughness value, it was decided to use the corresponding values of interlaminar fracture toughness for 8 mm of crack length. Table 9 presents the calculated interlaminar strengths for the different values of cohesive zone lengths modelled with four elements.

Tension loading
The numerical modelling of the ([0 • /90 • ] 4S ) tension specimen was confined to the gauge length in order to maintain the computational time within reasonable limits.

Boundary conditions
In terms of boundary conditions, one end of the specimen is constrained against translation longitudinally & transversely (U x = U y = U z = 0), while at the opposite free edge a tensile displacement rate, F of 3 mm/s is applied (Fig. 21).

Analysis parameters
A sensitivity analysis is employed to validate the mesh objectivity of the damage model. Four different in-plane mesh densities were examined of 8 mm, 4 mm, 2 mm and 1 mm, assigned with 1, 2, 4 & 8 elements through the thickness of the specimen respectively. In this analysis, no interface properties were implemented between the NCF blankets since the influence of possible delamination was not considered to be essential. The values of and introduced in this simulation are presented in Appendix A for each section and layer of the specimen (Fig. 21). Also a second FE analysis with a 1x1 mm mesh density is performed, in which the average values for misalignment are used across all layers & sections namely, = 2.1 • , = 3.4 • and compared against the detailed configurations.

Results
In Fig. 22 the numerical models are compared directly to the experimental testing through the stress-strain response, which reveals a mesh objective behaviour. The detailed configured FE models correlated very well with the experimental data, both in capturing the tensile stiffness and predicting the ultimate stress and strain with a less than 3% discrepancy. The numerical model with the average values of fibre misalignments, appears to overestimate the stiffness of the tested specimen by 6% which is almost double than the detailed configurations, but still within the range of acceptance. As expected though, the simplified FE model, captured accurately the tensile strength of the specimen, since the failure criteria in fibre and matrix tension are not affected by fibre misalignment. In Figs. 23 & 24, the failure indices are presented for the different mesh densities and correspond to the fracture band of the 0 • and 90 • layers at the midsection of the specimen respectively. The elements failed close to the load introduction points, probably because the end-tabs of the specimen are not included to the numerical modelling. An important parameter that could be associated with the location of fracture is that the gripping force on the end-tabs is replaced by a concentrated displacement. Nevertheless, the failure morphology correlates well with the fractured areas of the actually tested specimen as it is illustrated in Fig. 24. The fracture band is perpendicular to the loading direction and extended along the width of the specimen as exactly exhibited in the actual failed specimens.

In-plane shear loading
Due to the stable crack propagation during the in-plane shear test, which occasionally is initiated from the end tab-region, it was decided to numerically model a full specimen, including the end-tabs, with the lay-up [

Boundary conditions
The boundary conditions (Fig. 25) differ from the tension FE model because of the end-tabs modelling. Longitudinal (x-axis) translational constraint (U x = 0) is applied at the end-tabs opposite to the tensile displacement, F (3 mm/s). The gripping force on the end-tabs, generated by the hydraulic grips, is commonly assumed to be 10-15% larger than the maximum axial test load to be applied [74], which for the in-plane shear tests of this study corresponds to 12 kN. Hence the gripping force is estimated to be 13.8 kN per end-tab, resulting in a pressure load, P of 11 N/mm 2 based on the area of the end-tabs (25x50 mm 2 ). The translational constraints in the -direction represent the reaction at the opposite grip of the fixture generated by the applied pressure. Translational constrains across the specimen (U y = 0) are applied to both faces of the end tabs.

Analysis parameters
A mesh objectivity analysis is employed for the shear FE model too, applying the same in-plane mesh densities as for tension FE model and through-the-thickness discretisation. The values of and for each section and layer of the specimen (Fig. 25) are presented in Appendix A. Like in tension, a second FE analysis with a 1 × 1 mm mesh density is performed, in which the average values for misalignment are used across all layers & sections namely, = 1.8 • , = 3.4 • and compared against the detailed configurations.

Results
In Fig. 26, the shear stress-strain curves are shown for the different mesh densities compared with the test results. The numerical models predicted the experimental behaviour well. Somewhat higher stresses appear in the non-linear transition for the FE models meaning that the uniformly absorbed energy over all the specimen in that region is higher than in shear test. The stress-strain response of the FE model with the averaged values of misalignment is deliberately offset by 0.005 to avoid confusion with the other curves. It is apparent that the shear behaviour of that model is identical to the detailed configuration. Neither the shear modulus nor the shear failure parameters (stress & strain) changed. That was expected since in-plane shear is governed by a matrix dominated failure which is not affected by fibre waviness.
In terms of the crack prediction and failure morphology of the FE model, in Fig. 27 is noticed that only the FE case with 8 × 8 mm element size, did not capture successfully the crack morphology, although the respective stress-strain plot correlated very well with the experimental data. Concerning the simulations of the other three element sizes, shear failure is predicted as a localised crack and appears not to be affected by the size of the element, while they agree well with the failure pattern of the failed specimen.
A comparison between the shear strain maps of FEA and DIC at failure is illustrated in Fig. 28, which reveals that the FEA contours for the 4, 2 & 1 mm element sizes showed good agreement with the DIC observations and predicted realistically the failure pattern of the specimen. On the other hand, the 8 mm element configuration predicted the location of the failure but not the morphology of the localised strain. The evolution of the shear damage variable at four discrete points of the stress-strain curve is shown in Fig. 29. The successful correlation between the numerical and experimental results is obtained mainly due to the scissoring effect that was taken into account in the non-linear shear analytical formulation.
With regards to the prediction of delamination, in Fig. 30, the interface failure at the midplane of the FEA model (1 × 1 mm element size) for time step just prior to failure is illustrated and reveals that the delaminated areas are predicted to be close to the end-tabs and around the crack. The numerical model is compared with the actual failed specimen, and it is shown that there is a good phenomenological correlation with regards to the delaminated regions.

Compression loading
The numerical analysis conducted for the compression test, includes the modelling of the [0 • /90 • ] 4S lay-up configuration and the cases with deliberate in-plane misaligned blankets ( = 2 • , 4 • , 6 • , 8 • & 10 • ) and the corresponding . Due to the nature of the failure modes (i.e. delamination, possible under end-tab failure) during compression test, it was decided to proceed with a full specimen model, including the end-tabs and cohesive zones.

Boundary conditions
Boundary conditions were modelled according to the Imperial College compression test method [51] fixture and they are illustrated in Fig. 31. Longitudinal translational constraint (U x = 0) is applied at the end, opposite to the applied load, F (2 mm/s). A pressure load, P of 36 MPa is applied on each of the top end-tabs, simulating the clamping force created by two M10 bolts at each grip block of the fixture. The clamping pressure is obtained using the formula correlating the bolt clamping force, L with the tightening torque, T (15 Nm), = ⋅ ⋅ [75], where D, is the nominal bolt diameter and K, the torque-friction coefficient equal to 0.20. The translational constraints in the -axis correspond to the reaction at the wall of the fixture generated by the applied D. Gouskos and L. Iannucci   pressure. Translational constrains across the specimen (U y = 0) are applied to the one side of the specimen with the opposite side to be set free for translation in the y-axis, so that the xy-plane Poisson effect is accounted.

Analysis parameters
As for the previous numerical analyses, an objectivity study is performed for the case of [0 • /90 • ] 4S without deliberate in-plane misalignment, = 0 • , involving mesh density variation and different number of CZs. Four cases are simulated: 2 mm mesh density: (i) without CZ and (ii) using 3 CZs; 1 mm mesh density with: (iii) 7 and (iv) 15 CZs. The CZs in all related cases were symmetrically allocated through the thickness of the specimen. The calculated interface properties in Section 6.2 were used for modelling the CZs for the respective element sizes. Moreover, simulations were carried out for the compression tests with ( ) equal to 2 • , 4 • , 6 • , 8 • & 10 • , and compared against FE analyses without the novel terms of Eq. (11). Additionally, simulations of models with average values for misalignment across all layers and regions were performed. For all these cases the 1 mm element size with 7 CZs modelling approach was applied and 8 elements through-the-thickness were implemented, one for each NCF blanket. Similarly, the CZ between the end-tabs and the outer NCF blankets are modelled according to Huntsman Araldyte 2011 [44] adhesive mechanical properties. The end-tabs were modelled with continuum shell elements (SC8R) introducing the elastic mechanical properties of Glass-Fibre [76].  The values of and introduced in the numerical models, for each section and layer of the specimen (Fig. 31) as well as the average values of misalignment, can be found in Appendix A.

Results
The stress-strain response of the sensitivity analysis is presented in Fig. 32 and is in good agreement with the experimental data for the different configurations. In particular, the linear behaviour of the tested specimen was predicted accurately, proving that the applied stiffness knock-down factor performed well. For strains of 0.004 onwards, a discrepancy in the non-linear region is noted between the FE results and the experimental behaviour; this could be attributed either to manufacturing imperfections within the specimen such as voids, or to additional stiffness degradation of the 0 • layers if higher values of fibre misalignment exist within the specimen which could not be observed or named optically. Ultimate failure was overestimated by the FE models, ranging from 10%, for the 2 mm mesh density without CZs configuration, to 2.3% when 1 mm mesh density with 15 CZs was applied. The undulations of the plots in Fig. 32 for strain values higher than 0.006 are due to the mass scaling of the FE analyses.
In Fig. 33, the cohesive zone failure index of the FE models for 3, 7 & 15 CZs is presented. The illustrations were taken at strain equal to 0.0095 where interface failure index reached unity, prior to ultimate failure. A similar delamination pattern is revealed for all three configurations. In the 3 CZ model, decohesion was initiated from the two outer NCF blankets at the 90 • /0 • interface ( Fig. 33-left, point A) and then propagated at the midplane of the specimen (Fig. 33-left, point B), therefore ultimate failure did not occur simultaneously with delamination.
For the 7 & 15 CZs configurations, unlike the 3 CZ model, the amount of separation observed is not expanded along the gauge length (8 mm) of the specimen (Fig. 33 middle & right) but it is limited to 5 mm and 4 mm respectively. In particular, in the 15 CZs model it is noted that the interlaminar failures are further localised to the region where fibre kinking and matrix failure occurs. Moreover, decohesion due to peel is identified between the outer layers of the specimen (Fig. 33-middle, point C & Fig. 33-right, point E). Interlaminar shear failure is again predicted at the mid-section of the specimens and the neighbouring CZs ( Fig. 33-middle, point D & Fig. 33-right, point F). The discrepancy between the different number of CZ models in predicting interlaminar failure, can be attributed to the fact that the relative slippage between the layers is smeared between the NCF layers through the thickness of the specimen for a high number of CZs, which is expressed as a slight non-linearity in the stress-strain diagrams, whereas for FE models with fewer CZs, this relative motion only occurs to the first available CZ, concentrating higher amount of failure energy. Comparing the delaminated predicted areas of the 15 CZ model with the experimental observation of Fig. 10 a good qualitative agreement can be observed.
The damage sequence of the four FE models is illustrated in Fig. 34 (a, b, c, d), which depict that the shear driven intralaminar failure of the 90 • layers was predicted to be the first damage event, initiated at the midsection of the specimen and gradually propagated to the 90 • layers towards the surface of the specimen. The 0 • layers exhibited fibre kinking failure which led to the ultimate compressive damage as illustrated in Fig. 34 (e, f, g, h). The maximum compressive stress is reached and the kink-band is formed as all support is lost from the surrounding matrix material (decohesion) resulting in a Mode I delamination.
In Fig. 35, the numerical results for the compressive strength, X C and longitudinal modulus, E xx for different values of & are presented and compared against the test results. The results of the simplified numerical modelling (averaged values of fibre misalignment) are described with indication ''Avg'', whilst the detailed numerical models with specifically assigned fibre misalignment are indicated with the subscript ''Spec''. In principle, the FEA predictions for ''Spec'' simulations appeared to be slightly more optimistic for both magnitudes when compared with the mean experimental results but still well within the margin of the experimentally recorded values. Specifically, E xx was overpredicted in the range of 2.1-3.5% whilst X C discrepancy varied between 2.5-4.2%. This was expected since a significant amount of uncertainty exists in the internal structure of the specimens not only D. Gouskos and L. Iannucci  in terms of in-or out-of-plane fibre misalignment but also as voids or other manufacturing defects, not accounted by the present constitutive model. Moreover, instead of the maximum values of & , the mean values were introduced. As a consequence of this approach, possible premature failures of the actual specimen in those areas where the maximum fibre misalignment existed, could not be predicted. Notwithstanding the aforementioned observations, the FEA models captured successfully the trend with increasing of both the elastic modulus and compressive strength in compression. Concerning the ''Avg'' numerical model predictions, compressive strength and modulus are overestimated by 2-4.3% against the detailed numerical model. An explanation for this outcome is that the homogenisation of the angles eliminates the effect of excessive fibre waviness in the sections and layers of the specimen which would result in locally reduced stiffness and strength and thus to early failure of certain elements. Also, it worth mentioning that with the increasing values of in-plane misalignment, the deviation in stiffness and strength prediction is increasing too, starting from 2.2% for 0 • in-plane misalignment up to 4.7 • for 10 • . Although, the discrepancy is not prohibitively high when compared against the test results and the detailed numerical model, the simplified configuration lacks the capability to predict localised damage due to fibre undulation and forecast the failure sequence.
The contribution of each term of Eq. (11) to trigger kinking failure under pure compression ( 11 < 0), was also investigated, with particular focus on the newly added terms C and E. A side-by-side comparison was performed between numerical models with full version of Eq. (11) (terms C & E included) and without the novel terms. In the absence of terms C & E the predictive compressive strengths for varying in-and out-of-plane fibre misalignment is essentially higher by 12 to 20% than when these two terms are included as shown in Fig. 36. Specifically, the deviation becomes higher with increased fibre misalignment. This discrepancy is mainly attributed to the absence of term C, which is susceptible mainly to and the longitudinal compressive stress ( 11 < 0), whereas term E is more sensitive to and transverse tensile stress ( 22 > 0). The significance of each term of Eq. (11) is depicted in Fig. 37, in which the evolution of each term is quantified for different and values for both configurations. Term A appears to be proportional to , while it has a small influence in damage initiation (0-0.033) since it is primarily dependent on the transversely applied loads. It is apparent though that for the configuration without terms C & E, term A is more essential by obtaining almost twice the values of the full Eq. (11) criterion. On the other hand, term B and the newly added term C appeared to be the dominant contributors so that the failure index in fibre kinking reaches unity. At this stage, it becomes noticeable that the absence of term C in the simplified version causes the overestimated compressive strength, because although term B reaches unity, that occurs for higher compressive stress. Term D is proportional to both & , hence its contribution to the failure index is negligible for ranging from 0 • to 4 • (0-0.019), while for angles from 6 • to 10 • it gradually obtains stronger impact (0.06-0.21). As in term A, so in term D, the index values are twice as much for the simplified version of Eq. (11), which implies that the contribution of term E is D. Gouskos and L. Iannucci  smeared within terms A and D. The newly added term E due to the limited change of exhibits a relatively constant behaviour and, as expected, since no external transverse loads are applied, the influence to the failure index is low. Thus, term E does not crucially affect the prediction capabilities of the simplified version of Eq. (11). The conclusion from this comparative study, highlights the significance of primarily term C in Eq. (11) under pure compression, which undoubtedly improves the prediction capability of the failure criterion.

Discussion and conclusions
In the present work, a set of physically-based failure criteria for the prediction of the quasi-static in-plane mechanical behaviour of cross-ply NCF composites was developed. The failure criteria were formulated in an energy-based damage mechanics model, implemented in Abaqus/Explicit as a user-defined VUMAT subroutine. Experimental work was performed to characterise the material system and the internal structure of the composite.
The novelty in the present proposal regards the introduction of the in-and out-of-plane fibre misalignment to the failure criteria. Shear non-linearity and scissoring effect were also considered by the model. In this study, a knock-down factor composed by two separate sub-terms in function of in-and out-of-plane fibre misalignment to describe the longitudinal stiffness degradation was proposed.
The failure criteria were validated under tension, in-plane shear and compression loading at specimen level, against various inand out-of-plane fibre misalignment combinations implemented in detailed and simplified numerical models. Failure was predicted at D. Gouskos and L. Iannucci layer level assuming homogenised layers, without employing geometric features of the reinforcement to the model. The constitutive model was able to predict successfully the tensile behaviour as well as the damage band formation of the specimen. Similarly, for the in-plane shear loading case, the numerical results were in good agreement with experiments both in terms of predicting the transition of linear to non-linear response, the formation of shear bands and the ultimate shear strength. The compression numerical analyses for several combinations of in-and out-of-plane fibre misalignment overestimated the longitudinal modulus and compressive strength by an average factor of 4%, captured accurately the damage sequence of the specimen and revealed the strong influence of the newly added terms of the fibre kinking criterion as well as the higher influence of the out-of-plane fibre misalignment to fibre kinking than the in-plane.
Whilst simplified numerical models with average fibre misalignment values exhibited reasonably good agreement with the test results for tension and in-plane shear cases, that was not the case for compression which is more sensitive to fibre misalignment and strong discrepancies against detailed numerical modelling and experimental results were noted.
However, an uncertainty factor was identified with regards to the internal structure of the composite since it is not feasible to capture the entire amount of fibre misalignments of a specimen, therefore in a future study, a probabilistic approach might be advisable. Furthermore, material properties which are determinants for the initiation and propagation of failure, such as E 33 , G 13 , v 13 , G and Z T are not always available or readily measurable. Further research study to validate the model in excessive damage propagation such as compact tension and dynamic conditions (e.g. low velocity impact) is necessary. The current constitutive model, would require further adjustments to accurately model other types of NCF composites, such as triaxial or quadraxial where the interaction between the layers is more complicated than of a cross-ply.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A.3. Recorded values of in-plane fibre misalignment (in-plane shear)
See Table A.12.

A.4. Recorded values of out-of-plane fibre misalignment (in-plane shear)
See Table A.13.

A.6. Recorded values of out-of-plane fibre misalignment (compression)
See Tables A.20-A.26.   Table A. 26 The average values of in-and out-of-plane fibre misalignment of the compression specimens for the different nominal fibre disorientation.