DETERMINATION OF THE INITIATION OF DELAMINATION IN FIBER COMPOSITES

Predicting the initiation of delamination is essential for the design of composite structures, because delamination is a major failure mode of layered composites. The according delamination onset criteria can be evaluated on the basis of stress-strength relations, which requires an accurate representation of the through-the-thickness stress distribution, which is delicate for thin shell-like structures. Thus, in this paper, a solid-shell finite element is utilized, which allows for incorporating a fully three-dimensional, anisotropic, micromechanically motivated material model, still being suited for application to thin structures. Moreover, locking phenomena are cured by using both the EAS and the ANS concept, and numerical efficiency is ensured through reduced integration.


INTRODUCTION
Fiber-reinforced composites are gaining importance in technical applications.Their most beneficial characteristics -the very high Young's modulus and low density-are particularly leveraged in shell-like structures of lightweight constructions.The composites examined in this paper consist of multiple layers, each of which is composed of a woven fabric, with two families of fibers, embedded in a matrix material.Besides this anisotropic structure, the stress-strain behavior of fiber composite materials is highly non-linear.
The majority of models accounting for anisotropic material behavior at finite strains were developed in the field of biomechanics.For instance, axisymmetric orthotropic blood vessels were investigated in [12], whereas biological soft tissues were modeled in [38] on the basis of an incompressible transversely isotropic law for moderate deformations.A description of the transversely isotropic behavior of rubber was presented in [29], and orthotropic constitutive equations were provided in [4] for the simulation of human leg impact problems.More recently, anisotropic material models were proposed incorporating the micromechanical structure of wood, e.g.[21,22].In the present paper, however, the model proposed by Reese in [26] for fiber-reinforced rubber-like composites was adopted, in which the transition from the micro-scale to the macro-scale is formulated in a general manner.Therefore, this model is not restricted to rubber-like materials but also suitable for the carbon fibre-reinforced plastics (CFRP) considered here.
Structural collapse in fiber composite structures is caused by the evolution of either matrix transverse cracking, fiber fracture, or delamination.From these damage modes, the delamination is particularly important, because it drastically reduces the bending stiffness of the structure and promotes local buckling in case of compressive loads.Including delamination into the computation of composite structures requires the definition of an appropriate criterion for its onset as well as the prediction of its growth after an initial crack has evolved.
For the initiation of delamination, different criteria exist, formulated in dependence of stress-resistance relations, e.g.[5,9,11,36,41].After onset of delamination, the high stress gradients appearing at the crack front prohibit employing solely stress-based criteria.Thus, fracture mechanics approaches are often used for simulating the delamination propagation, such as the virtual crack closure technique, [18,19,24,35,42].As an alternative, delamination growth can be treated within the framework of damage mechanics using cohesive zone models, which are incorporated into the finite element simulation by interface elements, e.g.[2,8,10,39].However, in this paper, the onset of delamination is addressed based on stress-resistance relations.
Since fibre-reinforced composites are mostly applied in thin shell-like structures, the element formulation demands providing a suitable shape for thin structures while displaying realistically the three-dimensional stress states.Although shell formulations exist, which take into account the through-the-thickness stretching, see e.g.[6,13,16], the implementation of three-dimensional material models is much simpler in the context of solid elements.On the other hand, the latter typically provide a poor performance when being applied to thin shelllike structures.In particular, there are different locking phenomena to be coped with, which cause an overestimation of the stress state and an underestimation of the deformation.Using solid-shell elements represents one strategy to overcome this problem by combining the advantages of both solid elements and shell elements at the same time.Further, applying the enhanced assumed strain (EAS) concept eliminates the volumetric locking in case of (nearly) incompressible materials as well as the Poisson thickness locking, which occurs in bending problems of shell-like structures due to the non-constant distribution of transverse normal strain over the thickness.
In literature, one can find several solid-shell formulations incorporating the EAS concept, see e.g.[1,27,37], to name only a few.To cure the transverse shear locking, which is present in standard eight-node hexahedral elements, the assumed natural strain (ANS) method is applied.In the context of full integration formulations, the ANS can be found e.g. in [14,15,34], and for reduced integration solid-shell formulations e.g. in [7,[30][31][32].The formulation presented in this paper is based on the works of Schwarze and Reese [30][31][32].
For laminated layered composites, the accurate determination of the through-the-thickness stress distribution was recently investigated by several authors.For instance, in [28] an improved shell formulation was used for this, whereas in [40] and [23] the investigations were based on the solid-shell concept.For a more elaborate literature overview, the reader is referred to the review papers [17,20,25] and the references therein.However, to our knowledge no solid-shell formulations exist, which consider the orthotropic behavior of fiber composites with woven fabric accounting for different fiber directions.

Delamination onset criterion
The onset of delamination can be determined on the basis of stress-strength relations.In particular, delamination occurs in pure interlaminar tension (mode I), pure interlaminar sliding shear (mode II), and pure interlaminar scissoring shear (mode III), if the corresponding interlaminar stress component exceeds the respective maximum interfacial strength.Here, the interlaminar stress components are denoted by σ 33 , σ 13 , and σ 23 , respectively, where the 3direction is normal to the considered interface.Then, the respective interfacial strengths are Z 33 , Z 13 , and Z 23 .
To account for mixed-mode loading, the formulation of the onset criterion should incorporate the interaction of these modes.Due to the lack of available experimental data, failure criteria predicting the initiation of delamination have not been fully validated, and hence only few formulations exist.In this paper, the approach of Ye [41] is adopted, in which a quadratic interaction of modes is assumed: As the formulation presented in this paper is capable of taking into account finite strains, it is important to accurately represent the according stress components.First of all, the stresses calculated by the present solid-shell formulation are expressed by the second Piola-Kirchhoff stress tensor S, which has to be pushed to the current configuration where σ denotes the Cauchy stress tensor.From this, the interlaminar traction σ n and the interlaminar resultant shear τ n can be achieved by denoting the normal vector of the considered interface by n.For consistency, the maximum interfacial strength in tension and resultant shear are referred to as Z σ and Z τ , respectively.Consequently, the condition for delamination onset reads if σ n > 0 : This condition has to be checked in each loading step and in each interface of the laminated composite.Thus, the accurate prediction of the stress components in the interfaces is essential for a reliable prediction of the initiation of debonding.Since the layers are usually rather thin, solid elements are not suitable to achieve sufficient accuracy.To overcome this problem, spline approximations for the through-the-thickness stresses can be applied, as proposed in [5].Even so, using solid elements in the thin shell-like applications should be avoided, and solid-shell elements are preferable.Further, it should be noticed that this kind of criterion is suitable to predict the delamination onset, but delamination growth is not covered.

ORTHOTROPIC MATERIAL MODEL
The fiber composites examined in this paper consist of stacked layers, each of which is composed of a woven fabric embedded in a matrix material.The anisotropic material behavior of such composites is taken into account by using the micro-mechanically motivated model proposed in [26].The basics of the continuum model are summarized in the following.Therein, parameters are chosen to represent approximately the behavior of carbon fibers in an epoxy resin matrix.

Concept of structural tensors
Introducing the deformation gradient F, the deformation of a continuous body is represented by the right Cauchy-Green tensor The characterization of a hyperelastic body is then given by the existence of a scalar potential, which is the stored energy function W = W (C), such that is the second Piola-Kirchhoff stress tensor.In the case of orthotropic material behavior, the energy function W (C) reduces to an isotropic function of C and the structural tensors M 1 and M 2 , which are defined by where the vectors n 1 and n 2 are oriented in parallel direction to the fibers.For a discussion of the theoretical background, the reader is referred to [33] and the references therein.Then, the strain energy function W can be represented in dependence of the following invariants:

Strain energy function
In this work, the anisotropic model from Reese [26] is adopted, which assumes that the fibers do not carry any load in case of compression but only in tension, which is not realistic for the CFRP considered here.Therefore, this model is slightly modified, such that the matrix acts as an elastic continuous support for the embedded fibres.Moreover, the fibre volume fractions 0 ≤ ϕ 1 and 0 ≤ ϕ 2 for the two families of fibres are introduced, where ϕ 1 + ϕ 2 ≤ 1 holds.Except of this, we adopt the mentioned model and use the following strain energy function: Here, W NH denotes the Neo-Hookean part displaying the isotropic case in the small strain regime.The strain energy function is given by The anisotropic behavior of the fabric is introduced by the part Note that in [26] further coupling terms have been introduced, which have hardly influenced the results, and therefore are dropped here.The exponents α i , β i , γ i , (i = 1, 2), and ξ are chosen to be integers larger than 2.

SOLID-SHELL ELEMENT FORMULATION
Using standard solid elements in thin structures would require a very high mesh density to predict the stress distribution with sufficient accuracy.Hence, the solid-shell element formulation is used alternatively to avoid inefficient computations.Furthermore, the different locking phenomena are cured by application of the EAS and ANS concepts.The following is based on the works of Schwarze and Reese [30][31][32].

Finite element framework and interpolation
The solid-shell concept is predicated on the two-field variational functional where g ext denotes the virtual work of the external loading.For this formulation, the total Green-Lagrange strain tensor E is split additively into the compatible part E c and the enhanced part E e coming from the EAS concept: The second Piola-Kirchhoff stress tensor S(E) is a function of the total Green-Lagrange strain tensor E. Note that S depends additionally on a set of internal variables in the elastoplastic case.
In this work, isoparametric eight-node hexahedral finite elements are considered, such that both the position vector of the reference configuration X(ξ) = [X 1 , X 2 , X 3 ] T and the displacement vector U(ξ) = [U 1 , U 2 , U 3 ] T are approximated within the element by using tri-linear shape functions The position vector of the current (deformed) configuration reads Then, introducing D = ∂U/∂ξ, the Jacobian matrices J and J of the reference and the current configuration, respectively, can be written as follows: The columns of J and J represent the covariant base vectors with respect to the reference and current configuration, respectively.The polynomial form of the Jacobian matrices defined in (19) is obtained by means of the polynomials The contravariant base vectors with respect to the initial configuration and the current configuration are denoted by ∂X j e j = j ij e j and Hi = ∂ξ i ∂x j e j = jij e j These represent the rows of the inverse Jacobian matrices J −1 and J−1 , respectively, the coefficients of which are denoted by j ij = (J −1 ) ij and jij = ( J−1 ) ij .With this definition, the Green-Lagrange strain tensor can be written in terms of its Cartesian and covariant components E ij and Ēij = E ξ i ξ j , respectively, Denoting the Voigt notation by ( •) and exploiting symmetry as well as Γ ij = 2 E ij , the latter can be stored into the 6 × 1 vectors These can be transformed one to the other with the relation Ê = T Ê, where

Green-Lagrange strain field
Consequently, the compatible Green-Lagrange strain is written in the form Êc = T Êc .To get the components of the compatible Green-Lagrange strain tensor Êc in polynomial form, both the Jacobian matrix J and the displacement gradient D are split column-wise: Hence, the covariant compatible Green-Lagrange strain components are given by To cure curvature thickness locking, the ANS concept is implemented.For this, the covariant compatible strain terms E c ζζ | ξ K := E K c ζζ are evaluated at the sampling points K = A, ..., D, see Fig. 1, as proposed in [3].
Sampling points of the ANS concept at reference element.
Thereby, the covariant compatible strains can be interpolated within the shell mid plane of the reference element by means of bilinear ansatz functions to be evaluated in the points K = A, ..., D. In consequence, the assumed transverse normal strain distribution reads The present solid-shell formulation utilizes a reduced integration scheme within the shell plane (using one integration point), whereas a full integration is used in thickness direction, which allows for choosing arbitrary numbers of integration points (at least two), see Fig. 2. Thus, all integration points are located on the normal through the center of the element defined by ξ ⋆ := (0, 0, ζ) T .To cure volumetric locking as well as Poisson thickness locking, the EAS concept is adopted.These locking effects are treated on the level of the integration points, which can be expressed by Êe = Ê⋆ e , indicating values to be evaluated in the integration points by ⋆ .Since in (30) the assumed transverse normal strain E ANS e ζζ has been defined independently of ζ, the according value E e ζζ is constructed as being linear in ζ in order to overcome the locking effects, which reads Êe = Ê⋆ e = T 0 B⋆ e W e (32) in which the interpolation matrix B⋆ e = [0, 0, ζ, 0, 0, 0] T requires only one EAS degree-offreedom W e .Here, T 0 is the transformation matrix (26) evaluated in the center of the element.
In order to achieve a polynomial decoupling of the compatible Green-Lagrange strain tensor, Êc = T Êc , a Taylor expansion of the inverse Jacobian matrix is carried out with respect to the center of the element up to the linear terms.The goal is reached by exploiting the Taylor expansion of as well as JJ −1 = I, which must hold for arbitrary ξ = (ξ, η, ζ) T = (ξ 1 , ξ 2 , ξ 3 ) T .Introducing the coefficients of ( 33) into (26), the transformation matrix is given by in polynomial form, where terms of higher order have been dropped, as has been shown to be sufficient in [30].Nevertheless, using the approximation (26), the Cartesian compatible strain still represents a polynomial of high order.For this reason, a Taylor expansion of the Cartesian compatible Green-Lagrange strain is also carried out with respect to the center of the element.The hourglass strain term Êhg c represents an excellent basis for the construction of the hourglass stabilization.However, Ê⋆ c is no longer quadratic in ζ.Thus, the according quadratic term is added:

2nd Piola-Kirchhoff stress tensor
In this section, an efficient stress state Ŝ := Ŝ⋆ + Ŝhg is derived for the reduced integration, which is represented by the second Piola-Kirchhoff stress tensor.Here, Ŝ⋆ is computed at the integration points placed at ξ = ξ ⋆ and must be able to take into account the highly nonlinear stress distributions in thickness direction.The hourglass stress part Ŝhg has to guarantee that the element is free of hourglass instabilities.Following [27], a Taylor expansion of the stress field is carried out with respect to ξ = ξ ⋆ .
The tangent Ĉ⋆ is nonlinear in the thickness direction ζ but independent of ξ and η.In order to enable an analytical integration of the hourglass stabilization terms, Ĉ⋆ is replaced by the deviatoric part of the linear-elastic material tangent Ĉhg , which in Voigt notation reads and which only depends on the artificial hourglass shear modulus µ hg .In elastic problems µ hg is equal to the elastic shear modulus µ = E/(2(1 + ν)), whereas in elastoplasticity it reads as follows: Here, the superscript "dev" indicates the deviatoric part of the considered term.Summing up over the number of integration points (i = 1, ..., n ip ) leads to the following effective hourglass shear modulus where the weighting factors ω i are scaled, such that n ip i=1 ω i = 1 holds.

Discretized weak form
Introducing the virtual form of the Cartesian enhanced Green-Lagrange strain δ Êe together with (41), Eq. ( 14) reads Note that the infinitesimal volume element dV e = J dΩ e has been approximated by means of dV e ≈ J 0 dΩ e .Thereby, the integral Ω e δ Ê⋆ T e Ŝ⋆ J 0 dΩ e vanishes, because (32) is only linear in ζ.Further, the virtual forms of ( 37) and (38) together with (41) are incorporated into (13), which -on the element level-can be written as Again, all terms being linear in the natural coordinates drop out.This leads to the desired decoupling of the parts corresponding to the integration point and the hourglass stabilization.Since the hourglass stress has deviatoric character (see Eqs. ( 41) and ( 42)), the hourglass residual vector R hg u simplifies to in which the integration over the element domain is performed analytically.Further, terms being bilinear in ξη drop out, because Ŝhg does not include any summand depending on that.

NUMERICAL EXAMPLE
The proposed method was applied to a panel with a co-cured stiffener, which had already been investigated in [36].The panel's length and width were 203 mm and 25.4 mm, respectively, while the stiffener's length was 50 mm at the panel's skin and 42 mm at the uppermost layer.The flange was composed of 10 plies with a lay-up of (45 From these mechanical properties, the non-zero parameters for the described model were calculated, taking into account the fibre volume fraction For the interfacial strengths, the following values were adopted giving the following parameters for the delamination onset criterion (5): The panel was clamped at the left end, while it was loaded by a force in longitudinal direction at its right end.Since the onset of debonding was expected at the tip of the stiffener flange, the mesh was refined in this region as illustrated in Figures 3 and 4. In order to incorporate delamination, interface elements were located between all layers, which were furnished with material properties of the matrix alone.The latter were assumed to be ten times thinner than the layers.The described solid-shell element was used for both layers and interfaces, leading to a total number of 9825 solid-shell elements with 25152 nodes.
The location of delamination initiation is shown in Fig. 5, which corresponds to the time step, in which the delamination onset condition was met first.As one can see, delamination occurred at first at the tip of the stiffener flange, as was expected.This location of delamination is reasonable and in agreement with the results presented in [36].
Furthermore, experimental data for this problem can be found in [36].Therein, the maximum displacements from extensometer measurements corresponding to delamination onset are reported to be in the range of 0.11 mm to 0.15 mm.Unfortunately, the exact location of the measurement is not given, which does not allow a quantitative comparison.However, in the current calculation, the computed displacements at delamination onset are in a higher range up to 0.18 mm.This difference can be explained by the fact that residual stresses are present in the specimens, which have not been taken into account in the present calculation.

Conclusion
For many technical applications of fiber-reinforced composites, predicting the onset of delamination is essential for appropriately designing the considered structure.For this, a delamination onset criterion based on stress-strength relations has been suggested in this paper, which requires an accurate representation of the through-the-thickness stress distribution.The proposed solid-shell element is particularly suitable to achieve the required accuracy especially in the thin shell-like applications considered here.The formulation allows for including woven fabrics with two different families of fibers, incorporating a fully three-dimensional, anisotropic, micro-mechanically motivated material model.Concluding, the proposed method is capable of predicting the initiation of delamination of fiber-reinforced composites in shelllike structures accounting for the anisotropic material behavior.