A Plane-Stress Damage Model for Vectran Laminated Composite

The aim of this paper is to present a plane-stress damage model based on the Classical Lamination Theory (CLT), developed for polymer fibre-based composite. The proposed numerical model utilises a damage mechanics methodology coupled with fracture mechanics to predict composite failure, particularly under quasi-static and dynamic loadings. In addition, the proposed constitutive equations consider a single secant modulus to describe its tensile and compressive modulus, as opposed to the physically-proposed tier models for polymer fibres which possesses a ‘skin–core’ structure. The result of single element and coupon-level modelling showed excellent correlation with the experimental results. It is expected that the proposed numerical model will be able to predict, up to a considerable accuracy, the response of the composite under low and high velocity impact loadings.


Introduction
The problem of impact on composite structures has been a subject of review for more than three decades. To date, many review papers have been written by researchers [1][2][3][4][5][6][7][8][9] in reporting the advances observed in the field of impact mechanics on composite materials. These advances, which include damage prediction using numerical approaches such as the Finite Element Method (FEM), have strengthened our understanding to design a more damage tolerant structure for various applications. Mathematical models such as the spring-mass model, the energy-balance model and the Delamination Threshold Load (DTL) [10,11] have greatly helped in determining the performance of the composite structure under transverse impact loading, particularly under Low Velocity Impact (LVI).
The use of energy-based models has also seen an appreciable increase over the past three decades. Mathematical models such as the Classical Laminate Theory (CLT) are used in conjunction with energy-based theories such as fracture mechanics to describe laminate failure. For example, the energy required for fibre failure is taken as the energy obtained from standard fracture mechanics tests such as the Compact Tension (CT) or Compact Compression (CC). This energy is commonly known as the strain energy release rate, G C , and is assumed to be the area enclosed under the stress vs strain curve of the relevant modes. Upon reaching the failure strength (tension or compression), laminate stiffness is degraded to zero. The concept of degradation is essentially part of a more general Continuum Damage Mechanics (CDM) approach, which was first introduced by Kachanov [12] and later by Rabotnov [13] when attempting to describe the creep behaviour in metals. CDM is an attractive approach since it provides a method in which accurate determination of the material condition can be made -from a pristine condition (no damage) until final failure (full damage). Unlike earlier approaches to modelling composite laminate using stress (or strain) based criteria [14][15][16], in which the stress is immediately reduced to zero, upon reaching its threshold strength. This is a gross over-simplification which neglects the post-failure behaviour of a laminated composite.
The earliest implementation of energy-based damage mechanics approach is proposed by Ladeveze [17]. In-plane testing of various laminate orientations was simulated using the CDM approach and later compared with the experimental results. Damage onset was taken as the failure strength of relevant modes, and then linearly degraded until zero. An excellent correlation between the experimental and simulation results was obtained. Later, Matzenmiller et al. [18] utilised the CDM approach by considering the post-failure behaviour as a function of its Weibull distribution of strength. Williams et al. [19] implemented the approach suggested by Matzenmiller et al. [18] into LS-Dyna as a plane stress material model. Following this, Iannucci et al. [20][21][22][23][24] employed the CDM approach to model thin laminated composites (UD and woven) under LVI and HVI loadings. All numerical predictions including the load vs time/displacement response were in close agreement with the experimental results. Pinho et al. [25] and Donadon et al. [26] utilised the model first developed by Ladeveze [17] and Iannucci et al. [20] to develop a full 3D material model to enhance the accuracy and reliability of the model. Donadon et al. [26,27] have also implemented the smeared cracking approach to alleviate mesh dependency into irregular/ non-structured mesh types. The approach utilises some modifications on the element shape functions to accommodate non-structured mesh types used in FE analysis. Ehrich [28] has also implemented this modification into plane stress elements for the Abaqus FEM software package.
It is well established that Vectran fibres (commonly known as a type of Thermotropic Liquid Crystal Polymers, TLCP) possess the so-called 'skin-core' structure [3] [4], which determines the mechanical performance of the fibre. Differences in the processing conditions (i.e. draw ratio, temperature, etc.) will result in a change in the fibre structure, often manifesting into a different 'skin-core' ratio (i.e. the thickness of the 'skin' would increase or decrease), hence altering the mechanical properties of the fibre. For instance, a thicker 'skin' would result in an increase in the fibre tensile modulus, and a decrease in the tensile strength and failure strain. In general, the 'skin' structure is often identified by a smooth surface, often with a highly oriented crystalline domain. Moving in towards the 'core', the molecular structure was found to be less oriented, physically seen as fibrils in the fibre.
In this paper, a plane stress damage model established for Vectran laminated composites will be presented. The formulation follows previous works from Iannucci et al. [21][22][23], with some modifications derived specifically for the tensile and non-linear in-plane shear behaviour observed on Vectran composites [29].

