Meso-scale modeling of hypervelocity impact damage in composite laminates

Composites Part B 74 (2015) 95e103

Contents lists available at ScienceDirect

Composites Part B journal homepage: www.elsevier.com/locate/compositesb

Meso-scale modeling of hypervelocity impact damage in composite laminates Aleksandr Cherniaev*, Igor Telichev Department of Mechanical Engineering, University of Manitoba, E2-327 EITC, 75A Chancellors Circle, Winnipeg, MB R3T 5V6, Canada

a r t i c l e i n f o

a b s t r a c t

Article history: Received 30 April 2014 Received in revised form 25 November 2014 Accepted 12 January 2015 Available online 20 January 2015

Objectives: This paper presents an approach to numerical modeling of hypervelocity impact on carbon fiber reinforced plastics (CFRP). Methods: The approach is based on the detailed meso-scale representation of a composite laminate. Material models suitable for explicit modeling of laminate structure, including fiber-reinforced layers and resin-rich regions, are described. Two numerical impact tests with significantly different impact energies were conducted on thermoplastic AS4/PEEK materials with quasi-isotropic layups. Simulations employed both SPH and Finite element methods. Results: Results of simulations are verified against experimental data available from the literature and demonstrate good correlation with the experiments. Conclusions: Developed modeling approach can be used in simulations where post-impact damage progression in composite material is of the main focus. © 2015 Elsevier Ltd. All rights reserved.

Keywords: A. Polymer-matrix composites (PMCs) B. Impact behavior B. Delamination C. Numerical analysis

1. Introduction Composite materials are being used extensively for space applications. During operation they can be subjected to hypervelocity impacts (HVI) of orbital debris at an average velocity of 11 km/s. A number of experimental programs have been conducted over the last 25 years to study the behavior of composites under HVI [1e9]. However, due to limitation on the projectile speed which can be achieved using the existing laboratory equipment (typically, no higher than 8 km/s) and extreme expensiveness of the physical experiments, modeling became an important tool for study of hypervelocity impacts. Most papers related to modeling of HVI on fiber-reinforced plastics deal with particular type of application, such as spacecraft shielding (e.g., [10,11]). Composite materials in those works are represented as macroscopically homogeneous orthotropic media (referred as macro-scale or “homogenized laminate” approach in the following). This technique allows simulating adequately main HVI phenomena accompanied the shielding systems perforation and calculating the overall size of the impact hole or crater in composite. However, in order to obtain more detailed information about local damage and its progression in a composite, models of higher

* Corresponding author. Tel.: þ1 204 296 94 83. E-mail address: [email protected] (A. Cherniaev). http://dx.doi.org/10.1016/j.compositesb.2015.01.010 1359-8368/© 2015 Elsevier Ltd. All rights reserved.

resolution are needed. It is especially the case for spacecraft loadbearing structural members (e.g., composite truss tubes) and pressurized components of onboard systems (e.g., composite pressure vessels) where the impact damage formation is governed by superposition of impact and quasi-static loading. Here mechanisms of post-impact fracture progression and possible failure of impactdamaged composite structure is of the main interest rather than just evaluation of the size of a perforation hole. Thus, higher level of detail is required for modeling of impact on such components. Moreover, it was found in Ref. [12] that for filament-wound composites delamination propagation may be restrained by interweaving of filament bands and depends on the filament winding pattern used. Then adequate capturing of material damage and fracture requires explicit representation of complex winding patterns in numerical model. Macro-scale approach cannot be utilized in this case and, thus, is not applicable for simulating HVI-induced damage in general type of continuous fiber/polymer matrix composites. In this paper we present more general approach based on mesoscale representation of composite laminates, including modeling of individual plies and resin-rich regions between them. This technique can be classified under “continuum approach” in damage mechanics (e.g., [13]). Although we do not deal directly with filament-wound composites in the current paper, the technique we describe here can also be applied to this type of materials with no additional modifications required.

96

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

