Prediction of delamination crack growth in carbon/ﬁber epoxy composite laminates using non-local interface damage model

– The use of composite laminates is increasing in these days due to desired directional properties and low densities in comparison of metals. Delamination is a major source of failure in composite laminates where a crack like entity can initiate and propagate between diﬀerent layers of composite laminates under given loading conditions. Damage mechanics based theories are employed to simulate the delamination phenomena between composite laminates. These damage models are inherently local and can cause the concentration of stresses around the crack tip. In the present study integral type non-local damage formulation is proposed to avoid the localization problem associated to damage formulation. A comprehensive study is carried out for the selection of diﬀerent non-local variables. Finite element simulations based on proposed non-local damage models and classical local damage model are performed and results are compared with available experimental data for UD IMS/924 Carbon/ﬁber epoxy composite laminate.


Introduction
Composite laminates are used as a major load bearing structural parts in different industries like avionic and automobile. The driving force behind the increasing use of composite materials is high strength to weight ratio and directional properties. The composite laminates consist of different layers of fibers/fabric stacked together and cured with thermosetting resin [1]. Due to their inherent laminate nature, any two adjacent layers of composite laminates can separate due the presence of delamination crack and may fail under given loading conditions [2,3].
In order to simulate the delamination crack growth in composite laminates, interface modeling techniques are developed using damage and fracture mechanics failure theories under static and fatigue loadings [4][5][6][7]. A lot of work has been carried out at the meso-scale (layers and interfaces) to understand the delamination failure process. Meso-scale, which lies between micro and macro scale consists of two basic constituents, layers and interfaces. The interlaminar interface, which is a mechanical surface, connects two adjacent layers and depends on the relative orientation of their fibres [8,9].
Localization is a common phenomenon associated to already available classic damage models [10,11]. a Corresponding author: hassan605@yahoo.com Localization can cause the concentration of simulated results to certain regions rather a realistic and uniform distribution. One way to avoid localization is to introduce rate dependant damage modeling [12]. In recent times much attention is made to integral type non-local damage modeling in order to reduce the localization effects. This type of modeling is effectively discussed and employed to study the delamination crack growth in composite laminates [13,14].
In the present study a classical damage model proposed by Allix et al. [8] is further modified to as no-local damage model. In integral type non-local damage modeling, one variable is chosen for averaging of the results using Gaussian distribution function. For interlaminar interface damage modeling, three different types of variables like interfacial displacement "U ", damage variable "d" and damage energy release rate "Y d " are available for non-local damage formulation [14]. Finite element simulations on delamination crack growth with different types of non-local variables are performed and compared with each other and with classical local damage model in the present study. All the mathematical details regarding formulation of non-local modeling are comprehensively discussed and presented in this study. This paper is organized as follows: classical local damage model is presented in Section 2. Proposed nonlocal interface damage modeling is detailed in Section 3. Section 4 covers the finite element simulations of mode I delamination crack growth based on classical local and proposed non-local damage models for Carbon/fiber epoxy composite laminate. Finally, some concluding remarks are given in Section 5.

Local damage evolution law
The interface is a two-dimensional surface entity that ensures the transfer of stress and displacement between two adjacent plies [4,8]. This model coupled with the damage mechanics theory can take into account the phenomena of delamination that may occur during mechanical loading of structural parts. The relative displacement of one layer to other layer can be written as [8]: Here, U + and U − are displacement vectors on the top and bottom surfaces. Let the bisectors of the fiber direction be N 1 and N 2 . The direction N 3 is normal to the interface, Figure 1. The deterioration of the interface is taken into account by three internal damage variables d 1 , d 2 and d 3 . Here d 3 is associated with the opening mode (or mode I) and d 1 , d 2 are associated with two shearing modes (mode II and mode III) of the interface.
The relationship between stress and displacement in orthotropic plane of axis can be expressed as: where k 0 1 , k 0 2 and k 0 3 are interface rigidities associated to damage variables in orthogonal directions. Three different damage variables can be distinguished according to three modes of failure. The thermodynamic forces combined with the variables of damage and associated to the three modes of delamination are [8]: where X + represents the positive part of X. The energy dissipated in this model can be written as: Three different damage variables corresponding to three modes of failure are strongly coupled and are governed by equivalent single energy release rate function as follows [15]: where, γ 1 and γ 2 are coupling parameters and α is a material parameter which governs the damage evolution in mixed mode. The local damage evolution law is then defined by the choice of a material function as follows [15]: The damage function is selected of the form [8,15]: where, Y O is threshold damage energy, Y C is critical damage energy, n is characteristic function of material (higher values of n correspond to brittle interface) and Y R is energy corresponding to rupture and can be expressed as: A simple way to identify the propagation parameters (Y C , γ 1 , γ 2 ) is to compare the mechanical dissipation yielded by two approaches of Damage Mechanics and linear elastic fracture mechanics (LEFM). In the case of pure mode situations, when the critical energy release rate reaches its stabilized value at the propagation denoted by G C , comparison of dissipations between Fracture Mechanics and Damage Mechanics approaches leads to [15]: In order to satisfy the energy balance principle of LEFM, the area under the curve of stress-displacement curve for the whole debonding process (DP) obtained through Damage Mechanics formulation is set equal to critical energy release rate G iC , and the following relations for the modes I, II and III critical energy release rates can be written: The area under the curve will always equal to critical energy G C for any mode of delamination (see Fig. 2). For mixed-mode loading situation, a standard LEFM model is recovered as [15]:

