A methodology to simulate low velocity impact and compression after impact in large composite stiffened panels

Accepted Manuscript A methodology to simulate low velocity impact and compression after impact in large composite stiffened panels A. Soto, E.V. González, P. Maimí, J.A. Mayugo, P.R. Pasquali, P.P. Camanho PII: DOI: Reference:

S0263-8223(18)30728-1 https://doi.org/10.1016/j.compstruct.2018.07.081 COST 9999

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

20 February 2018 10 June 2018 20 July 2018

Please cite this article as: Soto, A., González, E.V., Maimí, P., Mayugo, J.A., Pasquali, P.R., Camanho, P.P., A methodology to simulate low velocity impact and compression after impact in large composite stiffened panels, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.07.081

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A methodology to simulate low velocity impact and compression after impact in large composite stiffened panels A. Sotoa,∗, E. V. Gonz´aleza , P. Maim´ıa , J. A. Mayugoa , P. R. Pasqualib , P. P. Camanhoc,d a AMADE,

Polytechnic School, University of Girona, Campus Montilivi s/n, 17071 Girona, Spain S.A., Av. Brigadeiro Faria Lima 2170, 12227-901 S£o Jos´e dos Campos, Brazil c DEMec, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal d INEGI, Instituto de Ciłncia e Inova£o em Engenharia Mec¢nica e Engenharia Industrial, Rua Dr. Roberto Frias, 400, 4200-465 Porto, Portugal b Embraer

Abstract Low velocity impact events significantly reduce the mechanical performance of composite structures even though the damage might be barely visible. Numerical simulations can be used to understand and improve the damage resistance and tolerance of composite structures. However, numerical simulations are usually computationally intensive and their application in large composite structures is limited. Furthermore, the numerical models require many parameters that affect their efficiency, accuracy, objectivity and robustness. The present work describes a methodology to simulate low velocity impact and compression after impact which is applied to a composite stiffened panel undergoing visible impact damage. The key definitions are discussed and special attention is devoted to the computational efficiency. The numerical results are compared with experimental data, and the suitability of the proposed methodology is discussed. Keywords: Low velocity impact, compression after impact, finite element analysis, composite structures

∗ Corresponding

author Email address: [email protected] (A. Soto )

Preprint submitted to Composites Structures

June 10, 2018

1

1. Introduction

2

Low Velocity Impact (LVI) events are a major concern in the design of aerospace structures. LVI can lead to Visible

3

Impact Damage (VID) or Barely Visible Impact Damage (BVID). VID involves delamination but, more specifically,

4

matrix cracking and fiber breakage are the dominant damage mechanisms, whereas BVID usually involves large de-

5

laminations with small surface indentation. BVID is particularly dangerous in structures made of composite laminates

6

because damage is difficult to observe through visual inspections and the strength of the structure could have been

7

significantly reduced. Assessing impact damage on composite structures requires extensive experimental campaigns

8

that are time-consuming and costly. Consequently, modeling impact events has been the focus of many researchers

9

over the years in an attempt to reduce the number of physical tests required. In the literature, the problem is treated

10

either analytically or by means of numerical methods.

11

Analytic impact models are a powerful tool that provide a rapid assessment of the type of damage. Some relevant

12

contributions are the work from Abrate [1] and Olsson [2, 3]. However, analytical impact models are usually restricted

13

to the elastic regime and specific geometry and the interaction among damage mechanisms and their progression is

14

not taken into account. For damage tolerance assessments such as Compression After Impact (CAI), initial damage

15

has to be assumed to analytically predict the residual strength [4].

16

Virtual testing by means of non-linear Finite Element Analysis (FEA) is an alternative that can tackle any geom-

17

etry and boundary conditions. The stress field takes into account the interaction between damage mechanisms and

18

their damage progression. In addition, the impacted model can be used for the subsequent damage tolerance analysis.

19

Numerical impact models can be classified as strength-based models [5], fracture mechanics based models [6] and

20

Continuum Damage Mechanics (CDM) based models. CDM based models proved their suitability in several studies

21

for the LVI prediction [7–10] and LVI with the sequential CAI [11–15] in composite laminates at the coupon level.

22

The CDM framework [16, 17] combines failure criteria with fracture mechanics to predict both damage initiation

23

and progression through coupling internal damage variables, which represent each damage mechanism. However,

24

these numerical analyses are computationally intensive which limits their application to structural components. Fur-

25

thermore, they require many parameters and definitions that can affect the numerical results as well as the model

26

objectivity and robustness, which sometimes are not clearly reported.

27

The finite element type together with the CDM model for ply modeling play an important role in the model’s

28

efficiency and accuracy. Different finite element types are available depending on the desired detail of the stress field

29

(e.g. solid elements or shell elements). The Cohesive Zone Model (CZM) is widely used to model delamination

30

and can be used by means of cohesive elements or cohesive surfaces. Cohesive elements use the nodal information

31

to calculate the openings, while cohesive surfaces calculate the openings through a contact algorithm. The selected 2

32

cohesive interaction technology and CZM are also important for the model’s accuracy and efficiency. Model ob-

33

jectivity can be ensured by using appropriate element sizes as well as defining a mesh regularization algorithm that

34

avoids mesh dependency due to strain localization. During the LVI and CAI simulation the continuum and cohesive

35

finite elements can suffer severe distortions due to damage which can eventually abort the simulation. The model’s

36

robustness depends on the type of finite element and cohesive interaction technologies but it depends on the element

37

deletion criterion in particular.

38

Undoubtedly, the simulation of both LVI and CAI tests for large structures is more interesting from an industrial

39

point of view. Some examples in the literature are found to deal with LVI simulation of composite laminate structural

40

components [18–21]. To the authors’ knowledge, Psarras et al. [22] is the only reported work on LVI and CAI

41

simulation of large composite structures.

42

Johnson et al. [21] performed LVI on WR E-glass/Derakane 8084 vinyl - ester marine composite panels. The

43

simulated impact energies ranged from 195 to 7200 J, which fall within the range of VID. Three different plates were

44

simulated: 228.6 × 177.8 × 6.35 mm3 , 1073.2 × 768.4 × 19 mm3 and 1370 × 1370 × 38 mm3 . Solid elements

45

were used with a CDM model inspired by the plane-stress model from Williams et al. [23], which was extended to

46

account for transverse shear damage. The damage model was fed with data obtained from cyclic tests. The damage

47

variable was limited to a value of 0.9 without element deletion. No mesh regularization algorithm was mentioned

48

to be used. The element sizes employed were 1.6 and 3.2 mm for the smallest and the largest plates, respectively.

49

Delamination was modeled by means of cohesive elements but not all the interfaces for delamination were considered

50

for the medium and large panel. The number and location of the interfaces for delamination can alter the predicted

51

damage evolution [24]. The stiffness and the dissipated energy were not well captured. However, the modeling

52

challenge due to the panel dimensions and the impact energies, which involved large numerical models with a large

53

amount of fiber breakage and delamination (VID), must be acknowledged.

54

Faggiani and Falzon [18] simulated an LVI of 15 J on a 450 mm by 375 mm Carbon Fiber Reinforced Polymer

55

(CFRP) panel with three I-stiffeners, which falls within the BVID range. The impact was located in between two

56

stiffeners. A refined mesh region at the impact location with all the plies and interfaces for delamination was modeled

57

with reduced integration 3-D solid elements and zero-thickness cohesive elements. The remaining of the panel did

58

not account for delamination. It was discretized with a single continuum shell element through the thickness which

59

was coupled to the refined region through tie constraints. The element sizes used in the model were not reported. The

60

CDM model from Donadon et al. [25] and the CZM from Camanho et al. [26] were used to model intralaminar and

61

interlaminar damage, respectively. A friction value of 0.5 was assumed for the ply-to-ply contact and 0.3 for the skin-

62

to-impactor contact. The damage variables were limited to a maximum value of 0.99 without element deletion to avoid 3

63

excessive element distortion issues and the crack band model was used [27] for mesh regularization. Delamination

64

was the main damage mechanism and only a small amount of fiber breakage was predicted due to the low impact

65

energy. The dissipated energy was under-predicted by approximately 12% but the projected delamination area and the

66

overall force - time response were in good agreement with the experimental ones.

67

Riccio et al. [19] published a numerical study on a 246 mm by 368 mm stiffened composite panel with two

68

omega shaped stringers together with experimental results. The impact was located close to a stiffener. The simulated

69

impact energies were 15 J and 25 J, which fall within the BVID range. Continuum shell elements were used with

70

the Hashin CDM model [28] together with the crack band model [27]. The element size employed was not reported.

71

Finite thickness (i.e. 0.01 mm) cohesive elements with quadratic nominal stress initiation criterion and a power law

72

energy based propagation criterion were used to model delamination. The elements were deleted even though the

73

criteria employed was not reported. A friction value of 0.5 was assumed for all contact pairs. The overall force - time

74

response was captured. However, the energy dissipated was over-predicted by 20 % and 50 % in the 15 J and 25 J

75

impacts, respectively.

76

Psarras et al. [22] investigated the compressive residual strength of large curved stiffened panels after multi-site

77

impacts. Two different laminate thicknesses and six multi-site impact scenarios were considered. The panel size was

78

1200 mm by 800 mm and the material used was T800/M21. The impact energies ranged from 25 to 58 J. Continuum

79

shell elements with the plane-stress intralaminar damage model proposed by Iannucci et al. [29] were used. Cohesive

80

elements were placed among sublaminates to account for interlaminar damage. No more details on the model’s

81

definitions such as the element size or element deletion were reported. While the reported numerical predictions of

82

the CAI strength were in good agreement for the thinnest panels (i.e. 1 - 3 %), the strength was over-predicted for

83

the thickest ones (i.e. 19 - 32 %). The experiments demonstrated the influence impact energy, location and number of

84

multi-site impacts have on the compressive residual strength.

85

Recently, Sun and Hallett [20] performed a 15 J LVI simulation on a 450 mm by 375 mm stiffened panel made

86

of HTA/6176C, which falls within the BVID range. Plies were modeled with reduced integration solid elements and

87

cohesive elements were inserted vertically into each ply to model matrix cracking. Cohesive elements were also used

88

to model delamination but fiber breakage was neglected. A shell to solid coupling was used to model the impact

89

region in detail for all the plies and interfaces, while single conventional shell elements through the thickness were

90

used out of the impact region for the sake of computational efficiency. The element size was 0.2 mm in the refined

91

region. The modeling approach was firstly validated with experimental data and compared with a full 3-D numerical

92

model for ASTM standard coupon specimens [30]. Computational time savings of up to 50 % were reported with the

93

shell-to-solid coupling approach in comparison with the full 3-D model. The results with the shell-to-solid coupling 4

94

approach tend to lead larger contact times and suffer more oscillations. However, the reported force-time response of

95

