Local zigzag effects and brittle delamination fracture of n-layered beams using a structural theory with three displacement variables

A BSTRACT . Equivalent single layer theories for layered beams effectively and accurately predict global displacements and internal force and moment resultants using a limited number of displacement variables. However, they cannot reproduce local effects due to material architecture or weak/imperfect bonding of the layers, such as zigzag displacement fields, displacement jumps at the layer interfaces and complex transverse stress fields, nor can they simulate delamination damage growth. In this work we will present some applications and discuss advantages and limitations of a recently formulated zigzag model. The model, through a modification of the equilibrium equations of an equivalent single layer theory, which maintains the same number of variables, reproduces local effects and delamination fracture under mode II dominant conditions. The approach is based on a local enrichment of the displacement field of first order shear deformation theory, the introduction of cohesive interfaces and homogenization. Advantages and limitations of the approach have been discussed. Local stress and displacement fields are accurately predicted also in very thick and highly anisotropic beams and wide plates with continuous imperfect interfaces and delaminations. Crack propagation under mode II dominant conditions has been successfully reproduced in various fracture specimens and bend beams: unidirectionally reinforced End-Notched Flexure and End-Load Split specimens; End-Load Split specimen made of a 46-layers cross-ply laminate; isotropic cantilever beams with multiple delaminations interacting during propagation. The solution is analytical or semi-analytical, requires only in-plane discretization of the problem and the displacement unknowns are only three as in first order shear deformation theory.


INTRODUCTION
he mechanical response of layered structures, such as laminates, sandwich composites, laminated glass, and laminated wood, is controlled by local effects due to the elastic mismatch of the layers, the presence of thin interlayers (adhesives, resin rich regions, polymeric layers), and interfacial imperfections, defects and delaminations, which are typical consequences of manufacturing processes. Zigzag displacement fields and complex transverse stress distributions, which are due to the multilayer architecture, and displacement discontinuities, due to the presence of interlayers or interfacial damage/delaminations, are not captured by the widely used equivalent single layer theories (ESLT), such as first order shear deformation theory (FSDT). These theories assume a C 1 continuous displacement field through the thickness which introduces discontinuities in transverse stresses and the impossibility to satisfy local equilibrium. Modeling local effects requires the use of layerwise, discrete-layer or 2D/3D finite-element models, where layers and interfacial regions between layers are distinctly modelled and continuity of displacements, but not of their derivatives (slopes), is imposed at the interfaces [1]. Modeling delamination evolution using these techniques requires in addition a T cohesive zone approach to capture the displacement discontinuities and an in-plane discretization to follow the evolution of the cracks. The main drawback of these approaches is that the number of displacement variables is related to the number of mathematical layers and to the kinematics assumed to describe the layers. To model a n-layered beam with perfectly bonded layers and transverse inextensibility using Timoshenko beam theory in each layer leads to a total of (n+2) displacement unknowns (namely (2n+1) initial displacements minus (n-1) displacement continuity conditions at the interfaces). Delamination damage evolution in the beam can be studied by releasing the displacement continuity at the layer interfaces and introducing cohesive interfacial tractions, normal and tangential, and cohesive traction laws which relate the tractions to the relative sliding and opening displacements [2][3][4]. This approach requires 3n unknowns for each beam region (or (2n+1) unknowns for mode II dominant problems, where crack opening is neglected and a single displacement variable describes the transverse displacements of the layers). An effective approach for modeling multilayered structures was proposed in the early works in [5,6] and is known in the community as "zigzag approach". The basic idea of the global-local zigzag approach is to enrich the displacement field of an equivalent single layer theory (ELST) by local zigzag functions [5,6, and subsequent works] in order to account for the effects of the multilayer structure (zigzag). This adds a large number of local kinematic variables to the global variables, which are then removed from the problem by imposing continuities conditions on interfacial tractions at the layer interfaces [6]. Variational techniques are then used to derive the equilibrium equations, which depend on the global variables only. Other refined approaches have been proposed, e.g. [7,8], which introduce additional global variables to try and overcome some drawbacks of the original theories. The presence of thin interlayers and interfacial damage has been modelled following two different techniques. The first, known as "compliant layer concept", is based on a description of thin interlayers or interfacial damaged regions as regular layers of the stack, with elastic constants apt to describe their response; the method then essentially follows the procedure described above for fully bonded structures [9,10]. The second technique, known as "spring-layer approach", introduces zero thickness interfaces and interfacial traction laws which describe the response of the interlayer/damaged regions [11][12][13][14]. In addition to the zigzag functions, the displacement field of the global equivalent single layer theory is then enriched by interfacial jumps. Considering the n-layered plate in Fig. 1a, the displacements ( ) k v  (=1,2,3) in an arbitrary layer k , with upper and lower coordinates  are displacement jumps at the layer interfaces [14]. Higher order equivalent single layer models can be used which increase the number of the global variables and better describe the transverse stresses in fully bonded problems [1]. Zigzag functions and displacement jumps are then defined by imposing continuity of transverse shear and normal tractions at the layer interfaces and the interfacial constitutive laws to relate tractions and jumps. This yields the macro-scale displacement field and homogenized equilibrium equations which depend on the global variables only. Using this second approach care must be put on properly accounting for the strain energy absorbed during deformation within the interface [15] in the variational equations used for the derivation of the equations of motion. This term was neglected in the early models in [11][12][13] and corrected only later in [14]. Zigzag models have been formulated to study interfacial damage evolution using the compliant-layer concept and continuum damage approaches in [9,16,17]. Recently, the model in [14], which uses a cohesive zone approach and the zigzag kinematic approximation in Eqn. (1), has been applied to simulate single and multiple delamination growth under mode II dominant conditions in wide plates within the framework of fracture mechanics [18].
In this paper the model in [14] is briefly reviewed, the equations which describe the response of n-layered beams with imperfect interfaces and delaminations under mode II dominant conditions, as shown in Figs. 1b,c, are recalled and some applications of the model, novel and from the literature, are presented and discussed. The applications are chosen in order to highlight the applicability of the homogenized approach to predict local stress and displacement fields in layered beams with imperfect interfaces and delaminations and to analyze single and multiple delamination crack propagation in homogeneous and layered beams with different boundary conditions. The results are used to discuss advantages and limitations of the homogenized approach.  v is the critical energy release rate or mode II fracture toughness of the interface (after [18]).