Constitutive Model
The damage model proposed in this paper is based on a plane stress formulation suitable for thin and thick shell element type. The assumptions made for the model include: • The maximum strain failure criterion was used as damage onset and growth for both tension and compression. The stiffness variable used in both tension and compression response is based on the 'skin' and 'core' structure found on the Vectran fibre, (Fig. 1) [30][31][32]. Although it is ideal to define the two (or more) stiffness variables explicitly for each structure (i.e. skin, core), it was found that the secant modulus approach was adequate to represent the material's behaviour; • Explicit damage variables were assigned to each failure modes (i.e. d 1 for tensile failure in the warp direction, d 5 for compressive failure in the weft direction, etc.); • Linear softening was assumed upon damage propagation for all damage modes. While the selection may seem arbitrary (i.e. no physical reasoning), the linear softening approach has been widely used by many researchers due to its simplistic nature. All cracks and damage types (i.e. fibre bridging and pull-out) were smeared into a. Representative Volume Element (RVE), analogous to the assumptions made by CDM; • Behaviour in the through-thickness direction (both normal and shear stress) was defined in a 'penalty' (linear) form in this UMAT. Through-thickness shear modulus was assumed to be similar to the in-plane shear modulus, whilst the through-thickness modulus was taken as the resin modulus (calculated using the rule of mixtures).

Governing Equations
From the initial assumptions made from CDM, the relationship between a damaged and undamaged material can be defined as: Where is the effective stress in the 'damaged' material, and σ is the pristine 'undamaged' material. d=1 represents full damage (i.e. = 0), and d = 0 represents no damage (i.e. = σ). The relationship defined in Eq. 1 can also be extended to define the material's stiffness, which is defined as: where E is the effective (current) stiffness, and E 0 is the pristine material's stiffness. From Eq. 2, the damage modes can be defined as: • d 1 : damage due to tensile stresses in the warp direction; • d 2 : damage due to tensile stresses in the weft direction; • d 3 : damage due to in-plane shear stresses; • d 4 : damage due to compressive stresses in the warp direction; 5 : damage due to compressive stresses in the weft direction.
All damage modes defined were implemented on the local element direction following Eq. 2 which is defined as: where t represents tensile failure, and c represents compressive failure. For tensile failure, Eq. 3 can be defined as: where i = 1, 2. For compressive failures: where j = 4, 5. Since the same stiffness was used for both tensile and compressive response, Eqs. 4 and 5 can be combined, yielding: where (1 − dx) and (1 − dy) are defined as: For in-plane shear, the stiffness-damage relationship is defined as The Poisson's ratio will also be degraded following Eq. 2 to prevent spurious energy generation. This is defined as: where the relationship between Poisson's ratio and the Young's modulus are defined as: For Vectran/epoxy composite, the general stress-strain relationship can be defined as: and where E 0 is the secant modulus, N and M (1 < N, M < 2 for an optimum curve fit) are constants which can be found by performing a least square fit of the tensile stress-strain curve.
As mentioned earlier, a secant modulus was used to represent the 'skin' and the 'core' structure, which was established for Vectran fibres. Finally, the 2-D stress-strain-damage relationship for each laminae can be defined as: where N is defined as: can also be written in an incremental form, which has been implemented in LS-Dyna, given by: where ∆E 0 is defined as: Therefore, Eq. 15 can also be re-written as: where the β term defines irreversible strain, ε p , from any unloading point during damage propagation. β = 1 will result in an unloading straight back to the origin, whilst β > 1 results in a positive irreversible strain, Fig. 2. The β term could also be used to define the amount of irreversible strain on the laminate -a common damage type usually observed on polymer fibre composites.