the stiffened panel modeled with the solid/shell approach was in agreement with experimental data.

96

The objective of the present work is to define an efficient modeling methodology able to predict both the LVI and

97

CAI response in large composite stiffened panels. Abaqus/Explicit [28] is used. The key parameters for defining the

98

numerical model are clearly described and justified. Firstly, the methodology is applied at the coupon level within

99

the range of BVID. The aim is to perform a numerical benchmark that compares the computational performance and

100

accuracy of different finite element types (i.e. solid, continuum shell and conventional shell) and cohesive interaction

101

technologies (i.e. cohesive elements and cohesive surfaces). In Section 3, an LVI and a CAI simulation of an aircraft

102

stiffened panel that undergoes VID is performed using the selected methodology and results in the better compromise

103

between accuracy and computational effort. The results are compared with the experimental data. Finally, some

104

conclusions and future work are discussed.

105

2. Modeling methodology

106

During an impact event both intralaminar (i.e. matrix cracking, fiber breakage) and interlaminar (i.e. delamination)

107

damage mechanisms take place. Several studies [7–14] support that the meso-scale provides enough detail on the

108

stress field for LVI and CAI predictions on composite plates. The aforementioned studies support the use of CDM and

109

CZM as frameworks to model intralaminar and interlaminar, respectively. However, numerical impact models require

110

the definition of several parameters that affect the robustness, efficiency, accuracy and objectivity of the numerical

111

model. In this section, a methodology to define the relevant modeling parameters is described and discussed.

112

2.1. Intralaminar damage modeling

113

There are several numerical studies on LVI predictions of composite laminates in which different CDM models

114

are used to model the intralaminar damage [8, 13, 21, 23, 25, 29, 31–35]. While it is out of the scope of this work

115

to review the existing CDM for composite laminates, the main differences among them stem from whether they are

116

three-dimensional or they assume plane-stress conditions, from the damage initiation and the damage propagation

117

criteria, the behaviour under shear loading or the mesh regularization algorithm used to address strain - softening.

118

In the literature, impact models assuming plane-stress conditions [14, 15, 29, 36] or three dimensional stress state

119

[8, 10, 11, 13, 18] are used depending on the desired compromise between computational effort and detailed stress

120

field. Out-of-plane stresses are neglected under plane-stress conditions. However, some examples in the literature

121

show that the structural response of both LVI and the sequential CAI can be predicted [14, 15].

122

The present work uses the CDM model proposed by Maim´ı et al. [33, 34] because it is a thermodynamically con-

123

sistent model with physically based damage activation functions [37]. The model assumes plane-stress conditions and 5

124

its predictive capability has been proven under different loading conditions [11, 34, 38]. In-plane isotropic plasticity

125

has been assumed for the sake of simplicity (see Fig. 1) to account for the matrix dominant behaviour under shear

126

loading, no experimental data was available to feed a kinematic plastic model. The fiber damage evolution has to be

127

defined in the model. The model allows the fiber traction separation law to be introduced (see Fig. 2). [Figure 1 about here.]

128

129

The model is used through an user-written subroutine (VUMAT) in Abaqus/Explicit [28]. The crack band model

130

[27] is used to deal with strain localization. Four surfaces for damage activation are defined as F N = φN − rN 6 0.

131

Each damage activation function (F N ) accounts for one damage mechanism: longitudinal (N = XT ) and transverse

132

(N = YT ) tension, longitudinal (N = XC ) and transverse compression (N = YC ). The loading functions (φN ) depend

133

on the effective stresses and material properties. The internal variables (rN ) set the maximum value that the loading

134

functions can achieve before damage propagation. They are equal to one while the material is undamaged. When

135

the loading functions are larger than one, damage begins and the internal variables are updated to the new damage

136

e threshold. Additionally, to account for plasticity under shear loading the yield function is defined by F P = |γ12 |−

137

i e i KP γ12 − S LP /G12 6 0, where γ12 , S LP , G12 , KP , γ12 are the elastic shear strain, the yield stress, the shear elastic

138

modulus, a plastic parameter and the isotropic hardening variable, respectively.

139

The model damage activation functions are based on the LARC03 failure criteria [37] but neglecting the out-of-

140

plane-stresses. Catalanotti et al. [39] proposed new 3-D failure failure criteria for longitudinal and transverse failure

141

accounting for the effect of ply thickness on material strength. The effect of ply thickness on material strength for

142

transverse tensile and compression is approximated here by the analytical expressions from Camanho et al. [40]

143

and Maim´ı et al. [41], respectively. The in-situ shear strengths from [40] were adapted to the linear elasto-plastic

144

behaviour considered in the constitutive model for the sake of consistency.

145

146

Eqs. (1) - (3) are the in-situ shear strengths for thick plies (S isLthick ), thin plies (S isLthin ) and outer thin plies (S isLout ), respectively. q S isLthick =

(KP + 1)(2S 2L KP + 2S 2L − S 2LP ) KP + 1

(1)

147

q S isLthin

=

h ply (KP + 1)(πS 2LP h ply + 8G12GS L KP ) √ (KP + 1)h ply π

(2)

S isLout

q h ply (KP + 1)(πS 2LP h ply + 4G12GS L KP ) = √ (KP + 1)h ply π

(3)

148

6

149

150

where h ply is the ply thickness. The in-situ strength of an embedded ply (S isLemb ) is the maximum value of Eqs. (1) - (2).

151

In Section 2.5 a benchmark is made in which different finite element types (i.e. solid, continuum shell and

152

conventional shell) are used to compare their performance. The constitutive model is adapted to be used with 3-

153

D solid elements by defining the damage variables d3 , d4 and d5 as a function of other damage variables as done in

154

previous works [10, 11, 38]: d3 (rXC , rYC ) = 1 − [1 − d1 (rXC )][1 − d2 (rYC )]

(4)

d4 (rYT ) = d6 (rYT )

(5)

d5 (rXT ) = d1 (rXT )

(6)

155

156

157

The damage variables (dN ) depend on the elastic and fracture properties as well as on the traction separation law

158

shape. D´avila et al. [42] found that linear softening was insufficient to correctly predict the damage initiation and

159

propagation of cross-ply Compact Tension (CT) specimens. Multiple damage mechanisms occur during the fracture

160

of a composite laminate (e.g. fiber-bridging, fiber pull-outs), which are embedded within the traction separation law

161

shape at the macro-scale. Matrix related damage mechanisms usually have a relatively small Failure Process Zone

162

(FPZ) in comparison to the specimen size in conventional composite laminates, while the FPZ related to fiber damage

163

mechanisms can span several millimeters [43]. Consequently, the fiber traction separation law shape has an effect

164

on the structural response, as reported in [14, 42, 44]. The model assumes bi-linear softening for the longitudinal

165

directions and linear softening for the transverse and shear directions (see Fig. 2). The parameters fXM and HXM define

166

the traction separation law shape (see Fig. 2) together with the strengths (i.e. X M , Y M , S L ) and fracture toughness (i.e.

167

G XM , GYM , GS L ). They can be obtained following the method presented in [44]. Alternatively, they can be fitted to the

168

load displacement curve of the CT and Compact Compression (CC) test specimens through FEA [14]. [Figure 2 about here.]

169

170

2.2. Interlaminar damage modeling

171

The CZM assumes that the entire FPZ is lumped into the crack plane. CZM are very efficient in situations where

172

the crack path is known beforehand and it is an accepted approach to model delamination. The main differences

173

between interlaminar CZM’s are the initiation and propagation criteria under mixed mode conditions and the cohesive

174

law shape considered.

175

The built-in CZM from Abaqus/Explicit [28] is used in all the models presented in this study. In Section 2.5

176

is used with the Abaqus built-in zero-thickness cohesive elements (COH3D8) and cohesive surfaces to compare the 7

177

computational performance and numerical results between both cohesive interaction technologies while in Section 3

178

zero-thickness cohesive elements are used to model delamination.

179

The quadratic traction criterion is used for damage initiation under mixed-mode loading. Orifici et al. [45] re-

180

viewed the existing formulations to model delamination based on the CZM approach. Later, Turon et al. [46] proposed

181

a thermodynamically consistent damage model in which the damage initiation criterion is formulated according to the

182

damage propagation (i.e. Benzeggagh-Kenane criterion [47]). The damage initiation leads to similar results as the

183

quadratic traction criterion [46]. Failure criteria that predict delamination initiation under mixed-mode conditions

184

have not been fully validated due to the lack of experimental data [46]. Nevertheless, the quadratic criterion is more

185

accepted than the maximum stress criterion. Benzeggagh-Kenane criterion [47] is used as a damage propagation cri-

186

terion under mixed mode loading. Camanho et al. [26] argue that the Benzeggagh - Kenane criterion [47] is suitable

187

for the critical energy release rate under mixed-mode ratio in epoxy composites.

188

Damage evolution is controlled through the cohesive law shape, which is assumed linear here. The FPZ in com-

189

posites delamination under pure mode I is usually small. Thus, the effect of the cohesive law shape on the structural

190

response is negligible as Alfano reported [48]. The FPZ in composites delamination under pure mode II is known

191

to be larger than in mode I [49] but few studies address what the cohesive law shape is under mode II delamination

192

[50, 51].

193

2.2.1. Density of cohesive interactions

194

The explicit dynamics procedure solves the boundary - value problem as a wave propagation problem. A bounded

195

solution is obtained when the simulation time increment (4t) is below the Stable Time Increment (STI), which is a

196

stability limit; the minimum time that a dilatation wave takes to move across any element of the model.

197

The STI of a continuum finite element (4telem ) is given by: r 4telem ≤ l α ∗

198

ρ E

where ρ is the material density, E the maximum component of the constitutive tensor, and α =

(7)

p

1 + ζ 2 − ζ with

199

ζ being the fraction of critical damping that prevents finite elements from ringing or even collapsing in the highest

200

element frequency [28].

201

The shortest characteristic length (l∗ ) controls the stable time increment. In composite laminates the smaller

202

characteristic length is related to ply thickness. When volumetric elements are used (e.g. solid and continuum shell

203

elements) ply thickness governs the stable time increment. In the case of conventional shell elements, l∗ is controlled

204

by the in-plane finite element dimensions only. 8

205

coh The STI of a zero-thickness cohesive element (4tint ) is:

r coh 4tint

≤

ρ¯ Kcoh

(8)

206

where ρ¯ and Kcoh are the surface density and the penalty stiffness value of the zero-thickness cohesive elements,

207

respectively. Eq. (8) also applies for Abaqus/Explicit [28] Surface Elements (SFM). In Abaqus/Explicit [28], cohesive

208

surfaces cannot be used with conventional shell elements in a straightforward manner when multiple plies are stacked

209

over each other due to poor kinematic description [15]. SFM allow a sheet of nodes without structural behaviour to

210

be created, which in turn, allows conventional shell elements to be used with cohesive surfaces as shown in Gonz´alez