DELAMINATIONS
n this section the equations of the model formulated in [14] for plates with weakly bonded layers, oriented along the geometrical axes, and delaminations shown in Fig. 1a are particularized to a beam (wide plate) with imperfect sliding interfaces and delaminations and used to describe the two schematics shown in Figs. 1b and 1c. Assuming transverse inextensibility and the absence of normal separation between the layers, in order to describe mode II dominant problems with layers in constrained contact, the displacement field in Eqn. (1) modifies since the dependence on 1 x is removed and n  layer interfaces (the shear tractions are calculated using the elastic constitutive laws in Eqn. (7) in the Appendix). The remaining local unknowns, which describe the interfacial jumps, depend on the mechanisms acting at the interfaces which are described through piece-wise linear interfacial traction laws relating interfacial shear tractions and jumps. To describe the beam in Fig. 1b, with continuous imperfect interfaces, such as those due to thin elastic interlayers (adhesives, resin rich layers), the linear traction laws shown in Fig. 2a,  (7), the macroscale displacement field then takes the form: (1) (1) 23 22 3 23 where ( ) 23 i G is the shear modulus of the layer i. In a beam with layers having the same elastic properties, as a unidirectionally reinforced laminate, ( 1) 22 3 ( ) 0 k S R x  and the displacement field coincides with that of first order shear deformation theory. If the elastic constants of the layers differ, in the intact portions of the plate where 1/ 0 k S K  the displacement field coincides with that assumed by the original zigzag theory in [6]. The low order of the global model and Eqns. (2), (7) imply that the shear strains obtained from (2) through compatibility are constants in and between the layers in the intact regions and vanish in the delaminated regions. Following the approach which is commonly used for the structural low order theories, accurate shear stresses and strains can be obtained a posteriori through the imposition of local equilibrium: ( ) On the other hand, the absence of shear strains in the delaminated regions induces under-predictions of the transverse displacements in shear deformable plates and the need for corrections or adjustments [19,20]. The problem will be discussed later.
Variationally consistent equilibrium equations and boundary conditions are derived through the principle of Virtual Works [14]. They are shown below for a beam subjected to distributed transverse load, 3 2 ( ) f x : where 22  , , , , C the coefficients of the 6×6 stiffness matrix in the layer k ; the remaining coefficients are given in (10) in the Appendix; 44 k is a shear correction factor which can be set equal to 5/6, as for a homogeneous beam, since the effects of the layered structure are already accounted for through the zigzag functions. The equilibrium equations in terms of global displacement variables (for a beam subjected to transverse surface tractions, 3 2 ( ) f x ) are [19]: They return the equilibrium equations of first order shear deformation theory in a n-layered beam with equal fully bonded layers, where all coefficients with upper index S vanish. And they return the equilibrium equations of the zigzag theory in [6] for fully bonded unequal layers.