Tensile and Compressive Direct Stresses
As mentioned earlier, damage propagation follows a linear softening curve, in which the stress is degraded as a function of strain. For simplification, Poisson's coupling is neglected prior to damage initiation. For both tensile and compressive component, damage propagation is defined as: where A is a constant given by: where ε max is the strain at zero stress (full damage), ε 0 N is the initiation strain, and ε is the current strain. Eq. 18 can also be written in an incremental form, which has been implemented in LS-Dyna, given by: In order to compute the strain at final failure, ε max , fracture mechanics is introduced whereby it is postulated that the energy dissipated during the damage process is equal to the strain energy release rate, G C , of the particular damage mode. For instance, the energy consumed during longitudinal (fibre) tensile failure is equal to that of the energy obtained from fracture mechanics testing, such as the Compact Tension (CT) test. The strain energy consumed during material failure is a volumetric energy and it is generally known as the area under the stressstrain curve (shaded area shown in Fig. 3). Since fracture mechanics considers the energy per unit area, kJ/m 2 , another length parameter has to be introduced, defined as: where l c is the characteristic element length describing the finite element mesh. Therefore, the size of l c can vary with the element size. This inherently alleviates the mesh dependency of the finite element model since the energy under the curve must equal to the G C to prevent any spurious energy dissipation.
Equation 21 is also the basis of the crack smearing methodology which will be discussed in the later sub-section. Finally, the 1-D stress-strain curve for tensile/compressive component is shown in Fig. 3.

Residual Strength for Compressive Failure
It should be noted that the compressive strength of Vectran laminates is considerably lower compared to its tensile strength. This is due to the fibrillar nature of the fibre, which defibrillate almost immediately when placed under a compressive load [30,31,33]. It was also observed that when Vectran/epoxy laminates were placed under uniaxial compressive loads, the laminate will exhibit a strain-hardening response, which allows the laminate to carry a further load with increasing strain. This is shown in Fig. 4, where it can be seen that non-linearity begins at approximately 40 MPa (ε = 0.021%). Upon passing its yield strength, the fibres will start to buckle, seen as micro-kinks at the fibre level [33]. Strainsoftening begins until approximately 80 MPa (ε = 4.3%) when finally, the specimen was observed to be 'crushing'.
It is common for brittle fibre -brittle matrix system such as CF/Epoxy and GF/Epoxy laminates that the residual compressive strength is estimated from the matrix compressive strength, which is considerably lower than the fibre compressive strength. However, for Vectran/Epoxy laminates, the residual compressive strength is taken as the ultimate strength before 'crushing' (i.e. approximately 80 MPa in this case). Therefore, upon reaching the maximum compressive strength of the laminate, no degradation of the compressive strength is being made, emulating an elastic-perfectly plastic response. Again, although this assumption could be seen as an oversimplification of the experimentally observed behaviour, it is assumed that the compressive response is significantly lower than its tensile component. Furthermore, if one were to fully capture the actual response of the laminate (i.e. kink band effects etc.), 3D-stresses are required which would necessitate a full 3D formulation. Several researchers have proposed a method to predict the kink-band formation for polymer fibres, especially the ones which possesses 'skin-core' structure [34,35]. The formulation was derived based on the assumption that the 'skin-core' structure was made of layers, with each layer connected with a spring (through-the-thickness), as well as frictional effects between each layer. Also, a separate spring is assigned to govern the kink-band formation when the 'damaged' layers start to rotate.

In-plane Shear
A simple approach in modelling the in-plane shear response, including the non-linear component normally observed in many in-plane shear test was implemented in LS-Dyna. The model was based on a 'piece-wise' approach with an exponential function to predict the nonlinear response. A linear softening follows upon reaching the ultimate shear strength until zero stress. For shear strain below the yield strain (strain at which non-linearity commences), a linear shear stress-shear strain relationship is proposed, given by: where G 12 0 is the initial shear modulus obtained from standard laboratory testing, and α relates to the strain-rate enhancement case defined by a sigmoid relation initially proposed by Iannucci et al. [21][22][23] given by: where κ is a material parameter defining the shear strain-rate, γ, with respect to the magnitude of the enhancement. Fig. 5 presents the predicted shear response under various strainrate. Enhancement can be observed in terms of the initial shear modulus, G 12 0 , the initial yield strength, τ 12 0 , and the ultimate with respect to the strain-rate. Also shown in Fig. 5 is the sigmoid relation, where it can be seen that the enhancement will always asymptote to a factor of 2 after reaching a strain-rate of 100/s and above.
Following this, the shear stress-shear strain response follows an exponential function, given by: where a, b, c, and d can be found be performing a least square fit on the non-linear component of the shear stress-shear strain curve. Finally, upon reaching the ultimate shear strength (damage onset), a linear degradation is proposed, given by: where τ 12 u is the ultimate shear strength, γ max is the maximum shear strain (damage = 1), and γ 12 u is the ultimate shear strain. The maximum shear strain, γ max can be obtained by using fracture mechanics, briefly discussed in the previous section. This can be estimated from the mode II interlaminar fracture test [36], in the absence of mode II intralaminar energy, shown by the area under the shaded curve in Fig. 5. Writing in terms of damage, Eq. 26 is given by: where the constant C is given by: Finally, converting Eq. 26 into an incremental form yields: Using the model defined above, the predicted shear stress-shear strain response is shown in Fig. 5.