211

et al. [15]. However, they require a mass definition (i.e. density by unit surface) because in a dynamics system all

212

the nodes need inertia and mass terms. In fact, their definition is analogous to zero-thickness cohesive elements. The

213

reader is referred to Gonz´alez et al. [15] for more details on the use of conventional shell elements together with

214

cohesive surfaces in Abaqus/Explicit [28].

215

con The STI of a contacting surface (4tint ) is controlled by:

r con 4tint

≤

ρhelem Kn

(9)

216

where ρ is the material density, helem the element thickness to which the contacting surface belongs and Kn the normal

217

contact penalty stiffness value.

218

coh con The computational time of an explicit simulation is governed by the STI; 4tmin = min(4telem , 4tint , 4tint ). The

219

coh con interaction (4tint , 4tint ) and element (4telem ) STI are not generally the same. It is one of these that controls the

220

model’s STI.

221

The mass of zero-thickness cohesive elements or surface elements should be zero because they have no volume

222

but a dynamic system to be solved requires all nodes in the model to be associated with inertial and mass terms. Thus,

223

the surface density (ρ) ¯ of zero-thickness cohesive elements or SFM, which is a numerical parameter, can be defined

224

to reduce the difference between the STIs and consequently reduce simulation time. For instance, this can be done by

225

distributing the mass of the laminate among the continuum elements and the cohesive elements [14].

226

The optimum density (ρ) and density by unit surface (ρ) ¯ is obtained by equaling the STI expressions:

coh con 4telem ' 4tint ' 4tint

227

(10)

This can be done with zero-thickness cohesive elements with any element type and SFM (i.e. cohesive surfaces) 9

228

with conventional shell elements. In the case of using cohesive surfaces with continuum shell or solid elements it

229

cannot be done because the density by unit surface from the cohesive surface is computed based on the density of the

230

underlying continuum element and its thickness/characteristic element length.

231

2.2.2. Penalty stiffness

232

The CZM requires defining a penalty stiffness value. The material is perfectly bonded before crack initiation.

233

Thus, the penalty stiffness value should be as large as possible so as not to contribute to the global compliance of the

234

structure.

235

Daudeveville et al. [52] defined the cohesive penalty stiffness (Kcoh ) as the ratio of the interface elastic modulus

236

with the interface thickness, which is considered to have a thickness of 0.005 - 0.01 mm [11, 38, 53]. Some authors

237

[52, 54] consider the interface thickness as one fifth of the adjacent sublaminate. Turon et al. [53] recommended

238

defining the penalty stiffness as Kcoh ≥

239

elastic modulus and t the adjacent sublaminate thickness. The aforementioned definitions lead to penalty stiffness

240

values in the order of 5 × 105 - 1 × 106 N/mm3 , which have been used with satisfactory results in [13, 26, 53]. It is

241

also worth noting that the penalty stiffness should have the same value for mode I and mode II to avoid changes in the

242

local mode ratio under mixed-mode damage evolution unless a consistent formulation is used [55].

50E22 t

to avoid affecting the compliance of the system; E22 is the transverse

243

In explicit simulations, the penalty stiffness does not only affect the accuracy, but also the computational perfor-

244

mance through the STI (see Eq. (8)). It is difficult to have a closed-form expression to estimate the penalty stiffness

245

that takes into account the number of interfaces considered, the elastic properties of the surrounding plies and the

246

loading conditions. Instead, a sensitivity analysis can be done on the elastic regime to choose the minimum penalty

247

stiffness value that converges to the same elastic response [14]. Thus, compliance is not added into the model while

248

mitigating the detrimental effect on the STI.

249

2.3. Element size

250

The element size of the model is important to ensure correct intralaminar and interlaminar energy dissipation.

251

On the one hand, CDM should avoid snap-back at the element level for correct intralaminar energy dissipation

252

[27] even when a mesh regularization algorithm is used. If a mesh refinement is unfeasible, the snap-back in the

253

constitutive model can be avoided by reducing the corresponding strength [27] but this can affect the accuracy of

254

the prediction in the damage initiation. In fact, snap-back should be avoided within the first branch of the traction

255

separation law because of its importance on the structural strength [56, 57]. Therefore, the criterion adopted is to

256

assess the maximum element size for each damage mechanism (N) that avoids snap-back within the first branch of

10

257

the traction separation law by means of Eq. (11) [35].

l∗ 6

EN HN

(11)

258

where l∗ , E N and HN are the characteristic element length, the elastic modulus, and the first branch slope of the traction

259

separation law for a given damage mechanism (N), respectively. In the case of the linear softening law, HN = XN2 /2G N .

260

On the other hand, the interlaminar mesh discretization has to be chosen according to the interlaminar FPZ size

261

because the damage initiation and propagation are not correctly computed if the FPZ is not discretized with a minimum

262

of three elements along the FPZ [49, 58, 59]. A conservative estimation of the interlaminar FPZ is obtained by using

263

the expressions from Soto et al. [49] and considering the smallest ply thickness modeled.

264

265

Finally, the element size selected is that which ensures correct intralaminar and interlaminar dissipation. 2.4. Finite element and interaction type

266

Different finite element and cohesive interaction technologies are possible in FEA. The CZM using cohesive

267

elements or cohesive contact surfaces is widely used in the literature to model delamination. Zero-thickness cohesive

268

elements [7, 10, 14, 18] and finite thickness cohesive elements [11, 19, 38] have been mostly used to model impact

269

in comparison to cohesive surfaces [10, 13, 15]. Cohesive surfaces naturally handle large deformations and non-

270

conforming meshes, however, they are more expensive computationally due to the contact tracking algorithm [10].

271

Three dimensional solid elements with reduced integration are mostly used in impact simulations [7–13, 18, 21,

272

25, 38] not only for their computational efficiency, but also to alleviate locking pathologies while taking into account

273

the three dimensional stress state. Mechanisms such as matrix cracking induced delamination or the in-situ effect can

274

be naturally captured with appropriate failure criteria and accurate stress fields. However, several elements through

275

the thickness of every ply as well as aspect ratios close to the unit may be needed for accurate stress fields.

276

Continuum shell elements have been used as an efficient alternative to solid elements for LVI simulation [19,

277

22]. Continuum shell elements (also known as solid-shell finite elements) have the same continuum topology as

278

three dimensional solid elements but with shell-like kinematics. Contact takes place on the actual shell surface and

279

thickness variations based on physical nodes are accounted for. Furthermore, they can be naturally connected with

280

solid elements because both have displacements as degrees of freedom. Continuum shell elements allow larger in-

281

plane to thickness aspect ratios than solid elements do. Nonetheless, the STI of continuum shell elements in explicit

282

simulations can be controlled by their element thickness (Eq. (7)).

283

Conventional shell elements are structural elements based on plate theory. Plate theory takes advantage of the

284

small thickness compared to the planar dimensions to reduce the full three-dimensional solid mechanics problem to 11

285

a two-dimensional one. This involves a reduction in terms of computational cost while being suitable for bending

286

problems because of the rotational degrees of freedom. Furthermore, in explicit simulations the element thickness

287

does not detrimentally affect the STI of the model because it is not accounted for in Eq. (7). However, conventional

288

shell elements cannot provide accurate results for the transverse shear and normal strains [60]. Despite the poorer

289

transverse description than solid or advanced continuum shell elements formulations, conventional shell elements are

290

considered a well-balanced combination of accuracy and computational efficiency [61]. However, contact related

291

issues and poor kinematic description to account for delamination is reported when stacking several conventional

292

shell elements over each other [15, 62, 63]. This has probably limited its application for simulating the LVI and

293

CAI simulation. Recently, conventional shell elements have been applied to predict the structural response of both

294

LVI and CAI tests in composite thin ply laminates [14] and standard laminates [15]. Gonz´alez et al. [15] employed

295

conventional shell elements with cohesive surfaces by using tie constraints between shell elements and SFM (Fig.

296

3a), which are a sheet of nodes without structural behaviour. Soto et al. [14] used tie constraints between shell

297

elements and zero-thickness cohesive elements [14] (Fig. 3b). Both strategies [14, 15] allow for meshing along the

298

fiber direction and account for friction effects. [Figure 3 about here.]

299

300

Another aspect to consider when using conventional shell elements is the shell contact thickness. For instance,

301

the general contact algorithm from Abaqus/Explicit [28] scales back the shell contact thickness automatically when

302

a certain fraction of the surface facet edge length to thickness is exceeded. This fraction generally varies from 20%

303

to 60% based on the geometry of the element [28]. This is not an issue for the case of Fig. 3a because SFM ensures

304

contact (Fig. 4a). However, the shell contact thickness reduction might lead to interpenetration or inaccurate contacts

305

in the case of deleting cohesive elements (Fig. 3b). This can be avoided by using SFM, as shown in Fig. 4b, to

306

ensure normal contact when the cohesive elements are deleted. Alternatively, Schwab et al. [36] scaled the ply

307

thickness within the contact formulation in regions of delaminations when using conventional shell elements with

308

zero-thickness cohesive elements.

309

[Figure 4 about here.]

310

A study that compares different finite element and cohesive interaction technologies in terms of computational

311

performance and accuracy is missing in the literature. Thus, a benchmark for the different finite element types (i.e.

312

solid, continuum and conventional shell) and cohesive interaction technologies (i.e. zero-thickness cohesive elements

313

and cohesive surfaces) is performed in Section 2.5. 12

314

2.5. LVI and CAI simulation in standard coupons

315

A benchmark for the different finite element types (i.e. solid, continuum shell and conventional shell) and cohesive

316

interaction technologies (i.e. cohesive elements and cohesive surfaces) is carried out for the unidirectional (UD) pre-

317

preg AS4/8552 carbon-epoxy previously studied by Gonz´alez et al. [11] using solid elements with finite thickness

318

cohesive elements. The aforementioned methodology from Section 2.1 - 2.4 is used in the model definitions.

319

The experimental drop weight impact tests were performed according to the ASTM D7136 standard [30] using a

320

CEAST Fractovis Plus instrumented drop-weight tower with an automatic anti-rebound impactor system. CAI tests

321

were performed according to the standard ASTM D7137 [64] on an MTS InsightTM Electromechanical tester with a

322

300 kN load cell.

323

The geometry and boundary conditions of the LVI model are shown in Fig. 5. 16 mm diameter and 5 kg hemi-

324

spherical impactor is modeled with rigid elements (R3D4) with a minimum element size of 0.5 mm. An initial

325

impactor velocity in the out-of-plane direction is assigned according to the impact energy tested. The remaining de-

326

grees of freedom are fixed. The rectangular 100 × 150 mm2 specimen is placed on a flat support with a 125 mm by

327

75 mm rectangular cut-out. The support, which has all the degrees of freedom fixed, is modeled with rigid elements

328