Our modeling approach is presented in Section 2 of the current paper. It is then used to simulate damage produced by aluminum projectiles in AS4/PEEK laminates under conditions similar to those used in experiments conducted at the University of Kent at Canterbury (UKC data) and NASA Johnson Space Center (JSC data) and reported in Refs. [5] and [7]. Two HVI tests used for simulations were chosen to represent impacts with different kinetic energies. Results of simulations compared with the experimental data in terms of predicted and measured delamination area. ANSYS AUTODYN-3D hydrocode was used for all simulations presented in this paper.

Fracture of a composite includes both intralaminar and interlaminar mechanisms. In the former case, it is confined to anisotropic linearly elastic fiber-reinforced layers of a laminate and is brittle in nature for the most of existing CFRP systems. Mechanisms of intra-ply damage include fiber breakage and matrix cracking along fibers. Interlaminar fracture or delamination is attributed to relatively tough isotropic polymer matrix concentrated in the resinrich regions (RRR) between each pair of fiber-reinforced layers. Deformation in these regions prior to failure initiation is determined by the properties of the polymer resin. Subsequent damage and fracture may be either due to deformation of the matrix within RRR (cohesive fracture), or due to debonding on the interface between resin-rich region and fibers adjacent to it (adhesive fracture) [14]. To reproduce the above mentioned fracture mechanisms more precisely, each composite laminate in our simulations has been represented as a structure consisting of alternating fiber-reinforced and finite-thickness resin-rich layers. By doing that, we uncouple the description of behavior of two physically different components of a laminate and employ more suitable material models for each of them. This introduces higher level of realism in modeling compared to the case of macro-scale representation. To be suitable for hypervelocity impact simulations, in general, materials being modeled has to be characterized in terms of 1) stress-strain relations; 2) equation of state (EOS); 3) failure initiation criteria; 4) post-failure response model. Two next sub-sections describe material models employed in our study to simulate behavior of fiber-reinforced layers (Section 2.1) and interlaminar regions (Section 2.2). 2.1. Modeling of fiber-reinforced layers

C12 C22 C32 0 0 0

C13 C23 C33 0 0 0

0 0 0 C44 0 0

0 0 0 0 C55 0

(2)

Defining the average strain increment as a third of the trace of the strain tensor

1 ðDε11 þ Dε22 þ Dε33 Þ; 3

(3)

and approximating, for small strain increments, the volumetric strain increment as

Dεvol zDε11 þ Dε22 þ Dε33 ;

(4)

the total strain increments can be expressed in terms of the volumetric and deviatoric strain increments, which results in the following orthotropic constitutive relation:

2

3 2 C11 Ds11 6 Ds22 7 6 C21 6 7 6 6 Ds33 7 6 C31 6 7 6 6 Ds23 7 ¼ 6 0 6 7 6 4 Ds31 5 4 0 Ds12 0 2

C12 C22 C32 0 0 0

C13 C23 C33 0 0 0

0 0 0 C44 0 0

0 0 0 0 C55 0

3 1 d Dε Dε þ 6 11 3 vol 7 6 7 7 3 6 6 d 7 1 0 6 Dε22 þ Dεvol 7 6 7 7 3 0 7 6 7 6 7 0 7 1 7; 7$6 d Dε Dε þ 6 7 0 7 6 33 3 vol 7 7 7 0 5 6 6 7 Dε23 6 7 C66 6 7 6 7 Dε 31 4 5 Dε12 (5)

Defining pressure increment as one third of the trace of the stress increment tensor

1 Dp ¼  ðDs11 þ Ds22 þ Ds33 Þ; 3

(6)

and substituting Dsii in (6) from (5), yields:

Each layer in our study was considered as a homogeneous orthotropic material. In this case, incremental stress-incremental strain relations can be represented as follows:

3 2 C11 Ds11 6 Ds22 7 6 C21 6 7 6 6 Ds33 7 6 C31 6 7 6 6 Ds23 7 ¼ 6 0 6 7 6 4 Ds31 5 4 0 Ds12 0

Dεij ¼ Dεdij þ Dεave

Dεave ¼

2. Modeling approach

2

this model, strain increments (Dεij ) in equation (1) can be split into their average (Dεave ) and deviatoric (Dεdij ) components in the following manner:

3 2 3 Dε11 0 6 7 0 7 6 Dε22 7 7 6 7 0 7 7$6 Dε33 7 6 7 0 7 6 Dε23 7 7 0 5 4 Dε31 5 Dε12 C66