Unloading Behaviour
Experimental evidence of the in-plane shear unloading response for Vectran/Epoxy shows an almost similar unloading stiffness with the initial modulus. This is shown in Fig. 6, where the unloading stiffness is nearly constant with increasing shear strain.
The total strain can be additively decomposed into the elastic and inelastic strain, given by: [37] where γ 12 e is the elastic (recoverable) strain, and γ 12 i is the inelastic (permanent) strain. This is schematically shown in Fig. 5, where unloading from point A will result in a residual permanent shear strain OB. The magnitude of permanent shear strain will increase with increasing strain -unloading from point C yields a permanent strain with magnitude OE. Upon passing the ultimate shear strain, γ 12 , stress will be degraded with a negative stiffness, defined in the previous section. Any unloading beyond this point will follow the damaged stiffness, defined as: It must be noted that unloading from any point before damage initiation will utilise the initial shear stiffness, G 12 0 , consistent with the results obtained experimentally.

Mesh Dependency
The proposed formulation utilises the established "smearing methodology" first proposed by Bazant and Oh [37], in which the fracture energy, G c is assumed to be smeared over a Representative Volume Element (RVE). This approach inherently alleviates the mesh dependency due to strain localization by coupling fracture mechanics and damage (26) mechanics in the constitutive equations. Upon failure, energy dissipated within an RVE can be defined as: where E e is the elastic strain energy, and E d is the 'damaged' strain energy. Eq. 31 can also be re-written in similar form following Eq. 21, given as: where V is the volume of the finite element, and Y is the thermodynamic force conjugated to the damage variable. This can be further deduced to define the characteristic element length, l c , which is an intrinsic requirement of the smearing methodology. Hence, Eq. 32 can be rewritten as: Figure 8 illustrates the smearing methodology employed in this UMAT. Square specimens of 1 mm 2 were discretised into 1 × 1, 3 × 3, 5 × 5, and 7 × 7 element sizes respectively -shown in Fig. 9. As expected, failure occurs in a band normal to the loading direction, without an outside trigger to initiate failure. Consequently, the load vs displacement curves shown in Fig. 7 clearly indicates a similar energy dissipation (area under the curve) for several different mesh sizes. The energy dissipation is a function of the element volume (i.e. strain energy density), therefore ensuring a constant and stable damage propagation.

In-plane Tests
The proposed plane stress damage model described in the previous sub-sections has been implemented into LS-Dyna for thin and thick shell formulations. The next sub-sections will present results of the in-plane modelling of standard coupon specimens for tension, compression, and ±45° in-plane shear. Finally, results from the FE modelling will be compared with the actual in-plane experiments therefore enabling an appropriate validation of the UMAT model.

Model Set Up
The model properties used in the FE model were summarised in Table 1, whilst the constants used were summarised in Table 2. Engineering properties (both in-plane and fracture toughness) were obtained by performing quasi-static mechanical characterization on the actual NCF laminate, shown in Fig. 10.
It must be noted that the Polyester stitches which are present in the fabric/laminate consists of approximately 7% of the total laminate volume fraction. Also shown in Fig. 10 is the Representative Unit Cell (RUC) of the fabric. Finally, one end of the specimen was constrained against translation along and across the specimen, whilst the other end was subjected to a ramped velocity boundary condition. The ratio between internal and kinetic energy were kept below 0.1% to ensure quasi-static conditions.