LOCAL STRESS AND DISPLACEMENT FIELDS AND DELAMINATION FRACTURE OF N-LAYERED BEAMS
n this section various applications of the model are presented in order to highlight its main features.

Stress and displacement fields in n-layered beams subjected to mechanical loadings
The equilibrium Eqns. (6) with boundary conditions in (4) are easily solved for n-layered beams and wide plates with continuous imperfect interfaces, shown in Fig. 1(b); the model is applicable to beams with various boundary conditions (simply supported, clamped) and loadings (distributed, concentrated loads). Exemplary closed form solutions have been presented in [19] for simply supported plates under thermo-mechanical loadings; in [21] closed form solutions have been derived for the wave propagation problem. The model accurately captures the transverse displacements in fully bonded beams while requires a correction in shear deformable beams with fully debonded or weakly bonded interfaces where the I effects of the shear strains on the global transverse displacement are not fully captured. The local fields are accurately captured also in very thick beams and wide plates. Predictions are less accurate near clamped edges, due to the presence of boundary regions which will be discussed later [19,20,22]. The diagrams in Fig. 3 show transverse shear stresses through the thickness in a unidirectionally reinforced orthotropic wide plate with fibers (L) along 2 x , two layers (n=2) and one interface. Results in the figure refer to a simply supported plate subjected to a sinusoidal transverse load. The solutions of the homogenized zigzag model are compared with exact elasticity solutions for perfectly bonded, partially bonded and fully debonded layers, after [14]. Geometry and elastic constants are defined in the caption. Fig. 4 shows axial displacements and transverse shear stresses in a thick wide-plate with three symmetrically oriented orthotropic layers (0/90/0) and two equally spaced weak interfaces. The boundary conditions are the same of the previous case. The homogenized model accurately captures the local effects due to the material architecture, e.g. the zigzag displacements, the jumps at the layers interfaces and their effects on transverse shear stresses.     (modified after [19]).

Single and multiple delamination fracture in n-layered beams
The equilibrium Eqns. (6) can be applied directly to study the response of beams (wide plates) with stationary delaminations, which are regions where the layers are fully debonded, as shown in Fig. 1c. For these problems a discretization in the axial direction is required, as shown in Fig. 1c2 since the coefficients in Eqns. (6) and (4) differ between regions. Continuity conditions between the different regions are then applied and the problem is solved analytically or semi-analytically. Displacement and stress fields accurately reproduce the 2D elasticity solutions but for boundary regions, which occur near clamped edges or crack tip cross sections, and will be discussed later in the paper. The energy release rate for the collinear propagation of the delaminations can be calculated using different techniques, namely the compliance method and the Jintegral along different paths [18]. Different fracture problems have been examined in [18]: bend beams with two and three layers (sandwiches), with single and multiple delaminations and various boundary conditions (simply supported and clamped). Here two relevant applications after [18] are recalled and novel applications are presented in order to highlight the capability of the model, which uses only three displacement variables as a classical single layer theory, to describe the discrete event of brittle delamination fracture, to reproduce the interaction effects of multiple delaminations and to follow the evolution of multiple cracks, in homogeneous and layered beams. . Exact solution in [27]. Fig. 5a shows the macro-structural response of the unidirectional End-Notched Flexure (ENF) specimen shown in the inset, through the critical load for crack propagation versus load-point deflection (after [18]). The layup and elastic constants are described in the caption. The energy release rate has been obtained using the compliance method and the transverse displacement calculated through Eqns. (6). Local effects in the layers due to the presence of the delamination are well captured by the model and the distribution of the shear stresses is qualitatively similar to that presented in Fig. 3. This is a pure mode II problem and the dimensionless critical load is obtained in closed form and is