(1)

Important aspect of the hypervelocity impact physics is the propagation of shock waves through the interacting materials. As, in general, response of composites to shock loading is nonlinear, incremental stress-incremental strain relationships described by equation (1) should be modified. The aim of such modification is to separate volumetric response of the composite from its ability to carry shear stresses and to describe the former in terms of nonlinear equation of state (EOS). In this part, we employed in our study constitutive model available in AUTODYN [15]. According to

1 Dp ¼  ½ðC11 þ C22 þ C33 þ 2ðC12 þ C23 þ C31 ÞÞDεvol  9 1 1  ½C11 þ C21 þ C31 Dεd11  ½C12 þ C22 þ C32 Dεd22  3 3 1  ½C13 þ C23 þ C33 Dεd33 ; 3 (7) Here multiplier of volumetric strain increment can be recognized as “effective” bulk modulus (K 0 ) of orthotropic material (following the analogy with isotropic bulk modulus):

1 K 0 ¼  ½ðC11 þ C22 þ C33 þ 2ðC12 þ C23 þ C31 ÞÞ 9

(8)

Finally, equation (7) can be modified to take into account nonlinear dependence of pressure on volumetric strain:

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

 1 h s2 ¼ $ ðs11  s22 Þ2 þ ðs22  s33 Þ2 þ ðs33  s11 Þ2 þ 6$ s223 2 i þ s231 þ s212 ;

Dp ¼ DpðDεvol Þ 1 1  ½C11 þ C21 þ C31 Dεd11  ½C12 þ C22 þ C32 Dεd22  3 3 1  ½C13 þ C23 þ C33 Dεd33 ; 3

(15) (9)

where the pressure contribution from volumetric strain pðεvol Þ can be expressed in the polynomial form:

p ¼ K 0 εvol þ A2 ðεvol Þ2 þ A3 ðεvol Þ3 þ ðB0 þ B1 εvol Þr0 e;

(10)

whereA2 , A3 , B0 , B1 are material constants, r0 is the initial density, e is internal energy per unit mass. As it has been noted above, intralaminar failure is brittle for most of CFRP systems, thus it was considered as instantaneous for the corresponding direction once one of the following criteria has been satisfied:

 e2 ¼

s11 XT

2 þ

s 2 12

S

þ

s 2 13

S

1

(11)

for the fiber failure, and

 e2 ¼

s22 YT

2 þ

s 2 12

S

1

(12)

for the transverse matrix cracking. Here XT and YT are the longitudinal and in-plane transverse strength of a composite, correspondingly; S is the in-plane shear strength. Failure due to throughthe-thickness loading was attributed to the interlaminar region and is described in the following subsection.

2.2. Modeling of interlaminar regions Behavior of the isotropic resin-rich regions prior to failure initiation was modeled as elasto-plastic, with the von Mises J2 function determining the yield surface. Due to software limitations, however, it is not possible in AUTODYN to combine directly von Mises strength model and Crack Softening algorithm (used to simulate interlaminar failure; described below) within a single material model. This problem was solved by employing generalcase anisotropic yield model available in AUTODYN [15] and reducing it to the von Mises yield function by the careful selection of plasticity parameters (aij ), as it was recommended in Ref. [16]. In this model, the yield surface is given as

  f sij ¼ a11 s211 þ a22 s222 þ a33 s233 þ 2a12 s11 s22 þ2a23 s22 s33 þ 2a13 s11 s33 þ þ2a44 s223 þ 2a55 s231 þ 2a66 s212 ¼ k

97

which represents the general form of the von Mises yield criterion. Following Refs. [15] and [16], relationship between plastic strain p increments dεij and stress increments in this plasticity model is defined in terms of associated flow rule, and can be written out explicitly as

2

3 p 2 3 dε11 2a11 s11 þ 2a12 s22 þ 2a13 s33 6 p 7 6 dε22 7 6 2a12 s11 þ 2a22 s22 þ 2a23 s33 7 6 p 7 6 7 6 dε 7 6 2a13 s11 þ 2a23 s22 þ 2a33 s33 7 6 6 33 7 7$dl ¼ 6 dεp 7 6 7 4a44 s23 6 23 7 6 7 6 p 7 4 5 4a s 55 31 4 dε31 5 p 4a66 s12 dε 12