(R3D4) that have a maximum element size of 2.5 mm. The specimen is restrained during the impact by means of four

329

cylinder-shaped clamps with a diameter of 14 mm. They are modeled with rigid elements (R3D4) with a maximum

330

element size of 1 mm and all the degrees of freedom are fixed.

331

To reduce the computational effort, a mesh refined region with regular elements, which ensures correct energy

332

dissipation, is only defined at the impact site. The model proved to be insensitive to any mesh bias effect [11, 65],

333

which is probably related to the small amount of fiber damage. The refined region is 84 × 90 mm2 to allow the damage

334

progression during the LVI and CAI simulation. [Figure 5 about here.]

335

336

The geometry and boundary conditions of the CAI model are sketched in Fig. 6. The specimen fixture assembly

337

was placed between plates and end-loaded under compression until specimen failure. A set of nodes are defined to

338

represent the experimental set-up, which restrains the out-of-plane displacement at the knife edges and the lateral

339

movement at the side clamping zones.

340

[Figure 6 about here.]

341

The case studied is the 19.3 J impact on the [454 /04 / − 454 /904 ]s laminate with a cured thickness of 5.8 mm from

342

[11]. The impact energy falls within the BVID. Ply clustering is modeled as one ply because delamination is not 13

343

expected to occur within ply clusters [11]. Previous studies used conventional shell elements to simulate the LVI and

344

CAI for coupon specimens of thin ply laminates [14] and standard ply laminates [15]. The selected laminate with an

345

unusual ply thickness due to ply clustering is an interesting scenario to assess the possible limitations of conventional

346

shell elements, which are less accurate in the transverse shear stresses and neglect out-of-plane-stresses.

347

The finite element and cohesive interaction types considered in the benchmark are summarized in Table 1. The

348

intralaminar and interlaminar models employed are described in Sections 2.1 and 2.2, respectively. SFM are used for

349

the S4R CON case as described in Gonz´alez et al. [15]. The shell contact algorithm scales the contact thickness due

350

to the used in-plane size (i.e. 0.5 mm) to thickness (i.e. 0.725 mm) ratio. Therefore, SFM are needed in the S4R COH

351

case to guarantee correct contact when the cohesive elements are deleted. [Table 1 about here.]

352

353

The material properties used for the intralaminar damage model are summarized in Table 2 [11]. As in Gonz´alez

354

et al. [11], no plasticity was considered. No experimental data was available to feed the fiber traction separation law

355

shape, thus, the longitudinal tensile traction separation law is selected so as to resemble the linear-exponential damage

356

evolution used from Gonz´alez et al. [11], while the longitudinal compression is a bi-linear traction separation law with

357

a plateau that represents the experimentally observed crushing stress as proposed by other authors [12, 18, 25]. The

358

compression fracture toughness before crushing is G XC = 106.3 N/mm [11], which is the area beneath the first branch

359

of Fig. 2b. [Table 2 about here.]

360

361

362

The in-situ strengths are summarized in Table 3. They are calculated according to [40, 41], Eqs. (1) - (3) and the material properties from Table 2. [Table 3 about here.]

363

364

The material properties used for the interlaminar damage model are summarized in Table 4 [11]. [Table 4 about here.]

365

366

The selected penalty stiffness value is K = Kn = Kcoh = 2.5 ×104 N/mm3 . The effect of the penalty stiffness on the

367