Tension
A tensile FE model was created with specimen dimensions defined as 250 mm x 25 mm x 2 mm (L x W x t) utilising quadrilateral thin shell elements. The total elements were 6,250, with 6 through-the-thickness integration points in each element, corresponding to 6 Vectran NCF plies. 5 damage variables were assigned to each integration points, corresponding to the damage variables described in the previous sub-sections. Figure 10 shows the FE prediction on the tensile stress-strain curve for Vectran/MTM57. Also shown in Fig. 9 is the experimentally obtained stress-strain curve from the quasistatic tensile testing of Vectran/MTM57. Excellent agreement can be seen on the experimental and FE predictions with similar ultimate stress and strain values obtained. Conversely, Fig. 10 shows the predicted failure in the FE model. As expected, failure occurs close to the constrained end of the specimen, where the stress concentration was highest. Damage initiates at the side of the specimen, progressing to the centre of the specimen, which finally failed when all integration points in the damage band have failed. Figure 11 Failure on the actual specimen was seen further away from the end tabs. Fig. 12b, and not very close on the end tabs, as predicted by the FE model. This can be explained from  the inherent defects which were present around the failed area, which initiate damage which finally leads to failure. FE model however does not consider any defects, and to account for this, more advanced theories such as the stochastic approach should be employed.

Compression
For compression, an FE model with specimen dimensions defined as 90 mm x 10 mm x 4 mm (L x W x t) was created using 900 quadrilateral thin shell elements. 12 throughthe-thickness integration points were assigned in each element, corresponding to 12 Vectran NCF plies. Similar to the tensile modelling, 5 damage variables were assigned to each integration points enabling a comprehensive investigation of damage in each point. As expected, failure occurs at the gauge region of the model, similar to the experimental observation -shown in Fig. 12a and b. In spite of this, the kink-bands observed in the actual experimental specimen were not be able to be simulated, since this requires a through-the-thickness formulation to be included in the UMAT. Figure 14 presents the stress vs strain response for Vectran/MTM57 under quasi-static compression testing. As mentioned before, the residual strength for Vectran/MTM57 compressive response were set at approximately 80 MPa, similar to the experimental observation. Although the FE model does not consider the strain-hardening response of Vectran/MTM57, the overall response, including impact was found to be able to predict the experimental response well. One possible approach to model fibre folding and kinkband formation is by using the stacked shell approach with each layer 'connected' using cohesive algorithms (cohesive surfaces/elements). One integration point will be assigned to each layer, resulting in a multiple 'membrane' ply stacked onto each other. Upon reaching the yield stress, the plies will start to deform forming a 'kink-band' similar to the experimental result. However, this approach is computational costly since each layer is required to be modelled, along with appropriate mesh sizes.
Numerically, the compressive strength could also be estimated analytically using the relationship first proposed by Fleck [38] assuming an initial misaligned angle in the composite given by: where R is the ratio of matrix tensile to shear strength, G is the matrix shear modulus, γ Y is the matrix shear yield strain, β is the kink-band angle, and is the initial misalignment angle. Taking G = 0.116 MPa, R = 1.335, and = 2°, and β = 0° yields σ c = 44.85 MPa. This is in excellent agreement with the experimentally obtained compressive yield strength

±45° In-plane Shear
Finally, the quasi-static ±45° in-plane shear test was modelled using 6,250 quadrilateral thin shell elements utilising similar dimensions with the actual specimen (250 mm x 25 mm x 2 mm). Six through-the-thickness integration points were assigned in each element, corresponding to 6 Vectran NCF plies. Finally, the orientation of the material axes in each element was rotated ±45° with respect to the loading direction -shown in Figure 15a so as to represent the rotated plies in the actual specimen  Figure 15 presents a comparison between the experimental and FE specimens. Consequently, Figure 16 shows the FE prediction for the ±45° in-plane shear stress vs strain behaviour, in which an excellent agreement can be observed between the experimental response. At first glance, a similar 'necking' behaviour can be observed between the FE model and the actual specimen. However, the 'necking' observed around the gauge region in the actual specimen was due to fibre rotation (i.e. scissoring effect), upon passing its 'yield' strength. In the FE model, 'necking' was due to the element co-rotational formulation (Belytschko-Lin-Tsay shell [39]), since the element local direction (fibre direction) was kept at a constant 90° relative to each other throughout the analysis. It was verified on a single element as well as coupon tests that the longitudinal and transverse strain in the ±45° in-plane shear were in good agreement with the experimental results for shear strains of up to 20%.

Tensile Behaviour
Behaviour under tensile loading showed excellent agreement with the experimental results, whereby both stress-strain response and damage pattern were closely predicted using minimal inputs. This is important since it allows users to use the UMAT with minimal experimental information which can be obtained from standard quasi-static in-plane tests. In spite of this, the physically-based two (or three) tier model could be implemented to enable a