where dl is the plastic factor, which can be determined from a simple tensile test stress-strain curve. Equation (16) can be reduced for the case of isotropy by substitution of aij given by (14). Effective plastic strain increment dεp (scalar parameter) can be introduced using concept of plastic work: p

dW p ¼ sij $dεij ¼ s$dεp ;

(17)

Using (13), (16) and (17) it can be shown that at any timestep effective plastic strain can be found from 2

ðdεp Þ ¼

8   f sij $ðdlÞ2 3

(18)

If stresses at any timestep are found to lie outside the yield surface defined by equation (13), the backward-Euler algorithm gets activated to return the stress state back to the yield surface (see Ref. [15] for details of its implementation). Once the yield function is satisfied, the plastic portion of the strain increment can be determined from equation (18). To predict failure initiation in the interlaminar region, the following interactive quadratic criterion has been used:

 e2 ¼

s33 ZT

2

 þ

s23 SILSS

2

 þ

s13 SILSS

2 1

(19)

As, in general, failure may initiate within matrix, as well as on the fiber-matrix interface, the strength constants in (19) correspond to the through-the-thickness tensile ZT and interlaminar

(13)

Here k is the hardening parameter depending on plastic deformation and is defined as k ¼ 23s2 (see Refs. [15], [16]), where s is the effective stress. The above function reduces to the von Mises yield criterion if the following values of plasticity parameters are used:

2 1 a11 ¼ a22 ¼ a33 ¼ ; a12 ¼ a23 ¼ a13 ¼  ; a44 ¼ a55 ¼ a66 ¼ 1 3 3 (14) Substituting definition of k and plasticity parameters given by (15) into (14), one can obtain:

(16)

Fig. 1. Schematic illustration of crack softening algorithm in AUTODYN.

98

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

shear SILSS properties of the composite, rather than the properties of the bulk polymer resin. Formation of the discontinuities in the interlaminar region following the failure initiation was modeled implicitly, according to the principles of damage mechanics. Technique employed here is analogous to Crack Softening algorithm available in AUTODYN [15] in the case when it is applied to the 33-plane only. It assumes linear softening of the RRR material due to microcracking, so that its behavior schematically can be represented as it is shown in Fig. 1. Here the area A under the softening portion of the stressestrain curve is related to the fracture energyGf , asA ¼ Gf =L, where L is a characteristic cell dimension in the direction of failure and is included to reduce the mesh dependency of the solution. The ultimate strain in this case can be calculated as

εu ¼

2Gf sfail L

Analogous expressions can be written for other damage variables ðD23 ; D31 Þ. Here C is the damage coupling coefficient, which represents the possibility that damage due to tensile loading may reduce material capability to resist shearing, and vice-versa. Its value may vary between 0 and 1. The value of C ¼ 0:2, recommended in AUTODYN [15], was employed in this study. As there is no direct translation of the fracture toughness of a bulk resin into the interlaminar fracture toughness of composite (see e.g. Ref. [14] for discussion on this topic), values of GIc and GIIc of composite should be used with the above damage model, rather than strain energy release rates corresponding to pure matrix. It was assumed that EOS of thin resin-rich regions would not have any noticeable effect on overall composite behavior under hypervelocity impact loading. As a result, simple linear EOS of the following type was employed:

(20)

Dp ¼ K$Dεvol ; where K is a bulk modulus of the polymer resin.

and gradient of linear softening can be found as

sfail L 2Gf

h¼

(21)

3. Validation cases

If at any timestep after failure initiation the state of stress is found to lie outside of failure surface defined by criterion (19), a backward-Euler algorithm is used to return the stress back to the surface. In this part, the procedure is analogous to the handling of material plastic behavior. The resulting inelastic increment of strain is then accumulated as crack strainεcr . Softening in this case may be expressed in terms of damage parameter D, such as

D¼

εcr $h sfail