STI is important because it can certainly speed up the simulation time (Eqs. (8) - (9). Thus, it is worth performing a

368

sensitivity analysis beforehand to obtain the optimum value in terms of accuracy (e.g. no effect on system compliance)

369

and computational performance as done in [14]. 14

370

Table 5 shows the maximum element size for each intralaminar damage mechanism based on Eq. (11). Based

371

on the material properties from Table 4 and considering a ply thickness of 0.725 mm, the interlaminar FPZ for pure

372

fracture mode I and mode II are, according to the expressions from [49], 0.73 mm and 2.06 mm, respectively. The

373

actual interlaminar FPZ measured in the LVI and CAI models were approximately 2 mm and 1.5 mm, respectively.

374

The lower interlaminar FPZ during the CAI simulation is due to a less dominant fracture mode II than in LVI. The

375

selected element size for the refined region is 0.5 mm because it ensures no snap-back at the element level for all

376

damage mechanisms and three elements along the interlaminar FPZ (Section 2.3). [Table 5 about here.]

377

378

The force - displacement and the energy absorbed by the plate is compared among the different finite element

379

types and cohesive interaction technologies to the experimental results in Figs. 7 - 8, respectively. Table 6 compares

380

the numerical results with the experimental ones in terms of delamination threshold (Fd ), maximum force (Fmax ),

381

maximum displacement (Umax ), energy dissipated (Edis ), projected delamination area (Ad ) and CAI strength (FCAI ).

382

The elastic regime and delamination threshold is well predicted by all the cases with relative differences below 5%.

383

However, conventional shell elements show a stiffer response when compared to solid and continuum shell elements

384

because they do not account for out-of-plane compliance. The maximum displacement is close to the experimental

385

one for all the cases with differences below 3%. However, larger differences in the prediction of the maximum force

386

are found. All the models tend to over-predict the maximum force by about 10%. Differences among modeling strate-

387

gies are found in the prediction of the dissipated energy and projected delamination area. The predictions are more

388

influenced by the chosen element type than for the chosen cohesive interaction technology. For instance, conventional

389

shell elements under-predict the dissipated energy by about 27.3% while solid elements by 15.8 %. This could be due

390

to the fact that conventional shell elements do not account for out-of-plane-stresses, which is important for thick plies.

391

The continuum shell element suffered from numerical issues which stopped the simulations prematurely regardless

392

of the cohesive interaction technology. Only for large residual matrix stiffness (i.e. Er = (1 − d2 )E2 = 400 MPa) could

393

the simulations be finished.

394

[Figure 7 about here.]

395

[Table 6 about here.]

396

[Figure 8 about here.]

397

The projected delamination area obtained using cohesive elements and cohesive surfaces are shown in Figs. 9 -

398

10, respectively. Conventional shell elements over-predict the projected delamination area by 25% (i.e. S4R CON) 15

399

and 14% (i.e. S4R COH). Solid elements are in better agreement with the experimentally obtained results (i.e. 3%).

400

However, it is worth noting that the case of conventional shell elements with cohesive surface was able to predict the

401

largest delamination which occurred in the 45o direction at the back - face. This could be explained by the fact that

402

conventional shell elements correctly capture bending with only one element per ply and cohesive surfaces are able to

403

better capture large deformation than cohesive elements can.

404

[Figure 9 about here.]

405

[Figure 10 about here.]

406

The CAI experimental test used the displacement measurement from the machine. Thus, the CAI set-up added

407

some compliance, which was also corrected as in [42]. The CAI force - displacement response is shown in Fig. 11.

408

All strategies predict rather accurately the experimental CAI stiffness and strength. CAI strength relative errors of

409

-3.6%, 2.2%, -3.6% and -1.9% were obtained for the cases S4R COH, S4R CON, C3D8R COH and C3D8R CON,

410

respectively. [Figure 11 about here.]

411

412

The computational time for each finite element type and cohesive interaction technology during the LVI is sum-

413

marized in Table 7. The features of the computer used are: two socket work stations (ASUS Z10PE-D16WS Moth-

414

erboard), with two CPU Intel Xeon processors E5-2687Wv3 with 3.1 GHz CPU frequency (10 cores, 160W power,

415

and 2133 DDR4 type memory); solid-state disc of 512 GB and 32 GB RAM for each CPU. All the simulations were

416

computed using 20 CPUs.

417

[Table 7 about here.]

418

The LVI simulation time of conventional shell elements with zero-thickness cohesive elements was 1.54 times

419

faster than when using them with cohesive surfaces. Similarly, solid elements with zero-thickness cohesive elements

420

was 1.3 times faster than when using them with cohesive surfaces. Nevertheless, larger differences in the computa-

421

tional time are found depending on the finite element type. Conventional shell elements were up to 3.3 times faster

422

than solid elements when used with zero-thickness cohesive elements. In fact, the simulation time difference between

423

solid elements and conventional shells should be larger because selective element mass scaling was applied for STI

424

lower than 1 × 10−8 s in order to finish the simulations of solid elements in a reasonable time (see Fig. 12).

425

Table 7 shows that solid and continuum shell elements have a larger initial STI than conventional shell elements

426

do. However, the computational time of conventional shell elements is the lowest. Conventional shell elements were 16

427

also faster than continuum shell elements before the simulation aborted. Conventional shell elements are the fastest

428

because they have a simpler constitutive behaviour and their STI is less affected during the simulation. As shown in

429

Fig. 12, volumetric finite elements suffer greater element distortion than plane elements do.

430

It is worth noting that in the case studies, the element STI was controlled by the in-plane finite element size while

431

for standard ply thickness the STI is usually controlled by the ply thickness. Thus, the computational benefit of using

432

conventional shell elements in situations where the element thickness is controlling the model’s STI is larger. [Figure 12 about here.]

433

434

2.6. Element deletion criterion

435

Excessive element distortion issues, which can eventually abort the simulation, are common in FEA that involve

436

severe damage such as LVI, CAI or crash simulations. This aspect can significantly compromise the accuracy and

437

robustness of the numerical models even though it is often not clearly reported.

438

Some authors [10, 11, 13, 36, 38] delete the continuum elements when a certain value of the fiber damage vari-

439

able is reached, while others [14, 15, 18, 21] limit the damage variable to a certain value without element deletion.

440

Depending on the author, the fiber damage variable can grow up to 0.9 [21], 0.99 [13, 18, 36] or 0.999 [10, 11, 38] in

441

UD laminates. However, this is not enough to avoid element distortion issues in UD laminates because the damaged

442

matrix undergoes severe distortions. In some studies, the transverse and shear damage variables are also limited to

443

grow up to 0.98 [11] or 0.99 [15, 38]. Tan et al. [13] deleted the elements when the shear strain reached 1% or when

444

the fiber damage variable reached 0.99. Lately, Chiu et al. [66] in crash simulations proposed deleting the elements

445

when the fiber damage variable reached 0.99 or when the determinant of the deformation gradient reached certain

446

values, which were based on a sensitivity analysis. Lopes et al. [10] reported for LVI with VID that the elements were

447

deleted when the fiber damage variable reached 0.999 or the matrix damage variable reached 0.99 and 0.99.

448

The damage variable usually depends on the material properties, the traction separation law shape and the element

449

size in the case of using a regularization procedure such as the crack band model [27]. Therefore, special attention

450

must be paid to the damage variable selected for element deletion because the energy dissipated could be far from the

451

actual material fracture toughness. If the element is not deleted but the damage variable is limited to a certain value,

452

the global response could be affected because the remaining element stiffness can absorb significant elastic energy

453

[14]. Furthermore, the residual stresses actually translate into a permanent bridging mechanism behind the crack tip

454

which affects crack progression.

455

In the results from Section 2.5, the continuum elements were deleted when the fiber damage variable was equal

456

to one in order to dissipate the whole fracture energy. The cohesive elements were also deleted when the damage 17

457

variable was equal to one. However, matrix degradation was limited to avoid element distortion issues. Different

458

approaches were considered. Finally, limiting the matrix damage variables, as done in [11, 38], was chosen. The

459

matrix damage variables were limited (i.e. d2 =0.987 , d6 =0.977) so that the matrix never had a residual matrix

460

stiffness (e.g. Er = (1 − d2 )E2 ) lower than 100 MPa. Figs. 13 and 14 show the effect of the residual matrix stiffness

461

value on the LVI and CAI response, respectively.

462

[Figure 13 about here.]

463

[Figure 14 about here.]

464

The results from Figs. 13 - 15 are for the S4R COH case. It can be seen that for a residual matrix stiffness lower

465

than 200 MPa, the LVI and CAI results converge to similar results. Fig. 15 shows the STI evolution during the LVI

466

simulation depending on the matrix residual stiffness. The lower the residual stiffness, the lower the STI during the

467

LVI due to element distortion. In fact, the STI for Er = 500 MPa is fairly constant throughout the LVI simulation

468

because there is little element distortion. [Figure 15 about here.]

469

470

The residual matrix stiffness should be the lowest possible (e.g. Er = 1 MPa). However, C3D8R COH and C3D8R

471

CON from Section 2.5 aborted the simulation for Er < 100 MPa. Thus, Er = 100 MPa was chosen.

472

3. LVI and CAI simulation in an aircraft stiffened panel

473

The most relevant experimental and numerical results from an LVI test and the sequential CAI test of a compos-

474

ite laminate stiffened panel are presented in this section. Conventional shell elements with zero-thickness cohesive

475

elements is the selected modeling strategy (see Fig. 3b), which provides a good balance between accuracy and compu-

476

tational efficiency. In this case, SFM are not needed to ensure correct contact because the in-plane finite element size

477

(i.e. 1 mm) to thickness (i.e. 0.1905 mm) ratio is not scaling the shell contact thickness. Correct contact is guaranteed

478

by the conventional shell elements.

479

3.1. Geometry and boundary conditions

480

The geometry of the composite panel with two stiffeners is shown in Fig. 16. All the plies and interfaces for

481

delamination are modeled in order to avoid affecting the damage sequence due to pre-defined potential interfaces for

482

delamination. Also, all the plies and interfaces are modeled in the stringers because they can play a role in the CAI

483

simulation. As in other studies [11, 13, 38], delamination is not considered in ply clusters. 18

[Figure 16 about here.]

484

485

The clamping system used aluminium blocks bonded to the specimen ends with a Hysol® EA 9394 (epoxy paste)

486

structural adhesive from Henkel® . The aluminium blocks were joined to the high stiffness steel frame of the impact

487

tester using two bolts at each end (see Fig. 17). Each aluminium block covered the whole width and 40 mm in

488

length of the panel (i.e. potting). The specimen was impacted after the aluminium blocks were bonded and cured. No

489

damage before testing was observed through visual inspection, neither at the interface between the aluminium blocks

490

and the specimen nor inside the specimen.

491

492

A 136 J impact was made with a 12.5 kg spherical shaped steel impactor with a diameter of 25.4 mm at the center of the skin, which, in turn, led to VID. [Figure 17 about here.]

493

494

Fig. 18 shows the boundary conditions of the LVI model. Two element sizes were considered in the model: a

495

finer mesh region (i.e. 1 mm) and a coarser mesh region (i.e. 5 mm). The entire panel is meshed with regular finite

496

elements. The refined mesh region is 50 × 240 mm2 to correctly capture the damage resulting from the LVI and CAI.

497

Intralaminar damage is only allowed to occur at the refined mesh region where the finite elements are oriented along

498

the fiber direction to avoid mesh bias effects [65].

499

Neither the structural adhesive nor the aluminium blocks are considered in the model for the sake of simplicity.

500

The boundary conditions are applied through node sets at the panel instead of modeling the whole set-up and taking

501

into account contacts. Fully clamped boundary conditions are only considered at the edges of the panel while the rest

502

of the pottings have all the degrees of freedom fixed with the exception of the longitudinal one in order to account for

503

some compliance added by the adhesive in a simple manner. [Figure 18 about here.]

504

505

The bonded aluminium blocks used for clamping the specimens during the LVI tests were also used for CAI

506

testing. Once the blocks were installed, the specimen was milled to achieve a tolerance of 0.1 mm in parallelism

507

between faces. Additionally, anti-buckling side clamps were applied along the lateral edges to avoid global instability

508

which could trigger compression failure (see Fig. 19). [Figure 19 about here.]

509

510

511

512

The CAI boundary conditions are introduced through node sets that represent the anti-buckling side clamps (Fig. 20). [Figure 20 about here.] 19

513

3.2. Material properties

514

The composite material used is an intermediate modulus carbon/epoxy pre-preg typically used in aerospace appli-

515

cations, with a nominal ply thickness of 0.1905 mm. The continuum damage model described in Section 2.1 is used

516

to model the UD plies. Unfortunately, the experimental test campaign was limited at the sub-component level. Thus,

517

the material properties to feed the models were taken from previous experimental campaigns and similar material

518

systems.

519

Table 8 shows the panel stacking sequence. The skin stacking sequence is according to the reference system from

520

Fig. 16a while the web stacking sequence is sketched in Fig. 21a. The manufacturing process assured the stacking

521

sequence symmetry in the web. The radius and noodle region, which is often a resin rich region, is not explicitly

522

modeled for the sake of simplicity (see Fig. 21b). However, a more detailed model with 3-D elements could be

523

needed in situations there is a significant tension (pull-off) load to the vertical web section [67].

524

[Figure 21 about here.]

525

[Table 8 about here.]

526

527

528

The corresponding UD ply material properties used to feed the intralaminar and interlaminar model are listed in Tables 9 and 10 in a dimensionless form for confidential reasons. [Table 9 about here.]

529

The penalty stiffness and friction coefficient are based on sensitivity analyses. A penalty stiffness value of 5 x 104

530

N/mm3 was enough to avoid adding compliance in the system as shown in Fig. 22. A friction coefficient value of 0.3

531

is used. The model was not sensitive to friction coefficients larger than 0.1.

532

[Table 10 about here.]

533

[Figure 22 about here.]

534

535

536

The in-situ properties are calculated according to [40, 41] and Eqs. (1) - (3). Fig. 23 shows a graphical representation of the in-situ properties of the material used. [Figure 23 about here.]

537

The element deletion criteria used is the same described in Section 2.6. The elements are deleted when the fiber

538

damage variable reaches a value equal to one, while matrix degradation is limited (i.e. d2 , d6 ) to leave a residual 20

539

matrix stiffness of 100 MPa in order to avoid element distortion issues. The cohesive elements are also deleted when

540

the damage variable is equal to one.

541

Table 11 summarizes the maximum intralaminar element size for each damage mechanism according to Eq. (11).

542

Using the expressions from [49] and considering the smallest ply thickness modeled (i.e. 0.1905 mm), the interlaminar

543

FPZ are 0.25 and 2.38 mm for pure mode I and mode II, respectively. However, the actual interlaminar FPZ was

544

measured in preliminary simulations in order to use a less conservative element size, which was about 3 mm. An

545

element size of 1 mm ensured correct intralaminar and interlaminar energy dissipation. [Table 11 about here.]

546

547

3.3. Low velocity impact results

548

The force - displacement response from the numerical model is compared with the experimental one in Fig. 24.

549

The overall response is well predicted, taking into account the complexity of the test and the large impact energy

550

which lead to a large amount of damage and skin perforation. [Figure 24 about here.]

551

552

The elastic response (point A) and damage initiation (point B) are fairly well captured. The maximum force is

553

reasonably well predicted but the numerical model predicts skin perforation at a much lower displacement (Point

554

C). The model shows a more brittle failure than the test. It is worth noting that only one panel was tested so the

555

possible experimental scatter is not known. The numerical model predicts a load drop that is mainly due to fiber

556

breakage rather than delamination. Afterwards, there is a stage of damage growth in which the model significantly

557

over-predicts the force. Fig. 25 shows that the predicted force is significantly improved when the residual matrix

558

stiffness (Er ) is reduced.

559

Initially, Er = 100 MPa was taken. The residual matrix stiffness should be the lowest possible but element distor-

560

tion issues had already aborted the simulation for Er < 50 MPa. Finally, a residual matrix stiffness of 50 MPa was

561

used which allowed the LVI simulation to be finished and the sequential CAI simulation to be performed, albeit Fig.

562

25 highlighting the effect on the results.

563

[Figure 25 about here.]

564

LVI simulations with VID [21] are more complex than simulations with BVID [18, 19, 22] because they suffer all

565

forms of damage (e.g. matrix cracking, delamination and fiber breakage), especially in the form of fiber breakage. In

566

LVI simulations with VID becomes very important to capture the interaction among damage mechanisms for accurate 21

567

results. For instance, the fiber traction separation law shape can play an important role. Unfortunately, the sub-

568

component experimental test campaign was not accompanied with a material characterization which makes it difficult

569

to identify the causes of disagreement between the model and the experimental test. Furthermore, LVI simulations

570

with VID are more complex from a numerical point of view because they suffer severe element distortion that can

571

abort the analysis. The model’s accuracy and robustness can be sensitive to the element deletion criterion. Soto et

572

al. [14] reported sensitivity of the results even for very low residual stiffness values when there is extensive damage

573

(VID).

574

The element deletion criterion represents a limitation of the approach presented here when applied to the simula-

575

tion of LVI and sequential CAI with VID, although there was no issue when applied for BVID in Section 2. It seems

576

an intrinsic issue of using a smeared damage approach. Some authors [7, 20] use cohesive elements to model matrix

577

cracking. This approach could circumvent numerical issues related to element distortion even though other aspects

578

arise (e.g. mesh bias effects, discretization criterion and computational time). [Figure 26 about here.]

579

580

Fig. 26 shows the experimental projected delamination area of the panel (Ad = 1575 mm2 ) together with the

581

numerical one (Ad =1236 mm2 ). The projected delamination area is under-predicted by 21.5 %. The model predicts

582

that the largest delamination appears at the back - face. Probably, the largest delamination in the test also occurred at

583

the back-face. However, the model is not able to capture the the extension of it. This could explain the under-prediction

584

of the delamination area.

585

Section 2 showed that it is difficult to capture the back - face delamination (see Figs. 9 - 10). Only the case of

586

conventional shell elements with a cohesive surface was able to predict the largest delamination, which occurred at

587

the back - face in the 45o direction. Thus, the under-prediction of the stripes and delamination could be due to the

588

modeling strategy that affects the interaction between interlaminar and intralaminar damage mechanisms. However,

589

the fiber traction separation law, which was not characterized, can also affect the interaction between interlaminar and

590

intralaminar damage but the discrepancies seem to be large to only attribute them to the modeling approach. The

591

material properties should be tested, especially those properties related to fracture propagation.

592

Fig. 27 shows the contour area of deleted elements due to fiber damage with the perforated skin region. The

593

perforated area is smaller than the area with deleted elements due to fiber damage. Thus, the fiber damage growth

594

seems to be over-predicted by the model while delamination was under-predicted.

595

[Figure 27 about here.]

22

596

597

Fig. 28 shows a view of the LVI simulation with the fiber damage field output where it can be seen skin perforation through the composite panel. [Figure 28 about here.]

598

599

The machine used to calculate the LVI of the composite panel was the same as for the numerical benchmark from

600

Section 2.5. The model had 11.45 million degrees of freedom and the computational time using 20 CPUs was 69

601

hours.

602

3.4. Compression after impact results

603

In this section, the CAI test is shown together with a pristine test that was performed on a panel with the same

604

configurations. Fig. 29 compares the numerical results with the experimental ones from both pristine and CAI tests.

605

The displacements were measured through a LVDT. Thus, no compliance calibration is needed as in Section 2.5 in

606

which the displacements were measured from the machine. [Figure 29 about here.]

607

608

The pristine panel failed at the bottom middle section as shown in Fig. 30. The pristine model had the same model

609

definitions as described for the CAI model but the mesh refined region was located at the bottom middle section to

610

capture the experimentally observed failure location. In the CAI model, the impactor was first removed to avoid

611

contact with the panel. [Figure 30 about here.]

612

613

The pristine simulation captures the stiffness up to F/Fmax = 0.6. Afterwards, the pristine test shows a loss of

614

stiffness. The model matches the stiffness again around F/Fmax = 0.9 when it buckles. The presence of some defect

615

in the panel could affect the buckling prediction that the model does not capture. The pristine model under-predicts

616

the compression panel strength by 6.2 %.

617

On the other hand, the CAI test failed successfully at the impacted region. Fig. 31 shows the experimental

618

and numerical failed CAI specimens, respectively. Fiber failure can be observed across the stringers width of the

619

numerical model as well as in the tested panel. A fracture that crosses the whole width is not seen in the numerical

620

model. Nevertheless, it is attributed to the fact that the model eventually aborts for excessive element distortion prior

621

to complete failure.

622

[Figure 31 about here.] 23

623

The CAI numerical model is able to predict the stiffness even the damage progression prior to failure with

624

fairly good agreement. The numerical model under-predicts the CAI strength by 5.36 % in the first load drop (i.e.

625

u/umax =0.6). Both pristine and CAI strength under-predictions are of similar magnitude what it could be attributed

626

to some discrepancy on the used fiber compression properties. Despite the discussed differences between the LVI

627

numerical model and the test results, the CAI strength is reasonably well predicted (relative error of 5.36 %). The

628

CAI model takes into account significant damage occurred at the impacted side during the impact. This is probably an

629

important aspect for the panel CAI strength prediction. The final failure is controlled by compression fiber breakage

630

at the stringers since the skin is severely damaged.

631

4. Conclusions

632

An efficient methodology to perform LVI and CAI in composite laminate sub-components was defined. Numerical

633

impact models require a number of several parameters that affect their efficiency, accuracy, robustness and objectivity

634

to be defined. Here, the key parameters have been described concisely, the selected values have been well justified

635

and their sensitivity discussed.

636

Firstly, the methodology was applied at the coupon level with the aim of performing a numerical benchmark that

637

systematically studies the computational performance and numerical predictions of different finite element types and

638

cohesive interaction technologies. The selected impact energy led to BVID. A laminate with unusually thick plies was

639

selected to highlight possible limitations of plane-stress based finite elements. The use of conventional shell elements

640

together with cohesive elements proved to be the fastest approach with relative errors on the CAI strength below

641

4%. However, conventional shell elements were less accurate than solid elements on the LVI prediction because of

642

poorer transverse shear description and neglect of out-of-plane-stress, both of which are more important for damage

643

predictions in thick ply laminates than in standard ply laminates.

644

Secondly, the methodology was applied for the LVI and CAI simulation of a composite stiffened panel. The se-

645

lected impact energy led to VID. Conventional shell elements together with zero-thickness cohesive elements were

646

used in order to make feasibly the simulation of the sub-component LVI and CAI. The impact prediction was not in

647

agreement with the experimental test in terms of projected delamination area and force-displacement response after

648

skin perforation. The material properties should be tested in order to identify possible limitations of the employed

649

approach. However, the predicted damage was good enough to obtain a reasonably good CAI prediction. The relative

650

error in the compression strength prediction of the undamaged and damaged panel were -6.2 % and -5.36 %, respec-

651

tively. The present work showed the potential the current numerical tools and material models to address problems of

652

industrial interest such as the strength prediction of both undamaged and damaged stiffened panels. The structure size 24

653

considered together with the studied impact energy, which led to skin perforation, represented a challenging case.

654

Nevertheless, the element deletion criterion proved to be the most important limitation of the presented approach

655

when there is VID. Further work is required to improve the FEA robustness of LVI and CAI simulations when there

656

is extensive damage (VID). It is thought that the use of discrete approaches (e.g. smoothed particle hydrodynamics)

657

could circumvent the reported limitations related to element distortion.

658

Acknowledgments

659

The authors gratefully acknowledge the support provided by Embraer in the form of all the experimental data of

660

stiffened panels, which were tested at TEAMS (Element facilities in Seville). The first author would like to acknowl-

661

edge the predoctoral grant BES-2014-070374 from the Spanish Governement. This work has been partially funded by

662

the Spanish Government through the Ministerio de Econom´ıa y Competitividad, under contracts MAT2013-46749-R

663

and MAT2015-69491-C3-1-R. The last author thanks FCT-Funda£o para a Ciłncia e a Tecnologia for the support

664

provided through the project MITPTB/PFM/0005/2013.

665

666

References

667

[1] Abrate S. Modeling of impacts on composite structures. Composite Structures 2001;51(2):129–138.

668

[2] Olsson R. Analytical prediction of large mass impact damage in composite laminates. Composites Part A: Applied Science and Manufacturing

669 670 671

2001;32(9):1207–1215. [3] Olsson R, Donadon MV, Falzon BG. Delamination threshold load for dynamic impact on plates. International Journal of Solids and Structures 2006;43(10):3124–3141.

672

[4] Naik NK, Ramasimha R. Estimation of compressive strength of delaminated composites. Composite Structures 2001;52(2):199–204.

673

[5] Hou JP, Petrinic N, Ruiz C, Hallett SR.

674 675 676 677 678 679 680 681 682 683 684

Prediction of impact damage in composite plates.

Composites Science and Technology

2000;60(2):273–281. [6] Lammerant L, Verpoest I. Modelling of the interaction between matrix cracks and delaminations during impact of composite plates. Composite Science and Technology 1996;(22):435–447. [7] Bouvet C, Castani´e B, Bizeul M, Barrau JJ. Low velocity impact modelling in laminate composite panels with discrete interface elements. International Journal of Solids and Structures 2009;46(14):2809–2821. [8] Raimondo L, Iannucci L, Robinson P, Curtis PT. A progressive failure model for mesh-size-independent FE analysis of composite laminates subject to low-velocity impact damage. Composites Science and Technology 2012;72(5):624–632. [9] Feng D, Aymerich F. Finite element modelling of damage induced by low-velocity impact on composite laminates. Composite Structures 2014;108(Supplement C):161–171. [10] Lopes CS, S´adaba S, Gonz´alez C, Llorca J, Camanho PP. Physically-sound simulation of low-velocity impact on fiber reinforced laminates. International Journal of Impact Engineering 2016;92(Supplement C):3–17.

25

685 686 687 688 689 690 691 692 693 694

[11] Gonz´alez EV, Maim´ı P, Camanho PP, Turon A, Mayugo JA. Simulation of drop-weight impact and compression after impact tests on composite laminates. Composite Structures 2012;94(11):3364–3378. [12] Rivallant S, Bouvet C, Hongkarnjanakul N. Failure analysis of CFRP laminates subjected to compression after impact: FE simulation using discrete interface elements. Composites Part A: Applied Science and Manufacturing 2013;55(Supplement C):83–93. [13] Tan W, Falzon BG, Chiu LNS, Price M. Predicting low velocity impact damage and Compression-After-Impact (CAI) behaviour of composite laminates. Composites Part A: Applied Science and Manufacturing 2015;71(Supplement C):212–226. [14] Soto A, Gonz´alez E, Maim´ı P, de la Escalera FM, de Aja JS, Alvarez E. Low velocity impact and compression after impact simulation of thin ply laminates. Composites Part A: Applied Science and Manufacturing 2018;109:413 – 427. [15] Gonz´alez E, Maim´ı P, Mart´ın-Santos E, Soto A, Cruz P, de la Escalera FM, de Aja JS. Simulating drop-weight impact and compression after impact tests on composite laminates using conventional shell finite elements. International Journal of Solids and Structures 2018;.

695

[16] Krajcinovic D. Continuum damage mechanics. Applied Mechanics Review 1984;37:1–6.

696

[17] Ortiz M. A constitutive theory for the inelastic behavior of concrete. Mechanics of Materials 1985;4(1):67–93.

697

[18] Faggiani A, Falzon BG. Predicting low-velocity impact damage on a stiffened composite panel. Composites Part A: Applied Science and

698 699 700 701 702 703 704

Manufacturing 2010;41(6):737–749. [19] Riccio A, Ricchiuto R, Saputo S, Raimondo A, Caputo F, Antonucci V, Lopresto V. Impact behaviour of omega stiffened composite panels. Progress in Aerospace Sciences 2016;81(Supplement C):41–48. [20] Sun XC, Hallett SR. Barely visible impact damage in scaled composite laminates: Experiments and numerical simulations. International Journal of Impact Engineering 2017;109(Supplement C):178–195. [21] Johnson HE, Louca LA, Mouring S, Fallah AS. Modelling impact damage in marine composite panels. International Journal of Impact Engineering 2009;36(1):25–39.

705

[22] Psarras S, Mu˜noz R, Ghajari M, Robinson P, Furfari D, Hartwig A, Newman B. Compression after multiple impacts: modelling and exper-

706

imental validation on composite curved stiffened panels. Smart Intelligent Aircraft Structures (SARISTU): Proceedings of the Final Project

707 708 709 710 711 712 713 714 715

Conference. Springer International Publishing; 2016, p. 681–689. [23] Williams KV, Vaziri R, Poursartip A. A physically based continuum damage mechanics model for thin laminated composite structures. International Journal of Solids and Structures 2003;40(9):2267–2300. [24] Johnson AF, Pickett AK, Rozycki P. Computational methods for predicting impact damage in composite structures. Composites Science and Technology 2001;61(15):2183–2192. [25] Donadon MV, Iannucci L, Falzon BG, Hodgkinson JM, de Almeida SFM. A progressive failure model for composite laminates subjected to low velocity impact damage. Computers & Structures 2008;86(11):1232–1252. [26] Camanho PP, Davila CG, de Moura MF. Numerical Simulation of Mixed-Mode Progressive Delamination in Composite Materials. Journal of Composite Materials 2003;37(16):1415–1438.

716

[27] Baˇzant ZP, Oh BH. Crack band theory for fracture of concrete. Mat´eriaux et Construction 1983;16(3):155–177.

717

[28] Abaqus Analysis User’s Manual. In: Abaqus documentation 6.11. Providence, RI 02909 2499, USA: Simulia World Headquarters; 2011,.

718

[29] Iannucci L, Willows ML. An energy based damage mechanics approach to modelling impact onto woven composite materials Part I: Numer-

719 720 721 722 723

ical models. Composites Part A: Applied Science and Manufacturing 2006;37(11):2041–2056. [30] ASTM D 7136/D 7136M-12. Standard test method for measuring the damage resistance of a fiber-reinforced polymer matrix composite to a drop-weight impact event. Standard; ASTM International; West Conshohocken PA, USA; 2012. [31] Ladev`eze P, LeDantec E. Damage modelling of the elementary ply for laminated composites. Composites Science and Technology 1992;43(3):257–267.

26

724 725 726 727 728 729 730 731 732 733

[32] Matzenmiller A, Lubliner J, Taylor RL.

A constitutive model for anisotropic damage in fiber-composites.

Mechanics of Materials

1995;20(2):125–152. [33] Maim´ı P, Camanho PP, Mayugo JA, D´avila CG. A continuum damage model for composite laminates: Part I Constitutive model. Mechanics of Materials 2007;39(10):897–908. [34] Maim´ı P, Camanho PP, Mayugo JA, D´avila CG. A continuum damage model for composite laminates: Part II Computational implementation and validation. Mechanics of Materials 2007;39(10):909–919. [35] Mart´ın-Santos E, Maim´ı P, Gonz´alez EV, Cruz P. A continuum constitutive model for the simulation of fabric-reinforced composites. Composite Structures 2014;111(Supplement C):122–129. [36] Schwab M, Todt M, Wolfahrt M, Pettermann HE. Failure mechanism based modelling of impact on fabric reinforced composite laminates based on shell elements. Composites Science and Technology 2016;128(Supplement C):131–137.

734

[37] D´avila CG, Camanho PP, Rose CA. Failure Criteria for FRP Laminates. Journal of Composite Materials 2005;39(4):323–345.

735

[38] Lopes CS, Camanho PP, G¨urdal Z, Maim´ı P, Gonz´alez EV. Low-velocity impact damage on dispersed stacking sequence laminates. Part II:

736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762

Numerical simulations. Composites Science and Technology 2009;69(7):937–947. [39] Catalanotti G, Camanho P, Marques A. Three-dimensional failure criteria for fiber-reinforced laminates. Composite Structures 2013;95:63 – 79. [40] Camanho PP, D´avila CG, Pinho ST, Iannucci L, Robinson P. Prediction of in situ strengths and matrix cracking in composites under transverse tension and in-plane shear. Composites Part A: Applied Science and Manufacturing 2006;37(2):165–176. [41] Maim´ı P, Camanho PP, Mayugo JA, Turon A. Matrix cracking and delamination in laminated composites. Part I: Ply constitutive law, first ply failure and onset of delamination. Mechanics of Materials 2011;43(4):169–185. [42] D´avila CG, Rose CA, Camanho PP. A procedure for superposing linear cohesive laws to represent multiple damage mechanisms in the fracture of composites. International Journal of Fracture 2009;158(2):211–223. [43] Catalanotti G, Arteiro A, Hayati M, Camanho PP. Determination of the mode I crack resistance curve of polymer composites using the size-effect law. Engineering Fracture Mechanics 2014;118(Supplement C):49–65. [44] Ortega A, Maim´ı P, Gonz´alez EV, Trias D. Characterization of the translaminar fracture Cohesive Law. Composites Part A: Applied Science and Manufacturing 2016;91(Part 2):501–509. [45] Orifici AC, Herszberg I, Thomson RS. Review of methodologies for composite material modelling incorporating failure. Composite Structures 2008;86(1):194–210. [46] Turon A, Camanho PP, Costa J, D´avila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mechanics of Materials 2006;38(11):1072–1089. [47] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites Science and Technology 1996;56(4):439–449. [48] Alfano G. On the influence of the shape of the interface law on the application of cohesive-zone models. Composites Science and Technology 2006;66(6):723–730. [49] Soto A, Gonz´alez EV, Maim´ı P, Turon A, de Aja JRS, de la Escalera FM. Cohesive zone length of orthotropic materials undergoing delamination. Engineering Fracture Mechanics 2016;159(Supplement C):174–188. [50] Dourado N, de Moura M, de Morais AB, Pereira AB. Bilinear approximations to the mode II delamination cohesive law using an inverse method. Mechanics of Materials 2012;49:42–50. [51] Catalanotti G, Furtado C, Scalici T, Pitarresi G, van der Meer FP, Camanho PP. The effect of through-thickness compressive stress on mode II interlaminar fracture toughness. Composite Structures 2017;182:153–163.

27

763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780

[52] Daudeville L, Allix O, Ladev`eze P. Delamination analysis by damage mechanics: Some applications. Composites Engineering 1995;5(1):17– 24. [53] Turon A, D´avila C, Camanho P, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engineering Fracture Mechanics 2007;74(10):1665–1682. [54] Iannucci L.

Progressive failure modelling of woven carbon composite under impact.

International Journal of Impact Engineering

2006;32(6):1013–1043. [55] Turon A, Gonz´alez EV, Sarrado C, Guillamet G, Maim´ı P. Accurate simulation of delamination under mixed-mode loading using a cohesive model with a mode-dependent penalty stiffness. Composite Structures 2018;184(Supplement C):506–511. [56] Maim´ı P, Gonz´alez EV, Gascons N, Ripoll L. Size Effect Law and Critical Distance Theories to Predict the Nominal Strength of Quasibrittle Structures. Applied Mechanics Reviews 2013;65(2):020803. [57] Kabeel AM, Maim´ı P, Gascons N, Gonz´alez EV. Nominal strength of quasi-brittle open hole specimens under biaxial loading conditions. Composites Science and Technology 2013;87:42–49. [58] Mi Y, Crisfield MA, Davies GAO, Hellweg HB. Progressive Delamination Using Interface Elements. Journal of Composite Materials 1998;32(14):1246–1272. [59] Falk ML, Needleman A, Rice JR. A critical evaluation of cohesive zone models of dynamic fracture. In: Proceedings of the 5th European mechanics of materials conference on scale transitions from atomistics to continuum plasticity. 2001, p. 43–50. [60] Tabiei A, Tanov R. A nonlinear higher order shear deformation shell element for dynamic explicit analysis:: Part I. Formulation and finite element equations. Finite Elements in Analysis and Design 2000;36(1):17–37.

781

[61] Noor AK, Burton WS, Bert CW. Computational Models for Sandwich Panels and Shells. Applied Mechanics Reviews 1996;49(3):155–199.

782

[62] Borg R, Nilsson L, Simonsson K. Simulating dcb, enf and mmb experiments using shell elements and a cohesive zone model. Composites

783 784 785 786 787

Science and Technology 2004;64(2):269 – 278. [63] D´avila CG, Camanho PMPRd, Turon Travesa A, et al. Cohesive elements for shells. © NASA TP Technical Reports, 2007, n´um 214869 2007;. [64] ASTM D 7137/D 7137M-12. Standard test method for compressive residual strength properties of damaged polymer matrix composite plates. Standard; ASTM International; West Conshohocken PA, USA; 2012.

788

[65] Laˇs V, Zemˇc´ık R. Progressive Damage of Unidirectional Composite Panels. Journal of Composite Materials 2008;42(1):25–44.

789

[66] Chiu LNS, Falzon BG, Boman R, Chen B, Yan W. Finite element modelling of composite structures under crushing load. Composite

790 791 792

Structures 2015;131(Supplement C):215–228. [67] Trask R, Hallett S, Helenon F, Wisnom M. Influence of process induced defects on the failure of composite t-joint specimens. Composites Part A: Applied Science and Manufacturing 2012;43(4):748 – 757.

28

σ

12

S S

L

KG 1+K P

LP

G

12

12 P

G l*

SL

g

12

Figure 1: In-plane shear stress - strain constitutive response with uncoupled linear hardening and linear softening damage. The parameters G12 , S L , S LP , KP , GS L and l∗ are the shear elastic modulus, the shear strength, the yield stress, the hardening plastic parameter, the shear fracture toughness and the characteristic element length, respectively.

29

σ

1

1

XT

2

M

C

XT

fX X T

T

GXT (a)

fX X C

GYM C

δ

δ

δ (b)

(c)

Figure 2: (a) Longitudinal tensile, (b) longitudinal compression and (c) transverse traction separation law shape where M = T, C. Stress vs. opening displacement.

30

S4R SFM SFM S4R

S4R

Tie

Tie COH Tie

Tie S4R

(b)

(a)

Figure 3: (a) Modeling strategy with conventional shell elements (S4R) linked with surface elements (SFM) or (b) zero-thickness cohesive elements (COH) through tie constraints for correct kinematic description of delamination. (The reader is referred to the web version of this article for interpretation of the colors).

31

S4R SFM SFM S4R

S4R Tie SFM SFM Tie S4R

Tie COH

Tie

contact thickness

contact thickness

(a)

(b)

Figure 4: (a) Sketch of modeling strategy with conventional shell elements (S4R) linked with surface elements (SFM) or (b) zero-thickness cohesive elements (COH) through tie constraints when shell contact thickness reduction occurs. (The reader is referred to the web version of this article for interpretation of the colors).

32

Figure 5: Sketch of the LVI numerical model boundary conditions. (The reader is referred to the web version of this article for interpretation of the color legend).

33

Displacement BC Fixture bottom BC Fixture top BC Fixture right BC Knife edge BC Coarse mesh Fine mesh

Figure 6: Sketch of the CAI numerical model boundary conditions. (The reader is referred to the web version of this article for interpretation of the color legend).

34

9 8 7

Force [kN]

6 5 4 3 Experimental S4R COH S4R CON C3D8R COH C3D8R CON SC8R COH SC8R CON

2 1 0 0

0.5

1

1.5

2

2.5

Displacement [mm]

3

3.5

4

Figure 7: Comparison of the force - displacement response for different finite element types and interaction technologies. (The reader is referred to the web version of this article for interpretation of the color legend).

35

20

Energy [J]

16

12

8

Experimental S4R COH S4R CON C3D8R COH C3D8R CON

4

0 0

0.5

1

1.5

2

2.5

3

Time [ms]

3.5

4

4.5

5

Figure 8: Comparison of the energy absorbed among different finite element types and interaction technologies. (The reader is referred to the web version of this article for interpretation of the color legend).

36

4016.54 mm2

(b)

(a)

Figure 9: (a) Projected delamination area of conventional shells and (b) solid elements with zero-thickness cohesive elements. The experimental projected delamination area is the dashed line.

37

3924.71 mm2

4885.19 mm2

(b)

(a)

Figure 10: (a) Projected delamination area of conventional shells and (b) solid elements with cohesive contact surfaces. The experimental projected delamination area is the dashed line

38

110 100 90 80

Force [kN]

70 60 50 40 30 Experimental S4R COH S4R CON C3D8R COH C3D8R CON

20 10 0 0

0.2

0.4

Displacement [mm]

0.6

0.8

Figure 11: Comparison of the force - displacement CAI results for different finite element types and interaction technologies. (The reader is referred to the web version of this article for interpretation of the color legend).

39

−5

x 10

S4R COH S4R CON C3D8R COH C3D8R CON SC8R COH SC8R CON

Stable time increment [ms]

6

5

4

3

2

1 0

1

2

3

Time [ms]

4

5

Figure 12: Comparison of the STI evolution during the LVI between different element and interaction technologies. (The reader is referred to the web version of this article for interpretation of the color legend).

40

9 8 7

Force [kN]

6 5 4 3 Experimental Er=1 MPa

2

Er=100 MPa Er=200 MPa

1

Er=500 MPa

0

0

0.5

1

1.5

2

2.5

3

3.5

4

Displacement [mm]

Figure 13: Effect of the residual matrix stiffness on the LVI force - displacement response. S4R COH case. (The reader is referred to the web version of this article for interpretation of the color legend).

41

110 100 90 80

Force [kN]

70 60 50 40 30

Experimental Er=1 MPa

20

Er=100 MPa

10

Er=200 MPa

0

Er=500 MPa

0

0.2

0.4

0.6

0.8

Displacement [mm]

Figure 14: Effect of the residual matrix stiffness on the CAI force - displacement response. S4R COH case. (The reader is referred to the web version of this article for interpretation of the color legend).

42

6.5

×10-5 Er=1 MPa Er=100 MPa Er=200 MPa

Stable time increment [ms]

5.5

Er=500 MPa

4.5

3.5

2.5

1.5

0.5 0

1

2

3

4

5

Time [ms]

Figure 15: Effect of the residual matrix stiffness on the STI evolution during the LVI simulation. S4R COH case. (The reader is referred to the web version of this article for interpretation of the color legend).

43

stringer web 3.05 mm

stringer base skin

m

0m 50

25 0

mm

150 mm 250 mm

(a) (b) Figure 16: (a) Dimensions of the stiffened panel and (b) its section.

44

Figure 17: Sketch of the clamping system.

45

Impactor BC Clamping BC Coarse mesh Fine mesh

Figure 18: Panel model assembly and boundary conditions for the LVI. (The reader is referred to the web version of this article for interpretation of the color legend).

46

Figure 19: View of the CAI test with the anti-buckling clamps highlighted in red. (The reader is referred to the web version of this article for interpretation of the colors).

47

Fixture top BC Fixture right BC Knife edge BC Coarse mesh Fine mesh

Figure 20: Sketch of the CAI model and boundary conditions. (The reader is referred to the web version of this article for interpretation of the color legend).

48

Ply orientations

R = 4 mm 33.5 mm

45º 64 mm Filler plies

(a)

(b)

Figure 21: a) Sketch of stringer lay-up and b) model cross-section. (The reader is referred to the web version of this article for interpretation of the color legend).