Crack propagation in unidirectionally reinforced End-Notched Flexure and End-Load Split specimens
with IIC G the critical fracture toughness [18]. The critical load coincides with predictions using classical structural mechanics assumptions and a discrete-layer approach. The exact solution of the problem [23] has two additional terms, which account for the effects of crack tip shear on the near tip deformations and are important when the crack is short. The response diagram is compared with the experimental results on a graphite/epoxy [0] 24 laminate tested in [24]. The two theoretical curves shown in the diagram have been obtained using the average values of the elastic constants and energy release rate, red dashed lines Model (a), and the maximum values, black solid lines Model (b). The experimental results, obtained under displacement, control show a load drop in the critical load at the onset of propagation; this is due to an unstable propagation of the crack which grows catastrophically and arrests near the mid-span. The homogenized model, which is under cracklength control [18,25], is able to capture the snap-back instability and follow the virtual branch where crack growth is associated to a reduction of the load-point displacement. Crack propagation is modelled also in the region beyond the midspan to show that the curve stably approaches the limiting solution (dotted line) corresponding to two fully delaminated layers. Fig. 5b shows the critical load for crack propagation in an End Loaded Split (ELS) specimen. The material is a unidirectional E-glass-epoxy laminate with elastic constants given in the caption of Fig. 5b. This is a pure mode II problem. Local effects are similar to those of the ENF specimen and described in Fig. 3. The critical load for crack propagation obtained using the homogenized model is / The solution coincides with predictions using discrete models and elementary beam theory. An accurate solution of the ELS specimen based on dimensional analysis, orthotropy rescaling and 2D finite element analyses has been obtained in [26] and has been recently confirmed by the analytical solution in [27]. For the E-glass-epoxy examined here, the exact solution is /