3.1. Experimental data A number of experimental programs have been conducted over the last 25 years to study the behavior of composites under HVI. Many experimental results from different sources were gathered and statistically analyzed by Tennyson and Lamontagne [7]. Their statistical analysis, carried out specifically for carbon fiber reinforced plastics, resulted in establishment of empirical correlations between projectile kinetic energy, as well as target and projectile material properties from one side, and resulting hypervelocity impact damage from the other. The correlations determine diameter of a perforation hole/crater ðDcr Þ and “equivalent” diameter of damaged zone obtained by ultrasonic c-scan inspection ðDcscan Þ, as follows:

(22)

Just before the failure initiation, when εcr ¼ 0, the damage parameter D ¼ 0; when εcr ¼ εu , the damage parameter takes value of unity ðD ¼ 1Þ. Here crack strain and, correspondingly, damage are tensors with each component of the crack strain tensor associated with corresponding damage variable. Maximum tensile stress at any time after failure initiation can be calculated as

smax ¼ sfail ð1  DÞ

Dcr ¼ 1:06$

(23)

Failure surface in the course of failure progression was updated according to the following criterion [15]:

 e233 ¼

Dcscan

2  2  2 s33 s23 s31 þ þ ZT ð1  D33 Þ SILSS ð1  D23 Þ SILSS ð1  D31 Þ (24)

The damage variable D33 for any timestep n þ 1 can be defined as

εcr εcr $h23 εcr $h31 33 $h33 þ C$ 23 þ C$ 31 s33fail s23fail s31fail

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t$rt 3 ; Ek Dp rp

(27)

sffiffiffiffiffiffiffiffiffiffiffi t 3 ¼ 4:04$ Ek ; Dp

(28)

where Ek e projectile kinetic energy, t e thickness of the target, Dp e diameter of the projectile, rt and rp are target and projectile densities. The Dcscan parameter should be understood as diameter of a circle encapsulating the area equivalent to irregular-shaped delamination area measured by c-scan, including the area of perforation hole. Correlations (27) and (28) have been made through a wide range of CFRP thicknesses and impact energies. Two validation cases from those analyzed by Tennyson and Lamontagne have been selected. They correspond to experiments conducted at the University of Kent at Canterbury (UKC) and Johnson Space Center (JSC). The cases were chosen in such a

1

n Dnþ1 33 ¼ D33 þ

(26)

(25)

Table 1 Experimental data. Shot ID from [Ten]

UKC#4 JSC#2

Target parameters

Projectile parameters

Results

No of plies

Thickn., mm

Layup

Diameter, mm

Velocity, km/s

Energy, J

Dcr, mm

Dc-scan, mm

16 16

1.83 1.96

[0, ±45, 90]s [0, ±45, 90]s

2 2

4.66 6.75

123 258

4.3 5.6

19.5 25.0

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103 Table 2 Mechanical properties of AS4/APC-2 thermoplastic composite (Vf ¼ 60%). Material property

Dimensionality

Value

Density, r Longitudinal Young's modulus, E1 Transverse Young's modulus, E2 In-plane shear modulus, G12 In-plane Poisson's ratio, n12 Longitudinal tensile strength, XT Transverse tensile strength, YT In-plane shear strength, S Flexural Strength

g/cm3 GPa GPa GPa e MPa MPa MPa MPa

1.60 138 10.3 5.7 0.30 2070 86 186 2000

manner to represent different impact energies. Test parameters and experimental results for both cases are shown in Table 1. 3.2. Material properties In both cases targets were made of AS4/PEEK (APC-2) composite, whereas projectiles were of aluminum. The original source presenting results of UKC experiments [5] gives mechanical properties of the composite exactly the same as those provided by the APC-2 resin manufacturer [17]. They are listed in Table 2. Table 3 provides properties of the bulk APC-2 resin. Values of fracture toughness obtained experimentally by Stanley and Pipes [18,19] for AS4/APC-2 composite under dynamic loading were employed in our study. They correspond to GIc ¼ 350 J/m2 and to GIIc ¼ 400 J/m2, and reflect dramatic reduction of the fracture toughness for both first and second opening modes, as compared with the static values provided by manufacturer (1700 J/m2 and 2000 J/m2, correspondingly [17]). A few assumptions regarding composite properties have been made to complete its characterization in terms of established material models. First, interlaminar shear strength, which value is not provided by manufacturer, was determined using the following correlation, given in Ref. [20] specifically for carbon fiber/PEEK composites: SILSS ¼ Flexural Strength=15. This gives the value of SILSS ¼ 135 MPa. Second, Ryan et al. [21] have shown that different types of CFRP systems reveal high level of similarity in their response to shock loading. Having this in mind, we adopted the values of A2 ¼ A3 ¼ 0 and B0 ¼ B1 ¼ 2:5 obtained experimentally in Ref. [21] for UMS2526/Krempel carbon fiber reinforced plastic. Behavior of projectiles has been modeled using combination of shock equation of state, JohnsoneCook strength and failure models with material properties of Al6061-T6 adopted from Ref. [24] (see Table 4). 4. Numerical models The Smooth Particles Hydrodynamics method (SPH) has been used to model adequately deformation and fragmentation of aluminum projectile in the course of impact event [22,23]. The 2 mm-diameter projectile was represented by 8063 SPH nodes (corresponds to “Particle size” of 0.08 mm). Table 3 Mechanical properties of APC-2 PEEK thermoplastic resin. Material property