49

Experimental Numerical. K=1e5 Numerical. K=5e4

1

F/ Fmax

0.75

0.5

0.25

0

0

0.25 u/umax

0.5

Figure 22: Penalty stiffness effect (K) on the elastic force - displacement response of the panel.

50

Yis T emb Yis T out YUD T Sis L emb

In-situ strength

Sis L out SUD L

0

0.25

0.5

0.75

1

1.25

1.5

Ply thickness [mm]

Figure 23: Graphical representation of the tensile and shear in-situ strength of the material used in Section 3. (The reader is referred to the web version of this article for interpretation of the colors).

51

Experimental Numerical

1

C

F/ Fmax

0.75 B D

0.5

A

0.25

0

0

0.25

0.5 u/umax

0.75

1

Figure 24: Comparison of the numerical force - displacement impact response with the experimental one.

52

Experimental Numerical. Er=100MPa

1

Numerical. Er=50MPa

F/ Fmax

0.75

0.5

0.25

0

0

0.25

0.5 u/umax

0.75

1

Figure 25: Effect of the residual matrix stiffness (Er ) on the LVI simulation.

53

Figure 26: Zoomed view from the impacted side after the LVI. The white dashed area is the experimental projected delamination area (Ad = 1575 mm2 ) while the dark dashed line is the contour of the numerical projected delamination area (Ad = 1236 mm2 ).