Non-local interface damage model
Damage models need to describe strain-softening behavior in order to characterize microscopic degradation processes. A phenomenon called "localization" occurs in structural analyses with strain-softening models. Regularization techniques are therefore essential in Damage Mechanics. The goal of using these techniques is to control the localization. In order to control the localization problem some authors proposed rate-dependent [12] and gradient enhanced [11] damage models. However, in this paper an integral type non-local damage model is proposed. Basically, three different types of non-local models are proposed in this section and their details are given below.

Averaging the damage variable "d "
In this type of non-local damage model, the damage variable d is made non-local by means of a spatial averaging over the surrounding domain V with respect to weight function α 0 [16]: where Here the damage variable d which is made non-local through above equations is calculated initially from equation (2). An isotropic function such as Gaussian distribution function is chosen as weight function over the domain V [13].
This weight function depends on the distance r = x − ζ between points x and ζ and on a parameter l called the internal length scale. This quantity represents a measure of the non-local averaging. The distribution of α 0 with respect to distance r, where r = x − ζ , is shown in Figure 3.

Averaging the variable Y di (i = 1, 2, 3)
Now take the damage energy release rates Y di (i = 1, 2, 3) correspond to three failure modes (see Eq. (3)). Now one can write Now V r (x) can be calculated using equation (12). Once the value of Y di (x) is known for three failure modes then non-local damage variable d can be calculated as a func- Here Y is the equivalent damage energy release rate. Once the value of equivalent Y is known using non-local Y di (x) variables, one can calculate damage variable d using equation (2).

Averaging the relative displacement U i (i = 1, 2, 3)
Another approach can be used by averaging the relative interface displacement expressed through equation (1): Now again V r (x) can be calculated using equation (12). Once the value of U i (x) is known then non-local damage variable d can be calculated as a function of U i (x) i.e. d = d(Y (U i (x))).

Finite element simulations
The proposed non-local damage models presented here are implemented in Cast3M (CEA) finite element code. Compute Compute damage variables d1, d2, d3 Updateandexit Table 2. Y d based non-local damage model implementation.
Updateand exit For the study of delamination, a special interface element of zero thickness has been used [17], which is already implemented in finite element software CAST3M [18]. A four node zero thickness interface element "JOI2" available in CAST3M library [18] is used in the finite element analysis. The "JOI2" element has 2dof/node and two Gauss points of integration per element. Implementation of three different non-local damage evolution laws is made in CAST3M via procedure PERSO1 and using subroutine UMAT [18]. The finite element implementations of all non-local damage models are given in Tables 1 to 3. Figure 7 shows the result of applied load versus delamination crack opening displacement for "d" based nonlocal formulations. The effect of internal length parameter l can be observed in the Figure 8. From figure it can be concluded that non-local damage evolution law behaves like local one for small values of l. This is because for small values of l there will be fewer elements available for making the average of damage variable d. The good value of l can be found from experiments by closely observing the delamination crack front because this material length l defines the size of interaction zone of damage Table 3. U based Non-Local damage model implementation. 1 Compute  localization [20,21]. As the value of l increases the crack initiation point moves farther away and then drops after the crack propagation has started but this crack propagation portion of the curve stabilizes a bit higher in comparison to smaller values of l (see Fig. 7). Figure 8 shows the result of applied load versus delamination crack opening displacement for "Y di " based nonlocal formulations. From Figure 8, it is clear that as the value of l increases, the crack initiation point also moves away but the crack propagation portion of the curve is not as stable as in previous case of "d" based non-local formulations. In fact, results can be considered quite erratic for "Y di " based non-local formulations because of ups and downs of curves for different values of internal length parameter l.
Similarly, Figure 9 shows the result of applied load versus delamination crack opening displacement for "U i " based non-local formulations. The crack initiation Experimental Results     phenomenon occurs in same way as for previous two cases i.e. shown increasing trend with the increase of internal length parameter l. But crack propagation portion of curve seems to be more stable than the previous two cases. In this case, crack propagation regions tend to conform each other for different values of length parameter l. Figure 10 shows the in-plane stresses obtained from finite element analysis of DCB for "d" based non-local damage model for different values of internal length parameter l during crack propagation process. From Figure 10, one can observe the higher amount of stresses with increasing values of l that is also in accordance with results shown in Figure 7.

Numerical Results
Similarly Figure 11 presents the in-plane stresses during crack propagation of DCB specimen for "Y d " based non-local damage model for different values of internal length parameter l. Contrary to "d" based finite element results, one can observe decreasing trend of stresses with increasing values of l. Figure 8 also shows similar trend for "Y d " based non-local formulation, but Figure 8 also shows that reaction force varies erratically as value of l increases.

Conclusion
In the present study, a classical local damage model presented by Allix et al. [8] is further modified to incorporate the non-local effects. Three different non-local interface damage formulations based on averaging of variables d, Y d and U are presented. Finite element simulations of IMS/924 composite laminate for mode I delamination crack growth are performed and results are compared with each other. From the finite element simulations it is clear that d and U based non-local damage models show much better results in comparison of Y d based non-local formulation. For Y d based non-local formulation, reaction force starts varying erratically with increasing values of l. From finite element simulations, it is also clear that non-local modeling is strongly influenced by internal length parameter l. For smaller values of l, the results of local and non-local interface modeling are close to each other because few number of elements are available for averaging in non-local modeling with smaller value of l.