Dimensionality

Value

Density, r Young's modulus, E Tensile strength, at break Tensile strength, at yield Failure strain, % Yield strain, %

g/cm3 GPa MPa MPa e e

1.32 3.6 100 95 70 5

99

Targets represented by 100  60 mm laminated plates have been modeled using Finite Element method in Lagrangian formulation. Each layer, as well as each resin-rich region, had been modeled as individual solid bodies with adjacent faces tied among themselves using simple “Bonded Contact with NO separation” feature in AUTODYN (see Fig. 2). As no separation was allowed on the interfaces, all of the through-the-thickness damage in the model, characterized by such macromechanical properties as interlaminar strength and fracture toughness, was confined to the RRR elements. In total, laminates consisted of 8 unidirectional plies and 7 resinrich layers. RRL thickness was assumed to be equal to 0.035 mm for all simulations, which is within the typical range of 0.01e0.1 mm. Uniform meshing of the composite plates resulted in 1,440,000 0.25  0.25 mm hexahedral finite elements (96,000 per ply, including resin-rich layers) with one element through the thickness of each layer. Preliminary studies showed that solution became almost mesh-independent if aforementioned element size was employed. To allow modeling of target perforation, numerical mechanism of erosion has been utilized. It automatically removed heavily distorted finite elements from the model when effective strain had reached a user-specified value. The threshold values of effective strain were determined by calibration of simulation results against experimental data from Ref. [7]. In the lack of direct experimental data for particular CFRP material system, the empirical relation (27) can be recommended to conduct such calibration. It was found in the preliminary studies that erosion strain of 0.55 for brittle UD plies used along with ultimate effective strain of 3.0 for ductile RRR material provided acceptable level of correlation between numerical and experimental results. 5. Results and discussion Impact holes predicted for UKC #4 and JSC #2 data sets along with surrounding delamination are exemplified in Fig. 3. Here redcolored elements indicate complete damage in the interlaminar region. Fig. 3 also shows secondary debris clouds resulted from fragmentation of target and projectile. It should be noted that fragments of projectile are modeled with SPH particles, whereas fragments of the composite plate are represented in numerical model by nodes of eroded elements. Therefore, the debris cloud is a

Table 4 Mechanical properties of Al6061-T6 [24]. Material property

Dimensionality

Density, Shock equation of state Guneisen coefficient, G Parameter C0 Parameter S Reference temperature Specific heat, JohnsoneCook strength model Shear modulus, G Yield stress, A Hardening constant, B Hardening exponent, n Strain rate constant, C Thermal softening exponent, m Melting temperature, JohnsoneCook failure model Damage constant, D1 Damage constant, D2 Damage constant, D3 Damage constant, D4 Damage constant, D4

g/cm3

Value 2.703

e m/s e K J/kgK

1.97 5240 1.40 293 885

GPa MPa e e e e K

26,000 324 114 0.42 0.002 1.34 925

e e e e e

0.77 1.45 0.47 0.00 1.60

100

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

Fig. 2. Discretization of the composite plate and Ø2 mm aluminum projectile.

Fig. 3. Predicted impact hole and surrounding interlaminar damage for a) UKC #4 and b) JSC #2 input data.

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