54

Figure 27: Zoomed view from the impacted side. The dark line represents the contour area of deleted elements due to fiber breakage (d1 =1).

55

Figure 28: View of skin perforation at instant u/umax = 0.8. Field output of the fiber damage variable (SDV8).

56

Pristine test Pristine simulation CAI test CAI simulation

1

F/Fmax

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6 u/umax

0.8

1

Figure 29: Comparison between experimental and numerical results of the force - displacement response of the pristine and CAI tests. (The reader is referred to the web version of this article for interpretation of the color legend).

57

Figure 30: Front view of failed pristine specimen.

58

Figure 31: a) Back - face image of the experimentally tested CAI panel and b) the numerical model.

59

Case SC8R CON SC8R COH C3D8R CON C3D8R COH S4R CON S4R COH

FE type Continuum shell (SC8R) Continuum shell (SC8R) Solid (C3D8R) Solid (C3D8R) Conventioanl shell (S4R) Conventioanl shell (S4R)

Interaction type Cohesive surfaces Zero-thickness cohesive elements (COH3D8) Cohesive surfaces Zero-thickness cohesive elements (COH3D8) Cohesive surfaces Zero-thickness cohesive elements (COH3D8)

Table 1: Finite element and interaction types studied in the numerical benchmark.

60

Density Elastic properties Strength