Compressive Behaviour
Due to the plane stress formulation implemented in the UMAT, it is inherently assumed that the through-thickness stresses are a constant which was defined using a 'penalty' approach utilizing a linear spring for both through-thickness direct (σ 3 ) and shear (τ 13 , τ 23 ) stresses. This place severe restriction in capturing the through-thickness effects of the laminate. In the case of compressive failure, the distinct kink-bands observed in the throughthickness direction, Figure 13c, cannot be predicted using the current UMAT framework since it requires additional constitutive formulation to describe the particular response. This could be represented using the approach proposed by Wadee et al. [5] [6], in which the kink-band phenomenon was described using explicit constitutive relations, including band width growth following kink-band growth saturation.

±45° In-plane Shear
Under in-plane shear loading, extensive 'necking' is observed on Vectran/MTM57 due to the extensive fibre-rotation phenomenon which occurs progressively with increasing load. This effect was not explicitly captured since it requires a specific formulation which links the rotating angle with the shear stress/strain. As mentioned before, the necking observed in the FE prediction was due to the element formulation, which considers co-rotational formulation at the local (element) level. In spite of this, excellent agreement can be observed between the experimental and numerical results, although the inputs required for the inplane shear is more involved (at least four additional constants are required to obtain the optimal correlation).

Critical Strain Energy
The UMAT framework outlined in this paper inherently requires additional information such as the strain energy release rate, G c , as input due to the coupling between CDM and fracture mechanics to describe failure. While G c for tensile (mode I) has been investigated [29], it is assumed that the G c is similar in compression. This is due to the lack of experimental data in characterising the fracture toughness in compressive failure mode for polymer fibre-based composites. Perhaps the closest work on characterising the compressive behaviour of polymer fibre composites was performed by Attwood et al. [41], on Dyneema HB26 and HB50 composites using a Compact Compression (CC) specimen. It was concluded that compressive failure of Dyneema composites was governed by its interlaminar shear properties and not fibre properties. This is contrary to Vectran fibres, or more generally polymer-based fibres possessing a 'skin-core' structure in which the composite compressive strength is governed by the fibre strength [11,40]. It must be also be noted that in this UMAT model, failure is governed by tension rather than compression or shear. Failure in each integration point occurs when full damage (d = 1) in the fibre tensile direction (d 1 , d 2 ) is reached. Consequently, the element in the FE model is deleted when full damage in either fibre direction (d 1 , d 2 ) has been reached -max(d 1 , d 2 ), since this option is seen to yield a much better correlation with the experimental results.

Conclusions and Perspectives
The FE model proposed for Vectran/MTM57 was found to be able to predict the response of Vectran/MTM57 with considerable accuracy despite several drawbacks with are inherent under plane stress assumptions. However, the UMAT developed is generally acceptable for thin laminates due to the flexural dominance in terms of impact deformation. At the coupon level (i.e. quasi-static characterisation tests), FE prediction was observed to be in close agreement with the experimental results, apart from compressive behaviour, in which the 'hardening' effect should be addressed using kink-band models which requires a full 3D formulation. While the FE prediction on the behaviour of Vectran/MTM57 correlates well with the experimental observation, the mathematical framework for Vectran fibre composite (or any fibres which possess a 'skin-core' structure) could be further improved to accurately represent, based on physical observations, the fibre's behaviour. First, a 3-D mathematical framework is necessary to predict the through-thickness permanent indentation often found on Vectran/Epoxy laminates. This could be incorporated in -conjunction with the modified Hertzian law [42] or coupled with CDM [43] hence would be able to accurately predict the performance of the composite. Secondly, the primary governing equation of Vectran's fibre could be further enhanced by using two (or more) non-linear springs, as well as a dashpot to represent the skin and the core. This is then consistent with the hierarchal model originally proposed for fibres possessing a 'skin-core' structure [31]. Thirdly, the in-plane compressive behaviour can be represented by the kink-band model proposed by Wadee et al. [44,45], therefore having the ability to predict the onset and growth of the 'kink-band' upon reaching the fibre's compressive yield stress. Lastly, enhancements can be made to modify the shape functions for the elements to accommodate the use of non-structured mesh (i.e. warped quadrilateral mesh) which can be done by including an objectivity algorithm as performed by Donadon et al. [27] and Ehrich [28] for Abaqus FE software package.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.