101

Fig. 4. Intralaminar and interlaminar damage predicted for UKC #4 data set.

mixture of both. This technique allows to retain inertia in the simulation upon erosion of finite elements and to use their nodes for modeling of subsequent impact events. Fig. 4 exemplifies predicted intra- and interlaminar damage for the UKC #4 simulation (only damage tensor components of main influence are presented, including D22 and D12 for intralaminar failure, and D33 for delamination). The contour plots for D22, D12 and D33 show the elements which are not capable to resist transverse tension, in-plane shearing and through-the-thickness tensile loading, correspondingly. JSC #2 calculations are not shown for the sake of brevity, as results are qualitatively the same due to similarity of layup. It should be noted that intralaminar failure has been modeled as a brittle phenomenon. As a result, variables D22 and D12 for an element were set to unity immediately after corresponding stress limit had been reached in it. Many experimentalists reported that in unidirectional composites exposed to hypervelocity impacts, perforation was accompanied by the front surface delamination and peeling (see e.g. original paper describing UKC experiments [5], as well as [1], [3] and [9]). This surface damage always had a form of elongated strips emanating from the impact crater in the direction of reinforcing fibers in external layer. Mechanism of such type of damage was also discussed in Ref. [12]. It was explained by shearing along the interface of adjacent plies with dissimilar orientation, produced by longitudinal waves traveling away from the area of impact. Damage pattern presented for the upper layer in Fig. 4 coincides with the description given above. It is characterized by a zone of intralaminar damage and delamination elongated in the direction of fibers in the outer ply. This evidences that the model is capable of

representing experimentally observed phenomenon in a qualitative manner. The quantitative assessment of the model's accuracy has been done by comparison of predicted and experimentally determined (by means of c-scan) delamination area. In order to make numerical results comparable to c-scan measurements obtained in UKC and JSC experiments, damage contours from all resin-rich layers were projected on a single plane. It was accomplished by setting up 50% opacity for the resin-rich layers while contours of D33 damage variable were being displayed. Resulting images are shown in Fig. 5. Damaged areas Adelam were carefully measured and used to calculate the “equivalent diameter of delaminated zone” as qffiffiffiffiffiffiffiffiffiffiffi Ddelam ¼ 4Adelam p . Only completely damaged elements with D33 ¼ 1 were taken into account when calculating the delamination area, as there is no absolute certainty that partly damaged areas could be detected by c-scan equipment. Results of calculations and their comparison with experimental c-scan measurements are presented in Table 5. Good correlation (within ~12%) between numerical and experimental results was found. Sensitivity of the solution to fracture toughness was explored by carrying out an additional pair of computations with the static values of GIc and to GIIc. Although overall damage pattern remained unchanged, the delamination area in this case was significantly underpredicted as compared with experimental data, resulting in equivalent diameter of delamination of 8.9 and 11.8 mm for UKC #4 and JSC #2 experiments, correspondingly. Thus, availability of dynamic fracture toughness data for the composite system of interest is important for accurate predictions of delamination.

102

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103

Fig. 5. Delamination patterns obtained from simulations with dynamic fracture toughness: a) UKC #4 and b) JSC #2 (All damaged surfaces projected on a single plane).

6. Conclusions Detailed progressive damage simulations of hypervelocity impacts on continuous fiber reinforced composites are presented. The

simulations employed combination of SPH and Finite element method in Lagrangian formulation to represent an aluminum projectile and a composite, correspondingly. The principal feature of numerical models is a meso-scale representation of composite

A. Cherniaev, I. Telichev / Composites Part B 74 (2015) 95e103 Table 5 Comparison of numerical and experimental results. UKC#4

JSC#2

Ddelam, mm Simulation Experiment

17.2 19.5

D, %

Ddelam, mm

11.8

Simulation Experiment

D, % 24.1 25.5

5.5