Fracture toughness

Traction separation law parameters

1.59 × 10−9 t/mm3 E1 = 128000 MPa ; E2 = E3 = 7630 MPa G12 = 4358 MPa ; ν12 = 0.35 ; ν23 = 0.45 XT = 2300 MPa ; XC = 1531 MPa YT = 26 MPa ; YC = 199.8 MPa S L = 78.4 MPa G XT = 81.5 N/mm ; G XC = 106.3 N/mm GYT = 0.28 N/mm ; GYC = 1.31 N/mm GS L = 0.79 N/mm fXT = 0.1 ; fXC = 0.1 HXT = 48681 N/mm3 ; HXC = 9922.7 N/mm3

Table 2: Material properties and model parameters of Hexply AS4/8552 used [11] for the intralaminar damage model.

61

Ply location Outer plies Embedded plies Embedded (symmetry) plies

h ply [mm] 0.725 0.725 1.45

YTis [MPa] 38.75 61.23 43.32

YCis [MPa] 199.8 199.8 199.8

Table 3: In-situ strength for the intralaminar damage model.

62

S isL [MPa] 78.4 78.4 78.4

G Ic [N/mm] 0.28

G IIc [N/mm] 0.79

τI [MPa] 26

τII [MPa] 78.4

η 1.45

Table 4: Interface material properties of Hexply AS4/8552 [11].

63

Damage mechanism XT XC YT YC SL

Maximum element size [mm] 2.37 11.6 1.14 0.5 1.12

Table 5: Maximum intralaminar element size for each damage mechanism.

64

Case Experimental S4R CON S4R COH C3D8R CON C3D8R COH SC8R CON SC8R COH

Fd [kN] 4.41 4.41 4.46 4.2 4.2 4.15 4.25

Fmax [kN] 7.74 8.66 8.65 8.65 8.34 -

Umax [mm] 3.72 3.66 3.64 3.81 3.81 -

Edis [J] 12.03 8.58 8.87 10.23 9.99 -

Ad [mm2 ] 3898.3 4885.19 4455.36 3924.71 4016.54 -

FCAI [kN] 105.32 107.7 101.54 103.29 101.51 -

Table 6: Comparison of experimental and numerical results in terms of delamination threshold (Fd ), maximum force (Fmax ), maximum displacement (Umax ), energy dissipated (Edis ), projected delamination area (Ad ) and CAI strength (FCAI ) for the studied cases.

65

Case S4R COH S4R CON C3D8R COH C3D8R CON SC8R COH SC8R CON

Initial STI [s] 4.51 10−8 4.75 10−8 4.73 10−8 5.01 10−8 4.73 10−8 4.99 10−8

LVI CPU time [h] 13 20 43 56 Not finished Not finished

CAI CPU time [h] 5.5 13 23 28 Not finished Not finished

Table 7: Comparison of the initial STI and computational time required for the LVI and CAI simulation for each modeling strategy.

66

Part Skin Web Base right Base left

Stacking sequence [45/90/ − 45/0/45/90/ − 45/0]s [45/0/ − 45/90/45/90/ − 45/0/0/ − 45/90/45]s [45/0/ − 45/90/45/90/ − 45/0/0/ − 45/90/45/90/ − 45/0/45] [45/0/ − 45/90/ − 45/90/45/0/0/45/90/ − 45/90/45/0/ − 45]

Table 8: Panel stacking sequence.

67

Density Elastic properties Strength

Fracture toughness

Traction separation law parameters

1.59 × 10−9 t/mm3 E1 /E2 = 15.78 ; E2 = E3 G12 /E2 = 0.478 ; ν12 /ν23 = 0.88 XT /XC = 1.6 YT /YC = 0.2 YT /S L = 0.59 G XT /G XC = 1.33 GYT /GYC = 0.1 GYT /GS L = 0.17 fXT = 0.4 ; fXC = 0.2 HXT /HXC = 1.6

Table 9: Dimensionless material properties and parameters used for the intralaminar damage model.

68

G Ic /G IIc 0.17

τIc /τIIc 0.59

η 1.9

K [N/mm3 ] 5 × 104

Table 10: Dimensionless material properties and parameters used for the interlaminar damage model.

69

Damage mechanisms XT XC YT YC SL

Maximum element size [mm] 10.9 17.2 0.94 0.95 1.12

Table 11: Maximum intralaminar element size for each damage mechanism to avoid snap-back at element level.

70