Crack propagation in End-Load Split cross-ply laminate [(02/90)7/02//02/(90/02)7]
The diagram in Fig. 6 shows the critical load for crack propagation versus load point displacement in the End Load Split (ELS) specimen tested in [28]. The material is a cross ply laminate with 46 unidirectional carbon plies of thickness 0.13 mm, layup [(02/90)7/02//02/(90/02)7] and elastic constants given in the caption of Fig. 6. The specimen has total length 160 T L  mm , width 20 B  mm and initial crack length 0 60 a  mm. The theoretical results have been obtained by assuming that a part of the specimen, of length 30 mm, has been clamped so that the free length is 130 L  mm (this is necessary to match 2D finite element results to the slope of the experimental linear branch; using T L would produce unacceptably large displacements). The diagrams show the critical load for crack propagation versus load-point deflection obtained using the homogenized model, the compliance method, the crack propagation criterion II IIC G = G and crack length control [18,25]. The results are compared with the experimental results in [28] and 2D finite element results obtained here. The finite element results have been obtained by using ANSYS, the VCCT technique, the crack propagation criterion II IIC G = G and a crack length control. The finite element model uses plane strain isoparametric and quadrilateral 8-noded elements to mesh each layer of the beam with elements of equal size 0.075 0.13  mm in order to ensure convergence on both load point displacement and energy release rate. The results of the homogenized model agree well with the finite element predictions and the relative error on the critical load is always less than 3%; this error is due to an underestimation of the compliance in the homogenized model which neglects shear deformations along the delaminated portion of the specimen. The theoretical curves are presented for two values of the mode II fracture energy 0.45 IIC  G N/mm and 0.5 IIC  G N/mm. The first value better matches the experimental results at crack initiation and the second the post-critical branches; this difference is probably due to the presence of nonlinear cohesive mechanisms which are not reproduced by the linear elastic fracture mechanics solutions. The theoretical curves, obtained under crack length control, show snap-back instabilities, with virtual branches having negative slope, which are also evident in the almost vertical drops of some of the experimental curves obtained using displacement control. In the simulations in [28] the fracture toughness was assumed as 0.7 IIC  G N/mm; according to the LEFM analyses, this large value of IIC G could be explained only if crack propagation occured under large scale bridging, with a large process zone, during the entire loading process. Further analyses are required to clarify this point; however propagation under extensive large scale bridging appears to be supported by the experimental measurements of the traction free crack lengths during propagation, which are shorter than those obtained using linear elastic fracture mechanics (not shown here).  Fig. 7 (after [18]) is used to highlight the capability of the homogenized model to analyze multiple delamination fracture and reproduce the effects of the interaction between delaminations. The dimensionless diagrams depict the critical load for the propagation of the cracks in the homogeneous and isotropic cantilever beams in the insets in Fig. 8 [18] are compared with the results of the discrete-layer cohesivecrack model with spring-contact in [4]. A local snap-through instability is observed in the diagram (a) when the upper crack starts to propagate in A and approaches the lower crack tip, in B. Then the load to propagate the crack must be increased, due to a shielding phenomenon, up to point C where the two cracks propagate together unstably. In the diagram (b) the lower crack, which is shorter, starts to propagate at the maximum load, point A; crack propagation is unstable and characterized by a snap-back instability up to point B. Then there is a sudden drop in the load, to point C, caused by a sudden amplification discontinuity. After point C the lower crack continues to propagate alone.

DISCUSSION
he homogenized model uses a zigzag kinematic approximation based on first order shear deformation theory (FSDT) which is enriched by local zigzag functions and interfacial jumps, as shown in Eqn. (1). The low order of the global model and the small number of unknowns, which are those of FSDT, implies that shear strains are constants in and between the layers in the intact regions of the beam and vanish in the delaminated regions. Following the approach which is commonly used for the structural low order theories, accurate shear stresses and strains can be obtained a posteriori through the imposition of local equilibrium. The diagrams in Figs. 3 and 4 show that the complex local fields which arise in multilayered wide plates with imperfect or fully debonded interfaces are accurately reproduced, also in very thick and highly anisotropic cases. However, stress recovery using the equation above is not very accurate in a finite element framework [8] and this limits the applicability of the model. The introduction of a shear correction factor, 44 5 / 6 k  , in the equilibrium Eqns. (6) allows to account for the shear strains in the solution of the problem and the displacements of fully bonded beams are accurately predicted. On the other hand, due to the vanishing of the shear strains in the delaminated regions of the beams, the transverse displacements in Eqn. (3) and the equilibrium Eqns. (6) will only account for the bending contributions and this may affect the load re-distribution between different regions in statically indeterminate systems. This problem cannot be solved by increasing the order of the ESL theory used as global model and, in our opinion, could be solved only by modifying the zigzag formulation and leaving the displacement jumps as part of the unknowns of the problem, as it has been done for instance in [29]. Another approach to solve this problem could be through the use of the refined zigzag theory [7]; however, preliminary analyses in [20,30] highlight other still unresolved difficulties of the approach in dealing with in-plane discontinuities. The effects of neglecting shear deformations along the delaminated portions of the beam, however, are not relevant in all problems dominated by bending defomations (as in the fracture analyses of the cross ply laminated ELS specimen in Fig. 6). Another limitation of the homogenized approach presented here is related to the imposition of continuity and boundary conditions in terms of global quantities, either displacement or force and moment resultants. This implies that continuity is not satisfied at the local level and displacements and force sub-resultants within the single layers may differ. This typically occurs near the crack tip and clamped ends, within boundary regions whose size depends on the mismatch of the elastic constants of the layers and on the stiffness of the interfaces and is very small when the mismatch is small and the interfacial stiffness is very large or very small [20,22]. The model and analyses presented here are limited to mode II dominated problems, were transverse extensibility of the layers and interfacial openings may be neglected. The homogenized approach is currently being extended (work in progress [31]) to general mixed-mode problems using the extended formulation in [14], which uses first order shear, first order normal deformation theory as global model.

CONCLUSIONS
homogenized structural model based on a zigzag approach has been presented for the analysis of beams with imperfect interfaces and delaminations and layers oriented along the geometrical axes. The model enriches the displacement field of first order shear deformation theory, in order to introduce zigzag effects due to the layered architecture and interfacial sliding jumps due to the presence of soft interlayers, imperfections or delaminations. Piece-wise linear cohesive traction laws describe the constitutive behavior of the interfaces. Homogenization and variational techniques are used to define the local variables in terms of the global variables and derive the equilibrium equations of the problem. The model has been validated through various novel applications and using experimental results and 2D elasticity solutions. T A Advantages and limitations of the approach have been discussed. Local stress and displacement fields are accurately predicted also in very thick and highly anisotropic beams and wide plates with continuous imperfect interfaces and delaminations. Crack propagation under mode II dominant conditions has been successfully reproduced in various fracture specimens and bend beams: unidirectionally reinforced End-Notched Flexure and End-Load Split specimens; End-Load Split specimen made of a 46-layers cross-ply laminate; isotropic cantilever beams with multiple delaminations interacting during propagation. The solution is analytical or semi-analytical, requires only in-plane discretization of the problem and the displacement unknowns are only three as in first order shear deformation theory.