laminates, including explicit modeling of fiber-reinforced plies and resin-rich regions between them. Results of simulations are verified against experimental data available from the literature. Good correlation with experiments has been revealed in terms of predicted delamination area and qualitative representation of external damage. Availability of dynamic fracture toughness data has been found to be essential for accuracy of numerical predictions. The developed meso-scale modeling approach is recommended for the use in simulations where post-impact damage progression in composite material is of the main focus. Acknowledgments This work was supported by a Discovery Grant No. 402115-2012 from the Natural Sciences and Engineering Research Council of Canada. References [1] Yew CH, Rodney BK. A study of damage in composite panels produced by hypervelocity impact. Int J Impact Eng 1987;5:729e38. [2] Schonberg WP. Hypervelocity impact response of spaced composite material structures. Int J Impact Eng 1990;10:509e23. [3] Christiansen EL. Investigation of hypervelocity impact damage to space station truss tubes. Int J Impact Eng 1990;10(1):125e33. [4] Silvestrov VV, Plastinin AV, Gorshkov NN. Hypervelocity impact on laminate composite panels. Int J Impact Eng 1995;17(1):751e62. [5] Lamontagne CG, Manuelpillai GN, Taylor EA, Tennyson RC. Normal and oblique hypervelocity impacts on carbon fibre/PEEK composites. Int J Impact Eng 1999;23(1):519e32.

103

[6] Schonberg WP. Protecting spacecraft against orbital debris impact damage using composite materials. Compos A 2000;31(1):869e78. [7] Tennyson RC, Lamontagne C. Hypervelocity impact damage to composites. Compos A 2000;31(1):785e94. [8] Lamontagne CG, Manuelpillai GN, Kerr JH, Taylor EA, Tennyson RC, Burchell MJ. Projectile density, impact angle and energy effects on hypervelocity impact damage to carbon fibre/PEEK composites. Int J Impact Eng 2001;26(1):381e98. [9] Baluch AH, Park Y, Kim CG. Hypervelocity impact on carbon/epoxy composites in low earth orbit environment. Compos Struct 2013;96(1):554e60. [10] White DM, Tylor EA, Clegg RA. Numerical simulation and experimental characterization of direct hypervelocity impact on spacecraft hybrid carbon fibre/Kevlar composite structure. Int J Impact Eng 2003;29(1):779e90. [11] Wicklein M, Ryan S, White DM, Clegg RA. Hypervelocity impact on CFRP: testing, material modelling, and numerical simulation. Int J Impact Eng 2008;35(1):1861e9. [12] Cherniaev A, Telichev I. Numerical simulation of impact damage induced by orbital debris on shielded wall of composite overwrapped pressure vessel. Appl Compos Mat 2014;21(6):861e84. [13] da Silva LFM, Campilho RDSG. Advances in numerical modeling of adhesive joints. Springer; 2012. [14] Bradley WL. Relationship of matrix toughness to interlaminar fracture toughness. In: Friedrich K, editor. Application of fracture mechanics to composite materials. Amsterdam: Elsevier; 1989. p.159e87. [15] Composite modelling. Revision 1.3. Release 14.0. Autodyn®; 2011. [16] Chen JK, Allahdadi FA, Sun CT. A quadratic yield function for fiber-reinforced composites. J Compos Mater 1997;31:788e811. [17] http://www.cemselectorguide.com/pdf/APC-2_PEEK_031912.pdf. [18] Smiley AJ, Pipes RB. Rate effects on mode I interlaminar fracture toughness in composite materials. J Comp Mats 1987;21:670e87. [19] Smiley AJ, Pipes RB. Rate sensitivity of mode II interlaminar fracture toughness in graphite/epoxy and graphite/PEEK composite materials. Comp Sci Tech 1987;29:1e15. [20] Hartwig G. Fiber composites in cryogenic applications. In: Proceedings of ECCM-8 Conference. Naples; June, 1998. p. 197e206. [21] Ryan S, Wicklein M, Mouritz A, Riedel W, Shafer F, Thoma K. Theoretical prediction of dynamic composite material properties for hypervelocity impact simulations. Int J Impact Eng 2009;36:899e912. [22] Monaghan JJ. Smooth particle hydrodynamics. Ann Rev Astron. Astrophys 1992;30:543e74. [23] Autodyn: explicit software for nonlinear dynamics. SPH user manual & tutorial. Century Dynamics Inc; 2005. [24] Corbett BM. Numerical simulations of target hole diameters for hypervelocity impacts into elevated and room temperature bumpers. Int J Impact Eng 2006;33:431e40.