Numerical validation of an improved model for the shearing and transverse crushing
of orthotropic composites
S´ergio Costaa,b,∗, Thomas Brua,b, Andr´e Portugalc, Robin Olssona
aSwerea SICOMP AB, Box 104, SE-431 22 M¨olndal, Sweden
bChalmers University of Technology, Department of Industrial and Materials Science, SE-412 96 Gothenburg, Sweden cINOV Contacto Programme
Abstract
This paper details a complete crush model for composite materials with focus on shear dominated crushing under a 3D stress state. The damage evolution laws and final failure strain conditions are based on data extracted from shear experiments. The main advantages of the current model are: no need to measure the fracture toughness in shear and transverse compression, mesh objectivity without the need for a regular mesh and finite element characteristic length, a pressure dependency of the shear response, account for load reversal and for some orthotropic effects (making the model suitable for Non-Crimp Fabric composites). The model is validated against a range of relevant experiments, namely a through-the-thickness compression specimen and a flat crush coupon with the fibres oriented at 45 and 90 degrees to the load. Damage growth mechanisms, orientation of the fracture plane, nonlinear evolution of Poisson’s ratio and energy absorption are accurately predicted.
Keywords: Crushing, Damage mechanics, Friction, Finite element analysis (FEA)
1. Introduction
The structural weight of vehicle body structures can be greatly reduced by implementing composite materials. Weight reduction has been found to be the most efficient non-engine related emission-reducing technology for auto-motive [1]. Besides their excellent specific stiffness and strength, composite materials can offer better specific en-ergy absorption than metals during crash events [2]. Crash events are structural phenomena while the term ”crush” here is used to describe continuous compressive loading of the material beyond the initiation of damage and peak stress. When applied to composite materials it is com-monly used to define the complex mechanisms of energy absorption taking place at different scales in the mate-rial [3]. The crush mechanisms are not fully understood, let alone included in a predictive model. The main factors responsible for the difficulties associated with crash mod-elling of composites have been identified in [4]: (i) the large strains and rotations involved, (ii) interaction of complex damage mechanisms, (iii) significant geometric transfor-mations in the crushing zone, (iv) material nonlinearity, (v) contact and friction and (vi) potential strain rate ef-fects. To ensure a predictive evaluation the energy ab-sorbed by a composite structure undergoing crushing, all the aforementioned aspects of composite crushing need to be accounted for in a material model without calibration
∗Corresponding author
Email address: sergio.costa@swerea.se (S´ergio Costa)
from experimental crush test data and with a step-wise validation. Starting with cases involving lower amount of mechanisms and low interaction between them and end-ing with more challengend-ing situations. In the current work, some of the tests used for material characterization were conveniently used to validate some of the assumptions used in the model. This could constitute a benchmark case for further crash model developments.
The continuum damage mechanics (CDM) approach has been generally accepted as a practical methodology for crash modelling [5, 6, 7]. A large number of CDM models for fibre-reinforced polymers (FRPs) are based on the work by Ladev`eze and Le Dantec [8]. The damage is handled at the macroscopic level using damage variables, which are associated with a particular failure mode and are used to degrade the material stiffness properties associated with that failure mode. The evolution of the Poisson’s ra-tio with damage has received low attenra-tion lately and it is still an unresolved issue although an important one [8]. The damage variables characterise microscopic cracks and voids that accumulate in the material before macroscopic failure occurs. Considering the case of the Matrix Com-pression (MC) damage mode, the damage mechanism up to final failure of an FRP ply is the nucleation and growth of microcracks in the polymer matrix until they eventu-ally coalesce to a macroscopic crack [9]. A similar failure sequence is also observed for the shear failure of a ply [10]. Friction occurs at the contact between microcrack surfaces and should be considered as a source of energy dissipation. Gutkin and Pinho [11] were the first researchers to couple
damage with friction in a progressive damage model. It was shown that internal friction effects could account for a part of the nonlinear response and the permanent strains observed in shear.
Most material models for composite crash applications available in commercial finite element (FE) software re-quire the use of tweaking parameters that cannot be mea-sured experimentally and need to be calibrated by trial and error (e.g. MAT54 in LS-Dyna [12]). Although exper-imental crash results can successfully be reproduced with MAT54, the model cannot be used as a fully predictive tool. The use of tweaking parameters, which may have no physical significance, also make the prediction with MAT54 extremely mesh sensitive, as observed by Deleo et al. [13]. Another source of mesh sensitivity of FE crash simulations is element deletion, which results in unrealis-tic peaks in the load–displacement curve. Post-processing filters are typically used to eliminate these peaks and ex-tract the crash response [14]. However, the use of filters should be avoided since it may affect the evaluation of the energy absorption and the value of initial peak in the load–displacement curve, two important parameters to as-sess the crashworthiness of a structure.
During FE simulations of composites undergoing crush-ing, elements adjacent to those which have failed are likely to experience unloading/reloading and/or load reversal. It is therefore important to account for the loading history during the simulation and to allow for the interaction be-tween the different failure modes [5]. Chiu et al. [5] ad-dressed load reversal for longitudinal failure modes (i.e. in the fibre direction) modelled without permanent strains but, to the knowledge of the authors, there is no reference to addressing load reversal for transverse behaviour where permanent strains are developed in compression. More re-cently, Shi and Xiao [15] demonstrated the importance of permanent strains and residual stiffness in energy absorp-tion predicabsorp-tions.
The present friction/damage model with focus on the shearing and transverse crushing is based on the work of Gutkin and Pinho [11] with the addition of major features to allow for orthotropy, interaction of damage modes and mesh objectivity. Additionally, four intralaminar damage modes are combined in a physically based way and the model response is verified in several situations. Finally the crush response is validated with two crush coupons. The MC damage mode accounts for pressure dependency and distinguishes between in-plane and out-of plane com-pressive response as well as in-plane shear and out-of-plane shear. Furthermore, intralaminar damage modes such as Matrix Tension (MT), Fibre Compression (FC) and Fi-bre Tension (FT) are also implemented and interactions between MC and MT are introduced to account for the load redistributions typically observed in the process of crushing composites. For MC, an experimentally based condition for final failure strain is introduced to repre-sent fully-developed macroscopic cracks. Using this fail-ure condition, it is possible to both avoid toughness
mea-surements and to differentiate between in-plane transverse and through-thickness behaviour. Therefore, the model can be used for plies that are not transversely isotropic, e.g. plies reinforced with unidirectional (UD) Non-Crimp Fabric (NCF). The interaction of matrix damage in the fibre compressive response is also included, although fibre kinking is simplified in the present model formulation as FC damage mode. In the future, the recent work of Gutkin et al. [16] and Costa et al. [17] on matrix-driven kinking failure will be included in the model.
2. Description of the model
2.1. Overview of the model
For simplicity the model is separated into four dam-age modes: Matrix Tension (MT), Matrix shear and Com-pression (MC), Fibre ComCom-pression (FC) and Fibre Tension (FT). The damage modes and their respective interactions are illustrated in Fig. 1(a). MC and MT are based on the fixed crack concept, i.e. the material constitutive law is de-fined on a potential fracture plane which is fixed when the failure index (fi ) reaches unity and remains fixed during damage growth, Fig. 1(b). MC is the result of combined shear forces and pressure acting on a fracture plane, so by definition it also includes the pure shear damage modes (longitudinal and transverse shear). MT and FT are based on the work of Pinho et al. [18] and a FC is modelled by maximum stress criteria with toughness-based degra-dation. A simple formulation was used for FC since the focus of this paper is on transverse failure and not on fibre kinking.
The calibration procedure for MC is performed us-ing experimental results of cyclic Iosipescu shear tests re-ported in [19]. The results of the calibration are shown in Fig. 2(a). In order to determine the damage evolution law of the material in shear, the average shear stress and shear strain in the specimen gauge section are monitored in cyclic tests with increasing stress levels. The nonlinear shear response of the material model is obtained by com-bining damage and friction (that occurs at the contact be-tween microcrack surfaces). The inelastic behaviour and the hysteresis loops are achieved using a stick/slip be-haviour rather than by plasticity. By doing so, it is possi-ble to account for the most important characteristics of the shear response in a physically based way. Virtual Iosipescu specimens were created in ABAQUS/Explicit 2016 to ver-ify the implementation of the MC damage mode. Fig. 2(b) shows that the predicted strain fields agree fairly well with the strain fields measured experimentally.
2.2. Constitutive relations 2.2.1. Matrix damage
Ply damage by matrix compression/shear (MC) is pre-sented in this section. The case of matrix tensile failure (MT) is described later on.
(b)
2 1 3 Nσ
Tτ
Lτ
α
(a)
3 2 3 2 1 2 1 2MT
FC
MC
FT
3 2 1 2Figure 1: (a) Intralaminar damage modes and respective interactions with double-headed arrows (red); (b) Coordinate
system (123) rotated to the fracture plane system (LNT).
In the current model friction is not included in the fail-ure criterion, but combined with damage in a physical way. Thus, we use a failure index for MC (fimc) to determine the onset of damage (i.e. onset of nonlinearity), as follows:
fimc= max τL τ0L 2 + τT τ0T 2 + hσNi+ YT 2 (1)
where YT is the transverse tensile strength, τ0L and τ0T are the stresses associated with the onset of damage in the case of longitudinal and transverse shear, respectively. The operator h·i+ in Eq. 1 is the Mc-Cauley bracket defined by hxi+= max{0, x}, meaning that the transverse term of the traction vanishes for normal compressive stresses. The traction vectors τL, τT and σN are acting on the fracture plane oriented with the angle α shown in Fig. 1(b). The fracture plane is searched for in the interval α ∈ [0, π] and is fixed to α0once the failure criterion (Eq. 1) is fulfilled. The stresses acting on the fracture plane are then degraded
according to σN τL τT = (1 − dmchσNi+/σN) ˜σN (1 − dmc) ˜τL (1 − dmc) ˜τT + 0 dmcτLfriction dmcτTfriction , (2) where dmc is the damage variable for MC, ˜σN, ˜τLand ˜
τT are the linear elastic components of the stress tensor in the fracture plane and τLfriction and τTfriction the shear stresses associated with friction.
The physical representation for coupling friction with the damage is that damage in the material assumes the form of microcracks, which when in contact generate stresses proportional to the amount of damage. Therefore, friction only affects the response after damage initiation (friction does not influence the evaluation of fimcwith Eq. 1). Asso-ciating friction with damage is more physically based than including friction already in the failure index to increase the strength as it is done in well-know criterion for matrix compressive failure, such as Puck and Sch¨urmann [20] and Pinho et al. [18] .
For σN ≤ 0, the friction stresses τLfriction and τ friction T are found using a Coulomb sliding criterion to account for the stick/slip behaviour [16]:
τfriction L τTfriction =±µL(|σN| + p0) ±µT(|σN| + p0) , if sliding τfriction L τTfriction =G12(γL− γ0L) G23(γT − γ0T) , if sticking (3)
where µL and µT are the friction coefficients in longitu-dinal and transverse direction, respectively. An internal pressure p0 is introduced, as suggested by Gutkin and Pinho [11]. γL and γT are the respective longitudinal and transverse components of the strain on the fracture plane. γ0Land γ0T are the strains at the onset of damage and are obtained by dividing τ0Land τ0T by the respective elastic shear modulus G12 and G23.
When unloading occurs the sticking behaviour takes place followed by a softening slip behaviour as described in [16]. This behaviour represents the physics well but does not allow for tensile loads. In the present work we allow tensile loads to occur after damage being introduced in compression. If tensile loading occurs after compres-sion, the stiffness of the slip behaviour is kept in tension. To ensure this happens, the remaining friction when ten-sion starts is kept constant, thus keeping the permanent strains developed during compression (visible in Fig. 6). Therefore, for σN > 0 the friction stress is defined as
τfriction L τTfriction =τ friction L τTfriction old , if σN > 0 (4)
Thus, this condition ensures the continuity of the stress state for σN > 0 and accounts for the permanent strains developed on the fracture plane.
0.0 0.2 0.4 0.6 0.8 1.0 0 20 40 60 80 100 0 2 4 6 8 10 12 dmc( ) τL(MPa) γL(%) Material point (FE) Experimental d το Final failure 9.1% τL=68 MPa FE Experiment CALIBRATION VERIFICATION (a) (b) γ cr 1 2 1.6 0.0 3.2 4.8 6.4 7.6 γL (%) ε22 (%) 0.4 0.2 0.0 -0.4 -0.2 τL≈ 68 MPa 9.1 23
Figure 2: Shear calibration and verification of the model using experimental data for HTS45/LY556 composite.
The damage variable for MC, dmc, represents the devel-opment of matrix microcracks during the damage process in compression and/or shear. The damage evolution law is similar to the one in [16]:
dmc=
|γp| − γp 0,mc
γf,mcp − γ0,mcp , (5) where γ0,mcand γf,mcare the strains at damage initiation and at full decohesion respectively. The exponent p is used to achieve a good agreement with the experimental shear stress–strain curve. The strain to drive the damage variable is given by γ = q γ2 L+ γT2 , if σN < 0 γ = q γ2 L+ γ 2 T+ ε 2 N , if σN ≥ 0 (6)
where εN is the strain normal to the fracture plane.
2.2.2. Final failure condition
The stress–strain shear curve extracted from a shear specimen is used as an input to the model. Once the speci-men breaks completely it is no longer representative of the material point. Although the response beyond the com-plete failure of the specimen is still necessary for crash ap-plications. The constitutive response is rapidly degraded at the same level as the specimen breaks. Since the shear stresses at this point are nearly constant a strain based fi-nal failure criterion is introduced to truncate the response as: f icrmc= s γL γcr L 2 + γT γcr T 2 + hNi+ cr N 2 (7) where γcr
L is the shear strain at failure of Iosipescu specimens measured in [21, 19] and γTcr is obtained from either in-plane or out-of-plane compression tests (already performed for material characterization). crN is obtained
from transverse tension tests and only actuates when in tension. Note that the current properties refer to a UD material, and may have to be adjusted for in-situ effects in a multidirectional layup. Once this failure index reaches the unity the damage variable dmc rapidly increases until one and the shear stresses reach the values τifriction (i = L or T ) given by Eqs. 2 and 3. The numerical handling of this transition is explained in the section “Mesh objectivity of MC”.
The influence of stress combinations is illustrated in Fig. 3 The model exhibits the correct trend for the
influ-0 20 40 60 80 100 0 2 4 6 8 10 τ12(MPa) γ12(%) σ22= 11.6 MPa (40% of Yt) Pure shear (Calibrated) σ22=-52 MPa (40% of Yc)
Figure 3: Model response for in-plane shear with combination of transverse compressive and tensile stresses.
ence of transverse compression and tension in the shear response, i.e. higher shear stresses when combined with transverse compression and lower when combined with transverse tension. It’s difficult to judge if the effects of friction are not overestimated so a full experimental vali-dation of the influence of transverse stresses in the in-plane and out-of-plane shear response is recommended but not performed here. In order to further validate Eq. 7 the whole failure envelope for a combination of shear and com-pressive stresses was plotted and correlated with results in the literature, Fig. 4. The failure envelope was obtained
by first applying a transverse load (compressive and ten-sile) followed by shear. Overall there is a very good match with the data in the literature indicating that using the strain based final failure criterion to truncate the response we are able to capture the final failure envelope well. The tensile part seems to be slightly unpredicted. It is worth to point out that the main ambition is to capture the right energy absorption and being able to predict the correct strength is not a guarantee that the right energy is also predicted, although it is a good indication. Note that all the material non-linearity in the current model is due to damage growth (micro-cracking) during shearing, and that any residual strains are due to stick-slip effects associated with friction on the micro-crack surfaces. Thus, τ0L and τ0T used in this paper are associated with the shear stresses at damage initiation, rather than with the maximum shear stresses commonly used to define ”shear strength”.
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 τ12/SL σ22/Yc Damage initiation Final failure
Final failure, Experiments E-Glass/LY556
AS4/55A 0
Figure 4: Model prediction of combined in-plane shear and transverse compression envelopes of HTS45/LY556 NCF composite
and the strength of E-Glass/LY556 [22] and AS4/55A20 [23].
This final failure criterion is also adapted to account for the orthotropic behaviour of NCF composites. For NCF composites both plane and out-of-plane compression in-volve a shear dominated failure, but the failure stresses and strains are significantly different. For instance the experi-mental values of the in-plane and out-of-plane transverse compressive strengths of a HTS45/LY556 composite are Yc= 129 and Zc= 203 , respectively [21]. To account for the difference in strength between these loading conditions the following two conditions were established. Thus, when compressive load is dominating γcr
T is defined as follows for in–plane dominated compression
|ε22| > |ε33| then γTcr= 2ε22cu (8) and for out-of-plane dominated compression:
|ε33| > |ε22| then γTcr= 2ε33cu (9) where ε22cu and ε33cu are the ultimate normal strains in uniaxial compression tests, and γTcrthe corresponding max-imum shear strain. For the current material ε22cu= 1.71%
and ε33cu= 3.81% as measured by Bru et al. [21]. The pre-dicted damage initiation and final failure envelope for the HTS45/LY556 NCF ply is shown in Fig. 5.
-500 -400 -300 -200 -100 0 -500 -300 -100 Damage initiation Final failure Final failure, Experiments
s
22(MPa)s
33 (M P a ) Yc Zc 0Figure 5: Model prediction of the σ22–σ33envelopes for the damage
initiation and for the final failure of HTS45/LY556 NCF composite.
The shape of the open σ22–σ33failure envelope in Fig. 5 is similar to the ones reported in the literature, e.g. [24]. Thus, orthotropy is included in a simple but effective way.
2.2.3. Transverse tensile damage without previous MC dam-age
For the situations in which matrix compression/shear damage has not occurred the matrix tension behaviour is treated as in [18]. The effect of pre-existing damage from MC on the MT failure mode is addressed in the section: ”Interaction between transverse compression and tension”. The stresses in the fracture plane for MT without previous MC damage are given according to [25] as follows,
σN τL τT = (1 − dmthσNi+/σN) ˜σN (1 − dmt) ˜τL (1 − dmt) ˜τT , (10)
and the damage evolution dmt is given by
dmt= 1 − γ0,mt γ γf,mt− γ γf,mt− γ0,mt . (11)
In Eq. 11 the strain at damage initiation γ0,mt and the strain at full decohesion γf,mt are based on the fracture toughness scaled to the element length.
2.2.4. Longitudinal damage
The longitudinal failure modes (FT and FC), or fibre failure modes, are responsible for a large amount of the energy absorbed in the crushing of laminated composite structures [26]. For the scope of this paper, which is trans-verse crushing, FT and FC have little or no influence, so
simplifications could be made when modelling these failure modes. If the present material model is to be used to sim-ulate the progressive failure of composites in which fibre failure is an important damage mechanism, more detailed models may be required to predict the correct energy ab-sorption (e.g. the kinking model for FC).
The material response in the longitudinal direction is assumed to be bilinear with the strength determined by a maximum stress criterion for both tensile and compres-sive loading and with a corresponding fracture toughness. Damage begins when the failure index reaches unity. The stresses are degraded with the damage variable according to Eq. 12. The damage evolution law for the damage vari-ables df tand df cis driven by the longitudinal strains and the translaminar fracture toughness of the material [25].
σ = (1 − df t) ˜σ Fibre tension σ = (1 − df c) ˜σ Fibre compression
(12)
2.3. Interaction between transverse compression and ten-sion
A common practice in CDM FE models is that a single damage mode is active for a given element, i.e. once a damage mode is activated it is locked for the rest of the simulation. This can be a reasonable assumption when the stress state only changes slightly during the simulation. However, local and rapid load redistributions are expected during crush simulations so it becomes important to allow every material point to experience different damage modes simultaneously, or allow them to occur in sequence.
In the current model the effects of damage in matrix tension do not affect the behaviour in matrix compression, representing thus the mechanism of crack closure in com-pression. However, the damage created in compression will cause a degradation of the properties in tension, as illus-trated in Fig. 6. The failure criterion proposed for MT failure after MC is as follows,
fimc−mt= τL SL(1 − dmc) + τT ST(1 − dmc) + σN Yt(1 − dmc) (13) where SLand ST are the longitudinal and transverse shear strength, respectively. When the material is not damaged in compression, i.e. dmc= 0, Eq. 13 is similar to the crite-rion for matrix tensile failure in [18]. Otherwise, when the material is fully damaged in compression, it cannot carry any tensile stresses. If some damage has occurred from MC it is possible that the resulting failure plane α0 (that verifies Eq. 1) would also be the final tensile failure plane when the load is reversed. The other possibility is that the compressive failure plane α0 is not the critical plane and that the tensile failure plane is oriented perpendicular to the load, as would be the case in pure transverse ten-sion without previous damage in compresten-sion. Therefore, Eq. 13 is only checked for α = α0and α = 0, which results in a good trade-off between computational speed and ac-curacy by ensuring that the two potential fracture planes are examined. -150 -125 -100 -75 -50 -25 0 25 50 -3 -2 -1 0 1 s22(MPa) e22(%)
Matrix tension after comp. Matrix comp. after tension
Figure 6: Loading-unloading curves in the transverse direction.
3. Model validation
3.1. Material properties
The material data of a uni-weave carbon/epoxy NCF composite was used in this work. The data was taken from published experimental results in [21, 19] and are gathered in Table 1 and Table 2. Besides the material elastic constants, the MC damage model requires cyclic shear tests from which the damage and friction parame-ters are extracted. Note that measuring the intralaminar fracture toughness in shear and transverse compression is not needed for the model. The only required toughness val-ues are for FT, FC and MT failure modes and are denoted Gf t
c , Gf cc and Gmtc , respectively. For matrix tension dom-inated cases the mixed-mode intralaminar toughness Gmt c is determined from the experimental inter laminar tough-ness values Gmt
Ic and GmsIIc (Table 2) using the stresses at failure initiation according to Pinho et al. [25],
Gmtc = GmtIc σN0 σ0 mat 2 + + GmsIIc τ0L σmat 2 + GmsIIc τ0T σmat 2 , with σmat = q hσNi2++ τL2+ τT2. (14)
In Table 1 SL and ST are the respective longitudinal and transverse shear strengths, while τ0L and τ0T rep-resents the onset of nonlinearity in shear. µL and µT are the friction coefficients on the corresponding fracture planes. The parameter p0accounts for the apparent inter-nal pressure in the material. The interinter-nal friction coeffi-cients and apparent internal pressure are extracted from experimental cyclic shear tests[19]. The values for µLand µT reported in Table 1 are within the range of the dif-ferent coefficients of friction measured by Sch¨on [28] on
Table 1: Mechanical properties of HTS45/LY556 NCF composite.
Table 1. Mechanical properties of HTS45/LY556 NCF composite.
Elastic properties
Moduli
(GPa)
Poisson’s ratios136
11
E
E
22
9.15
E
33
7.7
12
0
.
28
v
32
0
.
43
4.4
12
G
G233.02aG
13
3.7
13
12Strength properties (MPa) Damage
631
cX
Y
c
129
Z
c
203
X
t
1787
p
0
.
6
o,mc
0
.
6
23
0L
S
L
72
b T T S 34 0
Y
t
Z
t
29
Friction propertiesInternal pressure (MPa) Coefficients of friction
p
0
60
L
0.4
T
0
.
4
Critical strains (%) 1 . 9 L cr
Ncr 0.32 Compression dominated by: 7.62 42 . 3 33 22 cr T cr T
a Obtained using transverse isotropy
b
0 tan 2
c T Y S ,with 062 measured in [21]Table 2: Interlaminar, intralaminar and translaminar fracture toughness of HTS45/LY556 NCF composite.
Interlaminar toughness (kJ/m2) GIc = 0.149 GIIc= GIIIc= 0.69
Stress at interlaminar failure initiation (MPa) σ0
I = 14.7 σII0 = 56.7
B-K parameter for interlaminar failure η = 2.5a Intralaminar toughness (kJ/m2) GmtIc = GIc Translaminar toughness (kJ/m2) Gf t c = 67 Gf cc = 103 a Determined from G
Ic, GIIc and GI/IIc [27]
carbon/epoxy delaminated surfaces. The parameter p is calibrated from the shape of shear stress–strain curve [19]. Xc is the peak stress during compression in the fibre di-rection. In the model, the out-of-plane tensile strength Zt is assumed to be the same as the in-plane tensile strength Ytsince both values represent intralaminar failures. Note that experimentally, laminates under out-of-plane tension fail at a lower stress by an interface (inter laminar) fail-ure [21]. This orthotopy effect is not included in the in-tralaminar damage model, since it can be handled at the interface level in FE simulations. The in-plane and out-of-plane transverse compressive strength (Yc and Zc) are not a required input to the model since the critical shear strains γcrL and γTcr in Eq. 8 and Eq. 9 determine when final failure occurs.
3.2. Mesh objectivity of MC
For increasing shear strains the shear stress and dam-age parameter increase monotonously until the strain reaches an experimentally determined critical strain where the shear microcracks coalesce unstably into a macroscopic crack and the damage parameter rapidly reaches unity. The re-maining low shear stress is caused by the frictional forces on the fracture surface. This very steep strain softening be-haviour may cause instabilities and convergence difficulties when the damage is strain driven [29]. Sudden excessive strains start to localize resulting in oscillations and exces-sive element distortions. Artificial energy will increase and affect the reliability of the simulation. A new method is necessary that does not cause instabilities and still repre-sents the rapid softening.
For this rapid transition phase it was chosen to control the damage parameter by increases in time rather than by strains. This approach avoids the pathological mesh dis-tortions resulting from rapid strain increases during soft-ening, as all elements experience the same time increase although strains may be hugely different. When the final failure criterion in Eq. 7 is reached the damage variable dmc is redefined as:
dmc= 1 − 1 − d0mc (1 + δγ − t/t0) /δγ , (15) where d0
mc and t0 is the damage and the simulation time at the onset of ply rupture, and δγ is the variation of γ on the x–axis. In Eq. 15 the catastrophic failure progresses rapidly in order to have negligible influence on energy ab-sorption. Therefore, d0
mc it is kept within 2% of t0, i.e. δγ = 0.02. A sensitivity study was performed and it was shown that for values of δγ between 0.02 and 0.05 the simulation is stable and that no significant variations are observed in the crash response.
To validate the mesh objectivity a 1 3 cube with one element (1 × 1 × 1) was successively refined to a cube with 5×5×5 elements, Fig. 7. Pure transverse compression load cases along the 2– or 3–axis were applied to all the cubes and tested for mesh objectivity using ABAQUS/Explicit 2016, Fig. 7(a). Off-axis loading was also tested with the main purpose of verifying mesh objectivity for multi-axial stress conditions, Fig. 7(b). The response seems to exhibit the correct trend but was not further validated here.
For the FE results shown in Fig. 7(a) the peak stress corresponds to the point at which the final failure criterion is reached either using the transverse strains from Eq. 8 for out-of-plane loads or from Eq. 9 for in-plane domi-nated loads. The stress plateau is representative of a fully cracked ply in which only the friction acts between the fracture surfaces.
Using the time to drive the damage in situations when the energy absorption is negligible is clearly a robust method, it is able to capture the crack formation (Fig. 8) and re-sults in a mesh objective softening response. The rere-sults for off-axis compression in Fig. 7(b), besides being mesh objective, also look reasonable with higher stresses and lower strains at failure for lower off-axis angles. These results represent an additional interesting feature of the model but their validity is not of key interest in this pa-per.
3.3. FE analysis of a transverse compression specimen There is no test method available to determine directly the pure transverse shear response of composite materials up to final failure. Applying a pure 2–3 shear loading in a V-notched specimen (e.g. a Iosipescu test) usually re-sults in a tensile failure outside of the specimen gauge sec-tion [30]. Instead, the DERA through-thickness compres-sion specimens [31] were chosen to validate the pressure dependency on the out-of-plane shear response. Trans-verse compression failure is indeed a shear-driven phe-nomenon. The uni-axial stresses can be decomposed into
0 50 100 150 200 0 2 4 6 8 10 F/A (MPa) δ22/L andδ33/L (%) 1x1x1 2x2x2 3x3x3 5x5x5 Out-of-plane In-plane 0 100 200 300 400 500 0 1 2 3 4 5 6 7 8 F/A (MPa) δ11/L (%) 1x1x1 5x5x5 10° off-axis 20° off-axis 30° off-axis (a) (b)
Figure 7: FE model response showing the mesh convergence: (a) for the in-plane and out-of-plane compression; (b) for off-axis
compression
shear stresses and a normal compressive stress acting on a potential failure plane.
The through-thickness specimens were milled from a [0◦]187 HTS45/LY556 laminate. Details of the test spec-imens and the experimental setup can be found in [21]. The relatively large gauge section of the specimens allows monitoring the strain at the specimen faces with a Digi-tal Image Correlation (DIC) system. An FE analysis of the DERA specimen was performed in ABAQUS/Explicit 2016. The compressive loading was introduced by apply-ing a vertical displacement δ on a rigid body located on the top of the specimen. Another rigid body located on the base is totally constrained. A coefficient of friction of 0.16 [32] was chosen for the contact law discribing the in-teraction between the specimen and the rigid bodies. For the mesh solid 8-node linear elements with reduced inte-gration (C3D8R) are used, with a size of 0.2 and 0.4 for a finer and coarser mesh respectively. The material properties in Table 1 and Table 2 were used to predict the specimen behaviour.
The compressive stress–strain curves for the tests and the FE simulations are shown in Fig. 8. The strain fields and the fracture plane are also compared. The DIC strain maps show that the shear strains fields γ23 are well cap-tured in the simulation. The model response is elastic only at relatively small strains until ε33 = 1.2%, which corre-sponds to the onset of damage and its associated friction.
As the damage increases, the response becomes nonlinear up to the final failure of the specimen. The nonlinear be-haviour and the final failure predictions are exclusively the result of the parameters extracted from the cyclic Iosipescu shear tests, i.e. no tabular stress–strain was data needed neither for shear nor for compression. Since a reasonably good agreement with the experimental curves is achieved, the model parameters µT and p0in Table 1 can be consid-ered reasonable to account for the right amount of pressure on the failure plane. However, it would be of interest to perform biaxial compression σ22–σ33tests (Fig. 5) to thor-oughly validate the pressure dependency of the damage model.
The final failure criteria based on the experimental crit-ical shear strain predicts the compressive failure occurs 199 (1 in Fig. 8), which is in a good agreement with the experimental strength of Zc = 203 . A matrix crack prop-agates nearly instantaneously in a plane oriented at an angle slightly larger than 45◦ to the load, Fig. 9, as com-monly observed in transverse compression specimens [9] and attributed to friction effects [20]. However in the cur-rent model, and in contrast to Puck and Sch¨urmann [20], friction only exists after onset of damage (label a in Fig. 9). This approach yields a good prediction of the fracture mor-phology of the specimen (2 in Fig. 8).
Modelling the evolution of Poisson’s ratio for damaged materials has received limited attention in the literature since the work of Ladev`eze and Le Dantec [8]. A degra-dation of Poisson’s ratio with damage is not physical in compression, since the coupling between strains increases. In the current modelling the evolution of Poisson’s ratio is left for the model to handle. Fig. 10 shows the evolution of the secant Poisson’s ratio ν32 = −ε22/ε33 during the DERA compression test. The strains were calculated as the averages of the local strains in the entire gauge sec-tion. The Poisson’s ratio is constant for low strain levels and then increases as the compressive behaviour becomes nonlinear. A similar effect was observed for compression of a particle filled polymer inclusion [33]. The FE simulation predicts accurately the increase of ν32 with damage even though Poisson’s ratio in the material model formulation remains pristine.
3.4. FE analysis of crush specimen 3.4.1. Crushing simulation
A final verification of the presented shear dominated crushing model for orthotropic composites is performed by considering crushing of a small NCF coupon. The geome-try of the specimen and its supports in Fig. 11 were taken from Bru et al. [34, 35]. The chosen design is not meant to represent real energy absorbing structures which are gener-ally larger, self-supported (tube or corrugated structures) and with a multidirectional lay-up [2, 26, 5, 7, 6, 36]. The simplicity of the design aims to simplify the identification of the crush mechanisms. Flat specimens have the low-est geometric influence on the crushing behaviour making
them ideal to validate (or identify differences) a material model for crushing.
The composite plates used were [45◦]10HTS45/LY556 fabricated by resin transfer moulding according to the man-ufacturer’s specifications. The experiments were performed under quasi-static conditions using a displacement rate of 0.5 mm/min in an MTS universal testing machine equipped with a 10 load cell. Guiding rods were used to ensure the correct positioning of the crushing plate with the speci-men. The crushing distance δ was evaluated from opti-cal measurement with a laser and the force response was recorded directly from the load cell. Because of the short unsupported region available for crushing, the tests were terminated at δ = 5.6 (i.e. 1.5 after the trigger end). Vir-tual crush specimens were created in ABAQUS/Explicit 2016 with two different mesh sizes, as shown in Fig. 11. Only half of the specimen was modelled to reduce the com-putational time. The input ply properties are reported in Table 1 and Table 2. Rigid plates (analytical surfaces) are placed on the bottom and largest vertical sides of the specimen and fixed to simulate the constraints of the test rig. The friction coefficient between the specimen and the rigid bodies is estimated to 0.16 [32]. The loading was introduced by applying a vertical displacement δ in the crushing plate located at the top of the specimen. For simplicity, a linear elastic material and a bigger element size was used in the bottom of the specimens, i.e. away from the crushing zone.
3.4.2. Experimental and numerical results
The experimental and the predicted crushing response of the specimen are shown in Fig. 12. The crushing mode is mostly by fragmentation of the trigger, which results in a load increase as the width of the trigger increases. Delamination initiation and evolution was modelled as de-scribed by Tan and Falzon [26], using the properties pre-sented in table 2 and the standard cohesive surface con-tact available in Abaqus. Fig. 12 shows that the crushing response is well predicted by the model for the two differ-ent meshes. Using models based on CDM to model crash simulations usually results in sensitivity to element size. A correct crush response with low mesh sensitivity is not only achieved due to the mesh objectivity at the element level, but also as the current model does not rely on the element characteristic lengths. This means that it is able to capture the right energy absorption even for irregularly shaped elements. Delaying element deletion and not using it for model calibration also contributes to a response that is relatively mesh independent. A scalar isotropic damage is activated only when the shear strain reaches a value of 0.8 and subsequently all stress components are smoothly degraded to zero. At this point it is considered that the elements are not carrying any load and can be erased.
The crushing response is clearly shear dominated, as observed on the propagation of the compressive/shear dam-age variable, in Fig. 12. The fibre compressive damdam-age mode is also present and it is able to interact with the
0
50
100
150
200
250
0
1
2
3
4
5
Experiments
FE 0.2 mm
FE 0.4 mm
Compressive strain,
e
33(%)
Co
mpr
ess
iv
e str
ess
,
s
33(M
Pa)
3.810
50
100
150
200
250
0
1
2
3
4
5
Experiments
FE 0.2 mm
FE 0.4 mm
Compressive strain,
e
33(%)
C
o
mp
ressiv
e
str
ess,
s
33(M
Pa)
1 0 damage Just before specimen failure(%)
23
1.5 1.2 0.9 0.6 0.3 0.0 -0.3 -0.6 -0.9 -1.2 -1.5 FE DIC Artefact at specimen surface Just after specimen failure1
2
1
2
Figure 8: Experimental and numerical stress–strain curves of the DERA compressive specimens.
MC. Neither the FC damage mode nor its interaction have significant contribution to the global force response in this particular case. Delaminations between all the interfaces were also included, although neither the quantitative nor the qualitative results were significantly affected. More crash situations that challenge all the aspects of the model will be investigated.
An important metric in crashworthiness assessment is the Specific Energy Absorption (SEA), which can be es-timated from the crush stress of the specimen undergoing crushing, σcr ≈ F/A(δ), and the density of the material ρ = 1509 . The SEA was calculated in the 0.5 ≤ δ ≤ 3 interval for the two FE simulations and is compared to the experimental values of the four specimens tested in Fig. 13. A good agreement is achieved between the ex-periments and the model predictions. Furthermore, the coarse and the fine mesh yield very similar SEA predic-tions. Predicting the correct energy absorption indepen-dently on the mesh size is very important for crash sim-ulations. To the knowledge of the authors, there is no other evidence of mesh convergence crash simulations of composites reported in the literature.
The transverse crushing of a specimen with similar ge-ometry but with fibres oriented transversely with the load was also carried out. This time only the course mesh was evaluated and delaminations were not included. The stress
response correlates very well with the two available exper-iments, Fig. 14.
4. Conclusions
This paper has presented a material model which incor-porates all the intralaminar damage modes necessary for reliable predictions of shearing and transverse crush mech-anisms of orthotropic composites. The added final failure condition distinguishes between in-plane and out-of-plane loading. An additional damage variable guarantees a mesh objective response even at catastrophic failure and enables relatively fast simulations. Conventional material charac-terisation tests were used to validate some of the assump-tions used in the model. The predictive capabilities of the model have been evaluated for the crushing of a flat unidirectional specimen. The results were in accordance with experiments in terms of stress response, energy ab-sorption and damage morphology, and showed a very low mesh dependency.
The presented model can be used to drive the design of crashworthy composite structures, reducing the devel-opments costs considerably. Future work will focus on combining the current model with a fibre kinking model developed by the authors as well as the study of more complex structures where delaminations and large defor-mations play an important role. The strain rate effects
0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 Compressive strain,
e
33(%) D am age, Onset of damage 3 2 a a b c Immediately after specimen failure ~ 51° 2 3 c Immediately before specimen failure 3 2 b 3 2Figure 9: Distribution and evolution of the matrix damage variable dmcin the DERA compressive specimens for an element at the
centre of the specimen.
were not considered in the present model but are being addressed in a parallel project.
5. Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publica-tion of this article.
6. Acknowledgments
This work was supported by the Swedish Energy Agency [project number 34181-2]; and FFI/VINNOVA [dnr 2016-04239]. The discussions with Renaud Gutkin at Volvo Car Corporation, Martin Fagerstr¨om at Chalmers University of Technology, and Gaurav Vyas at Swerea SICOMP AB are gratefully acknowledged.
References
[1] T. Jones, Assessment of fuel economy technologies for light-duty vehicles. Non-Engine Technologies, Washington, DC: The National Academies Press., 2011. doi:10.17226/12924. [2] J. J. Carruthers, A. P. Kettle, A. M. Robinson, Energy
absorp-tion capability and crashworthiness of composite material struc-tures: a review, Applied Mechanics Reviews 51 (1998) 635–649. [3] G. Barnes, I. Coles, R. Roberts, Crash safety assurance strate-gies for future plastic and composite intensive vehicles (PCIVs), Tech. Rep. DOT-VNTSC-NHTSA-10-01, U.S. Department of Transportation (2010). 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 Experiments FE Specimen failure Compressive strain, e33(%)
n
32 Damage No damageFigure 10: Evolution of the Poisson’s ratio ν32with matrix damage.
[4] S. T. Pinho, P. P. Camanho, M. F. De Moura, Numerical sim-ulation of the crushing process of composite materials, Interna-tional Journal of Crashworthiness 9 (3) (2004) 263–276. [5] L. N. S. Chiu, B. G. Falzon, R. Boman, B. Chen, W. Yan, Finite
element modelling of composite structures under crushing load, Composite Structures 131 (2015) 215–228.
[6] C. McGregor, N. Zobeiry, R. Vaziri, A. Poursartip, X. Xiao, Calibration and validation of a continuum damage mechanics model in aid of axial crush simulation of braided composite tubes, Composites Part A: Applied Science and Manufacturing 95 (2017) 208–219.
[7] Jackson KE, Fasanella EL and Littell JD. Development of a con-tinuum damage mechanics material model of a graphite-kevlar hybrid babric for simulating the impact response of energy ab-sorbing subfloor concepts. In: AHS International Annual Fo-rum & Technology Display; 73rd, Fort Worth, TX, United States, May 9–11 2017, 12p.
[8] P. Ladev`eze, E. Le Dantec, Damage modelling of the elementary ply for laminated composites, Composites Science and Technol-ogy 43 (3) (1992) 257–267.
[9] C. Gonz´alez, J. LLorca, Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: Micro-scopic mechanisms and modeling, Composites Science and Tech-nology 67 (13) (2007) 2795–2806.
[10] L. G. Melin, J. Neumeister, K. B. Pettersson, H. Johansson, L. E. Asp, Evaluation of four composite shear test methods by digital speckle strain mapping and fractographic analysis, Journal of Composites Technology & Research 22 (3) (2000) 161–172.
[11] R. Gutkin, S. T. Pinho, Combining damage and friction to model compressive damage growth in fibre-reinforced compos-ites, Journal of Composite Materials 49 (20) (2015) 2483–2495. [12] LS-Dyna Keyword User Manual – Volume II. Material Mod-els., Tech. Rep. Version R7.0, Livermore Software Technology Corporation (2013).
[13] Deleo F, Wade B, Feraboli P et al. Crushing of compos-ite structures: experiment and simulation. In: Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Struc-tural Dynamics, and Materials Conference, Palm Springs, CA, United States, 7 May 2009, 22p.
[14] H. A. Israr, S. Rivallant, C. Bouvet, J.-J. Barrau, Finite ele-ment simulation of 0/90 CFRP laminated plates subjected to crushing using a free-face-crushing concept, Composites Part A: Applied Science and Manufacturing 62 (2014) 16–25.
[15] D. Shi, X. Xiao, An enhanced continuum damage mechanics model for crash simulation of composites, Composite Structures 95 (2018) 208–219.
Crushing plate (with friction) A0/2 Bottom support Side support (with friction) Symmetry plane Coarser mesh (0.3 mm) Finer mesh (0.17 mm) Y X Z
14.6 1.8 22 4.1 45° specimen 90° specimen fibresFigure 11: Geometry with dimensions (in millimetres) and FE model of the crush specimens.
[16] R. Gutkin, S. Costa, R. Olsson, A physically based model for kink-band growth and longitudinal crushing of composites un-der 3D stress states accounting for friction, Composites Science and Technology 135 (2016) 39–45.
[17] S. Costa, R. Gutkin, R. Olsson, Mesh objective implementa-tion of a fibre kinking model for damage growth with fricimplementa-tion, Composite Structures 168 (2017) 384–391.
[18] S. T. Pinho, L. Iannucci, P. Robinson, Physically based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development, Com-posites Part A: Applied Science and Manufacturing 185 (2006) 774–785.
[19] T. Bru, R. Olsson, R. Gutkin, G. M. Vyas, Use of the Iosipescu test for the identification of shear damage evolution laws of an orthotropic composite, Composite Structures 174 (2017) 319– 328.
[20] A. Puck, H. Sch¨urmann, Failure analysis of FRP laminates by means of physically based phenomenological models, Compos-ites Science and Technology 58 (1998) 1045–1067.
[21] T. Bru, P. Hellstr¨om, R. Gutkin, D. Ramantani, G. Peterson, Characterisation of the mechanical and fracture properties of a uni-weave carbon fibre/epoxy non-crimp fabric composite, Data in Brief 6 (2016) 680–695.
[22] S. Swanson, M. Messick, Z. Tian, Failure of Carbon/Epoxy lamina under combined stress, Journal of Composite Materials 21 (7) (1987) 619–630.
[23] P. D. Soden, M. J. Hinton, A. S. Kaddour, Biaxial test results for strength and deformation of a range of E-glass and carbon fibre reinforced composite laminates: failure exercise benchmark data, Composites Science and Technology 62 (12–13) (2002) 1489–1514.
[24] H. Molker, D. Wilhelmsson, R. Gutkin, L. E. Asp, Orthotropic criteria for transverse failure of non-crimp fabric-reinforced com-posites, Journal of Composite Materials 50 (18) (2016) 2445– 2458.
[25] S. T. Pinho, L. Iannucci, P. Robinson, Physically based fail-ure models and criteria for laminated fibre-reinforced compos-ites with emphasis on fibre kinking. Part II: FE implementa-tion, Composites Part A: Applied Science and Manufacturing 37 (2006) 766–777.
[26] W. Tan, B. G. Falzon, Modelling the crush behaviour of ther-moplastic composites, Composites Science and Technology 134 (2016) 57–71.
[27] M. Benzeggagh, M. Kenane, Measurement of mixed-mode de-lamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus, Composites Science and Technology 56 (4) (1996) 439–449.
[28] J. Sch¨on, Coefficient of friction of composite delamination sur-faces, Wear 237 (1) (2000) 77–89.
[29] Z. Baˇzant, M. Jir´asek, Nonlocal Integral Formulations of Plas-ticity and Damage: Survey of Progress, Journal of Engineering Mechanics 128 (11) (2002) 1119–1149.
[30] L. Jia, L. Yu, K. Zhang, M. Li, Y. Jia, B. R. K. Blackman, J. P. Dear, Combined modelling and experimental studies of failure in thick laminates under out-of-plane shear, Composites Part B: Engineering 105 (2016) 8–22.
[31] R. F. Ferguson, M. J. Hinton, M. J. Hiley, Determining the through-thickness properties of FRP materials, Composites Sci-ence and Technology 58 (9) (1998) 1411–1420.
[32] C. Stocchi, P. Robinson, S. Pinho, A detailed finite element investigation of composite bolted joints with countersunk fas-teners, Composites Part A: Applied Science and Manufacturing 52 (2013) 143–150.
[33] R. Olsson, J. Iwarsson, L. G. Melin, A. Sj¨ogren, J. Solti, Exper-iments and analysis of laminates with artificial damage, Com-posites Science and Technology 63 (2) (2003) 199–209. [34] T. Bru, P. Waldenstr¨om, R. Gutkin, R. Olsson, G. M. Vyas,
Development of a test method for evaluating the crushing be-haviour of unidirectional laminates, Journal of Composite Ma-terials 51 (29) (2017) 4041–4051.
[35] Bru T, Olsson R, Vyas GM et al. Validation of a novel model for the compressive response of FRP: experiments with differ-ent fibre oridiffer-entations. In: Proceedings of ICCM-21 Conference, Xi’an, China, 20–25 August 2017, 9p.
[36] L. Grauers, R. Olsson, R. Gutkin, Energy absorption and dam-age mechanisms in progressive crushing of corrugated NCF laminates: Fractographic analysis, Composite Structures 110 (2014) 110–117.
dam ag e
0
20
40
60
80
100
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
δ
(mm)
F/
A
0(M
Pa)
FE 0.17 mm
FE 0.30 mm
FE 0.30 mm + delam.
Experiments
Experiment:
FE:
0.00 0.17 0.33 0.50 0.67 0.83 1.00Figure 12: Experimental and numerical crushing response of the [45◦]10crush specimen.
0 10 20 30 40 50 60 FE 0.3 mm FE 0.17 mm Exp. Max. Min. Sp eci fi c en er g y ab so rp ti o n (kJ/ kg )
Figure 13: SEA of the experiments vs. simulations of the [45◦]10
specimen . 0 10 20 30 40 50 60 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 FE 0.30 mm + delam. Experiments d(mm) F /A 0 (M Pa)
Figure 14: Experimental and numerical crushing response of the [90◦]
