Numerical Study of the Mixed-Mode Delamination of Composite Specimens

The present research deals with the delamination process in multi-layered composite specimens, with a reduced computational effort. The adhesive interface between sublaminates is represented as a continuous distribution of elastic-brittle springs in the normal and/or tangential direction depending on the interfacial mixed-mode condition. Each composite adherend, instead, is modelled according to the Timoshenko’s beam theory. The proposed formulation is here enhanced through the Generalized Differential Quadrature (GDQ) method, where the differential equations of the problem are solved directly in a strong form. Thus, the possibility of tracking the delamination response of the specimens is provided locally in a numerical sense, in terms of interface stresses, internal forces and displacements but also in terms of critical fracture energies and mode mixity angles. A further check of the proposed formulation is performed with respect to some standard solutions available in literature. The good agreement between numerical and theoretical predictions verifies the efficiency of the proposed GDQ approach for the study of complex mixed-mode delamination phenomena in composite materials and joints.


Introduction
In a context where a lot of engineering components and material systems are layered and made of high performance composites, an increased attention in the last decades has focused on the development of appropriate numerical tools to simulate the complex damage process known as delamination.This phenomenon can be approached at different scales, from the microscale up to the macroscale.At the micro-level, indeed, the cracks form within the matrix-reach regions, affecting the interfacial weakness of the microstructure.At a macro-level, laminated models are appropriate when the structure is subjected to in-plane, bending or combined loading conditions, as considered in the parametric analysis of the present work.
Different methods based on strength of materials or on fracture mechanics are applied in literature to detect the onset of damage, according to experimental observations, as well as to analytical or numerical approaches.In the first case, several test configurations have been developed to measure the delamination strength and/or toughness, see for example, [1][2][3][4], among others.The mixed-mode delamination can either be modelled through analytical approaches, accounting for possible geometrical and mechanical non-linearities of materials and interfaces.
In the pioneering work by Williams [5], the beam theory was applied to determine the mixed-mode energy release rate G and its pure-mode contributions G I and G I I , for homogeneous beams with arbitrary crack depths.Suo and Hutchinson proposed a Local Method (LM) to determine the mode-mixity associated to the growth of a semi-infinite crack between two homogeneous and isotropic layers [6] or within an infinite orthotropic elastic strip [7].Further studies include the effects of shear forces on the delamination process of isotropic or orthotropic materials, as proposed by Li et al. [8] and Andrews and Massabò [9].Subsequent elastic interface models, based on a continuous distribution of linear springs with appropriate stiffness parameters, where proposed in [10][11][12][13][14][15] to study the delamination process of homogeneous or bimaterial beams.More specifically, Bruno and Greco [10] proposed a simple approach in which the straight delamination growth was studied between layers, modelled as Kirchhoff or Reissner-Mindlin plates, while giving the closed form solutions for the energy release rate in the limited case of rigid interface.A comparative evaluation between rigid, semi-rigid and flexible interface models was also proposed by Qiao and Wang [11], who provided the closed-form solutions of fracture parameters, useful for practical applications.Alfredsson and Högberg [12] introduced a theoretical model for studying the mixed-mode behaviour of adhesive joints, including the Euler-Bernoulli beam hypotheses for the adherends and the elastic connection for the adhesive layer.The algebraic expression for the mode mixity was also presented by the same authors in the limit case of semi-infinite asymmetric double cantilever beam (ADCB) test.An improved version of the Timoshenko's beam theory, named as Enhanced Beam Theory (EBT), was applied by Bennati et al. [13][14][15] to study the standard ADCB test [13], as well as the mixed-mode bending (MMB) test [14,15].Based on the EBT model, the adherends were modelled as extensible, flexible and shear-deformable beams, whereas the adhesive interface was regarded as a continuous distribution of linearly elastic-brittle springs acting along the normal and tangential directions with respect to the interface plane.The same model has been recently extended in [16,17] to laminated orthotropic beams, including the effect of general stacking sequences, the bending-extension coupling and the shear deformability of the adherends.Starting with the main idea provided recently by Valvo [16] and Dimitri et al. [17], we apply the EBT to study the mixed-mode delamination of composite specimens, made by extensible, flexible and shear-deformable adherends, partly jointed by a deformable interface, under different axial, shear and bending loading conditions.The EBT is solved, first, in a closed form, for the simplest cases, while selecting the interface stresses as the main unknowns.This generalized formulation could be applied for the design of non-standard delamination toughness tests and may include different possible combinations of standard specimens, such as the symmetric or asymmetric DCBs, the end load split test (ELS), a peel test, as well as the same moment-loaded DCB (MLDCB), among others.This would enable a full account of mixed-mode effects, both on the static, kinematic, strain and energy response.
The same problem is further solved numerically via the Generalized Differential Quadrature (GDQ) method, in line with the approach proposed in the recent work by Dimitri et al. [17].Thus, the main solution of the problem is found in terms of generalized displacements of the sublaminates, internal forces, local interface stresses, mode-mixity angle and critical energy.The feasibility of the EBT for generalized mixed-mode conditions is checked against the LM, as applied by Hutchinson and Suo [18] in a closed-form.At the same time, the accuracy of the GDQ approach is verified with a preliminary convergence analysis of the numerical solutions with respect to the analytical ones, with excellent results also for a reduced computational cost of the problem.In the recent years, it has been demonstrated that the GDQ is a fast and reliable approach for solving in-plane problems with discontinuities, for which a very fine mesh would be required to capture well the results around discontinuities (see e.g., the parametric investigations in [19,20] compared to the Finite Element Method (FEM), eXtended FEM (XFEM), or isogeometric approximations).The GDQ solves the strong form of the mathematical problem and needs to enforce the boundary conditions a posteriori, whereas FEM solves the weak problem and set a priori the boundary conditions.In this framework, we want to apply the GDQ approach to delamination problems in mixed-mode conditions, as recently done in the companion work [17], for a moment-loaded DCB and here extended to include different possible combinations of standard specimens.A further development of the proposed model will account for the presence of a softening stage in the stress-separation laws, as typically considered for cohesive crack models and recently applied by the authors in [21,22].
The manuscript is organized as follows: in Section 2 the analytical formulation of the problem based on the concepts of the EBT is provided.Section 3 describes a possible analytical solution strategy, while Section 4 focuses on the main fundamentals of the GDQ approach, here suggested to solve also complex cases.An extended parametric analysis is addressed in Section 5, whereas the main conclusions and possible developments are drawn in Section 6.

Mechanical Problem
Consider a layered beam, with total length L, precrack length a, such that b = L − a corresponds to the uncracked length.The specimen is made of two orthotropic sublaminates, here labelled as 1 and 2, with layers made of different materials and thicknesses, as well as fibre-reinforced laminae with different orientations.General stacking sequences are allowed by the proposed model, including those with bending-extension coupling but neglecting the out-of-plane effects (namely the torsion, the out-of-plane shear, etc.).The sublaminate thicknesses are indicated as H 1 and H 2 , respectively, whereas the global thickness of the specimen is H = H 1 + H 2 .The width of the specimen is generally defined as t 1 and t 2 for sublaminates 1 and 2, respectively, such that t = t 1 = t 2 .Each arm i of the specimen is loaded at the left side with an external normal force N i , shear force T i , or bending moment M i (see Figure 1), in quasi-static conditions, such that the dimensionless loading ratios α N = N 1 /N 2 , α T = T 1 /T 2 and α M = M 1 /M 2 can be suitably varied from −1 to 1 in order to embrace the whole range of mixed-modes between pure-modes I and II.The manuscript is organized as follows: in Section 2 the analytical formulation of the problem based on the concepts of the EBT is provided.Section 3 describes a possible analytical solution strategy, while Section 4 focuses on the main fundamentals of the GDQ approach, here suggested to solve also complex cases.An extended parametric analysis is addressed in Section 5, whereas the main conclusions and possible developments are drawn in Section 6.

Mechanical Problem
Consider a layered beam, with total length L , precrack length a , such that b L a = − corresponds to the uncracked length.The specimen is made of two orthotropic sublaminates, here labelled as 1 and 2, with layers made of different materials and thicknesses, as well as fibre-reinforced laminae with different orientations.General stacking sequences are allowed by the proposed model, including those with bending-extension coupling but neglecting the out-of-plane effects (namely the torsion, the out-of-plane shear, etc.).The sublaminate thicknesses are indicated as t t t = = .Each arm i of the specimen is loaded at the left side with an external normal force i N , shear force i T , or bending moment i M (see Figure 1), in quasi-static conditions, such that the dimensionless loading ratios  The two sublaminates are here modelled as Timoshenko beams, partly connected by an elastic-brittle interface (see Figure 1).Each sublaminate i is constituted by i l laminae, with the th k layer enclosed generally within the ( ) i k z + coordinates from the centreline, such that its thickness is ( ) . This means that the total thickness for each arm is equal to The two sublaminates are here modelled as Timoshenko beams, partly connected by an elasticbrittle interface (see Figure 1).Each sublaminate i is constituted by l i laminae, with the kth layer enclosed generally within the z i(k) and z i(k+1) coordinates from the centreline, such that its thickness is h i(k) = z i(k+1) − z i(k) .This means that the total thickness for each arm is equal to as visible in Figure 1.
The coordinate s denotes the general distance of cross sections from the crack tip along the longitudinal direction.Thus, O i x i z i stands for the local reference system for each sublaminate centred at its mid-plane.Accordingly, we define the mid-plane displacements u i , w i for each sublaminate, along the axial and transverse directions, respectively and we denote the cross-section rotation φ i , positive if counter-clockwise (Figure 1).
Based on classical lamination theory, the sublaminates are here modelled as homogeneous equivalent beams, while introducing the equivalent stiffnesses, 44 , to define their extensional stiffness, bending-extension coupling stiffness, bending stiffness and shear stiffness, respectively, as follows: where κ stands for the shear correction factor, here assumed equal to 5/6.The equilibrium equations for sublaminates in the unbroken part (i.e.,s N i , T i , M i being, respectively, the axial and shear forces and bending moments within sublaminates.These actions are strictly related to the normal and tangential interface stresses, defined as: where ∆w = w − 2 − w + 1 , ∆u = u − 2 − u + 1 are the relative interface displacements in the transverse and axial directions, respectively, whereas w + 1 , u + 1 are the kinetic quantities at the bottom surface of sublaminate 1 (z 1 = H 1 /2) and w − 2 , u − 2 are the kinetic quantities at the top surface of sublaminate 2 (z 2 = −H 2 /2).The mechanical formulation proposed hereafter, is based on the Timoshenko's beam theory, according to which the transverse displacement maintains constant along the thickness (i.e.,w + 1 = w 1 ,w − 2 = w 2 ) and the axial displacement varies linearly with the thickness coordinate (i.e., u ).Hence, the relative separations are The constitutive laws for each arm can be expressed as where By combining Equations ( 5), ( 6) and ( 2) we obtain the following set of differential equations which is associated to the following boundary conditions (BCs) These last relations can be expressed in terms of the kinematic quantities u i , w i , φ i and their derivatives, by adopting Equations ( 5) and (6).Please, note that Mi = M i + aT i in Equation (8).
The differential problem ( 7) is solved, first, in closed form, for homogeneous, isotropic and unidirectional specimens, after a suitable change of variables, as proposed by Bennati et al. [13] for a similar problem.More in detail, the interfacial stresses are selected as primary unknowns for the analytical solution of the problem.The same is then solved numerically in a discretized GDQ setting, which enables a direct solution of Equation (7) in the kinematic unknowns with a reduced computational effort.Once the accuracy and stability of the proposed GDQ approach has been checked, we perform a parametric study of the delamination response in a static, kinematic and energy sense, for isotropic and orthotropic specimens under different mixed-mode conditions.

Analytical Solution Strategy
The closed-form solution of the differential problem (7) is not straightforward and it is here found only for the specimen with homogeneous, isotropic and unidirectional sublaminates.In this case, the sublaminate stiffnesses

get as follows
where E i is the longitudinal Young's modulus, G i is the shear modulus and I i = t i H 3 i /12 is the second moment of inertia for each sublaminate.Thus, the problem is governed by the following differential equations which can be redefined as function of the interface stresses, by combining Equations ( 3) and (4).
After some mathematical manipulation, we obtain the following differential equation for the normal interface stress and The general solution of Equation ( 11) reads where constants C 1 , C 2 , ..., C 6 are computed by enforcing the BCs (8), whereas λ 1 , λ 2 , ..., λ 6 are the roots of the characteristic equation associated to Equation ( 14).In addition, the tangential stress τ is determined by solving the following 1st order differential equation with σ defined according to Equation (14).By integration of Equation ( 15), the tangential stress is found to be equal to where an additional constant C 7 must be determined.By substitution of Equations ( 14) and ( 16) into Equation ( 2) and by integration, we determine the explicit expressions for the internal forces, namely for the axial forces, (18) for the shear forces and finally for the bending moments.With a similar procedure, we can determine the explicit expressions for the kinematic quantities of each arm.More specifically, the substitution of the internal forces ( 17)- (19) into Equations ( 5) and ( 6) and the integration of the final relations with respect to s, yields to the expressions below for the axial displacements, for the transverse displacements and lastly for rotations.
It is worth noticing that the static and kinematic relations in Equations ( 17)-( 22) depend on 12 additional constants C 8 , C 9 , ..., C 19 , with a total number of constants that has to be determined.In addition to the 12 BCs defined by the relations (8), other 7 relations are required to determine all the unknown constants.These remaining relations are determined by substituting Equations ( 14), ( 16), ( 20)- (22) into Equation (3).More specifically, the following relations can be applied to evaluate the first six constants independently, namely with âi = 1/(E i I i ), for (i = 1, 2), whereas the following further constants have been introduced Upon mathematical manipulation, the remaining 13 constants are computed as a function of the first six ones, namely After defining the analytical expressions for the kinematic quantities, we can determine the energy release rate as in the following where σ c , τ c denote the normal and tangential interface stresses at the crack tip, defined as Moreover, we can compute explicitly the mode-mixity angle, as shown below which is here adopted in the crack-growth criterion of the type [18] in order to control the delamination stage of the specimen.

Numerical Strategy: The GDQ Approach
In this section, we show the main fundamentals about the GDQ approach, as here applied to solve numerically the system of Equation ( 7) in a discrete form.Based on the GDQ approach, the nth order of a derivative at a coordinate s = s i can be computed as a weighted linear sum of the functional values collocation points as follows where N is the total number of collocation points for the discretization of domain along the s−directions, whereas ξ (n) ij denotes the weighting coefficients determined with the formulae by Shu and Richards [23,24], for i, j = 1, 2, ..., N and n = 1, 2, ..., N − 1 ij stands for the weighting coefficients for the first-order derivatives defined below A key aspect of the GDQ approach is dictated by its efficient and accurate application, which is in turn related to the appropriate selection of grip points within the domain.Among different possible choices of grid point distributions, non-uniform discretizations represent an efficient way to ensure more accurate results compared to the uniform option, in agreement with findings by Shu [25] and some comparative evaluations, as found in [26,27].In this context, a preliminary convergence analysis is performed for varying grid distributions, in terms of the traction vector t along the adhesive interface.This is defined by its normal and tangential components σ, τ, respectively, as follows More in detail, the following grid distributions are evaluated comparatively, that is, Che-Gau-Lob, Chebyshev I (Cheb I), Chebyshev II (Cheb II), Chebyshev III (Cheb III), Chebyshev IV (Cheb IV), Legendre (Leg), Chebyshev-Gauss (Cheb-Gau), Legendre-Gauss (Leg-Gau), Lobatto (Lob), Lobatto-Gauss (Lob-Gau), Legendre-Gauss-Lobatto (Leg-Gau-Lob), Chebyshev-Gauss-Radau (Che-Gau-Rad), Legendre-Gauss-Radau (Leg-Gau-Rad), Jacobi (Jac), Jacobi-Gauss (Jac-Gau), see [28] for more details.
As visible in the curves of Figure 2, the rate of convergence is very fast and almost unaffected by the selected discretization.The method seems very stable and limit the error to 10 −6 ÷ 10 −8 with a low number of grid points ranging between N = 75 and N = 101.It follows a round-off error phenomenon for higher discretizations, as already observed in the companion work [17].Based on the results, a Che-Gau-Lob interpolation with a total number of collocation points N = 81 will be selected, hereafter, for the numerical investigation.Thus, each collocation point along the interface is defined by its coordinate s k from the crack tip, as in the following L error norm of the traction vector at the interface for varying grid distributions.

Numerical Examples
This section is devoted to many illustrative examples aimed at investigating the mixed-mode structural response for varying loading, geometry and material conditions, either in an independent or a combined form.Comparison with existing formulations from literature helps to highlight the capability of the proposed mixed-mode EBT to capture the fracturing response.An axial force 2 500 N N = is applied on the left side of the specimen, while embracing different possible mixed-modes by varying N α from 1 (pure mode-I) to −1 (pure mode-II).In Figure 3 we plot the structural response in terms of normal and tangential stresses σ and τ along the coordinate s of the specimen.As visible in Figure 3, the peel stress σ in the normal direction maintains always null for all mixed-modes, whereas the shear stress τ reaches the maximum value at the crack tip and decays within a short distance of about 5 mm with a monotonic behaviour.As also expected, the magnitude of the peak values of τ , reduces gradually moving from pure mode-II to pure mode-I.Only in the last case, τ vanishes along the whole specimen, which means that the adhesive is completely unloaded.

Numerical Examples
This section is devoted to many illustrative examples aimed at investigating the mixed-mode structural response for varying loading, geometry and material conditions, either in an independent or a combined form.Comparison with existing formulations from literature helps to highlight the capability of the proposed mixed-mode EBT to capture the fracturing response.

Loading Effect
A first parametric investigation aims at analysing the mixed-mode response induced by a varying loading condition for a symmetric specimen with length L = 150 mm, precrack length a = 35 mm, width t i = 20 mm, thickness H i = 2.5 mm(i = 1, 2).Each sublaminate is made of a 16-ply unidirectional E-glass/epoxy material with elastic moduli E i = 25.7 GPa and G i = 2.5 GPa.
An axial force N 2 = 500 N is applied on the left side of the specimen, while embracing different possible mixed-modes by varying α N from 1 (pure mode-I) to −1 (pure mode-II).In Figure 3 we plot the structural response in terms of normal and tangential stresses σ and τ along the coordinate s of the specimen.As visible in Figure 3, the peel stress σ in the normal direction maintains always null for all mixed-modes, whereas the shear stress τ reaches the maximum value at the crack tip and decays within a short distance of about 5 mm with a monotonic behaviour.As also expected, the magnitude of the peak values of τ, reduces gradually moving from pure mode-II to pure mode-I.Only in the last case, τ vanishes along the whole specimen, which means that the adhesive is completely unloaded.As far as the structural behaviour of the sublaminates is concerned, Figures 4-6 present the parametric results in terms of internal forces, kinematic and strain quantities, respectively.More in detail, the axial forces in the two sublaminates 1 2 , N N are perfectly the same and maintain constant in the whole specimen when 1 N α = (pure mode-I) and assume a symmetric behaviour (i.e., equal in magnitude and opposite in sign) when 1 N α = − (pure mode-II).In mixed-mode conditions (i.e., , N N feature a monotonic behaviour with a decreasing or increasing trend up to a threshold value within a short distance from the crack tip (see Figure 4a).As also visible in Figure 4b, the internal shear forces 1 2 , T T are null everywhere within the specimen, while the bending moments 1 2 , M M are always coincident and different from zero, for both sublaminates, except for pure mode-I for which 0 τ = and do not trigger any bending moment within the specimen (Figure 4c).
A set of parametric investigations is also repeated in a kinematic sense, as visible in the results plotted in Figure 5.More specifically, Figure 5a represents the axial displacements 1 2 , u u as computed at the centrelines of sublaminates.These kinematic quantities are always different from zero and vary almost linearly along each sublaminate between a maximum value (for 0 s = ) and zero at the clamped side of the specimen (for ).In Figure 5b,c, the transverse displacements 1 2 , w w and rotations 1 2 , φ φ are exactly the same in magnitude and sign for all mixed-modes, with the maximum deflection reached at the crack tip of the specimen.By means of the kinematic relations ( 16), we compute the axial strain 1 2 , ε ε , shear strain 1 2 , γ γ and curvature , χ χ of the two sublaminates, whose results are depicted in Figure 6.As expected, for an isotropic specimen under an axial loading, the only deformations involved are the axial and flexural ones.
In mixed mode I/II fracture problem, it is well known that the mode-mixity can be characterized by the phase angle of the complex stress-intensity factor.Thus, we compute the mode-mixity angle ψ according to Equation (45) and plot this quantity against the mixed-mode loading ratios N α (see Figure 7a).A comparative evaluation of the results is also provided with respect to the LM, as proposed in literature by Suo and Hutchinson [18].According to this last formulation, the shear deformability of sublaminates is not taken into account and does not include any possible dependence of the mode-mixity upon the delamination length a .As drawn in Figure 7a, the mode-mixity angle ψ seems to be unaffected by the selected N α and maintains constant for all mixed-modes (at least for an isotropic specimen, as studied in this first example).The EBT and the LM, in addition, give almost the same results, which proves the validity of the proposed numerical formulation to solve a fracturing problem.As far as the structural behaviour of the sublaminates is concerned, Figures 4-6 present the parametric results in terms of internal forces, kinematic and strain quantities, respectively.More in detail, the axial forces in the two sublaminates N 1 , N 2 are perfectly the same and maintain constant in the whole specimen when α N = 1 (pure mode-I) and assume a symmetric behaviour (i.e., equal in magnitude and opposite in sign) when α N = −1 (pure mode-II).In mixed-mode conditions (i.e., −1 < α N < 1) N 1 , N 2 feature a monotonic behaviour with a decreasing or increasing trend up to a threshold value within a short distance from the crack tip (see Figure 4a).As also visible in Figure 4b, the internal shear forces T 1 , T 2 are null everywhere within the specimen, while the bending moments M 1 , M 2 are always coincident and different from zero, for both sublaminates, except for pure mode-I for which τ = 0 and do not trigger any bending moment within the specimen (Figure 4c).
A set of parametric investigations is also repeated in a kinematic sense, as visible in the results plotted in Figure 5.More specifically, Figure 5a represents the axial displacements u 1 , u 2 as computed at the centrelines of sublaminates.These kinematic quantities are always different from zero and vary almost linearly along each sublaminate between a maximum value (for s = 0) and zero at the clamped side of the specimen (for s = b = 115 mm).In Figure 5b,c, the transverse displacements w 1 , w 2 and rotations φ 1 , φ 2 are exactly the same in magnitude and sign for all mixed-modes, with the maximum deflection reached at the crack tip of the specimen.By means of the kinematic relations ( 16), we compute the axial strain ε 1 , ε 2 , shear strain γ 1 , γ 2 and curvature χ 1 , χ 2 of the two sublaminates, whose results are depicted in Figure 6.As expected, for an isotropic specimen under an axial loading, the only deformations involved are the axial and flexural ones.
In mixed mode I/II fracture problem, it is well known that the mode-mixity can be characterized by the phase angle of the complex stress-intensity factor.Thus, we compute the mode-mixity angle ψ according to Equation (45) and plot this quantity against the mixed-mode loading ratios α N (see Figure 7a).A comparative evaluation of the results is also provided with respect to the LM, as proposed in literature by Suo and Hutchinson [18].According to this last formulation, the shear deformability of sublaminates is not taken into account and does not include any possible dependence of the mode-mixity upon the delamination length a.As drawn in Figure 7a, the mode-mixity angle ψ seems to be unaffected by the selected α N and maintains constant for all mixed-modes (at least for an isotropic specimen, as studied in this first example).The EBT and the LM, in addition, give almost the same results, which proves the validity of the proposed numerical formulation to solve a fracturing problem.Similar considerations can be also repeated in terms of critical energy which dictates the starting point for the delamination process.Also in this case, the EBT agrees very well with the LM, whose estimations are always lower than predictions based on the EBT.This would have a direct consequence on the global load-displacement response, as will be investigated in a further work.
As second example, we consider a more complex orthotropic specimen, with the following sequence of laminae for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness 1 0.75 mm h = , 2 1.3 mm h = , 3 0.45 mm h = (and total thickness 1 2 3
The Graphite-Epoxy material features a Young's modulus 137.9 GPa and shear modulus ( ,1) ( ,3) 13 13 7.1 GPa , whereas the Glass-Epoxy material has a Young's modulus Although a similar problem would be cumbersome to be solved in a closed form, the numerical solution provides an easy and efficient strategy to check for the local response of the specimen, for all mixed-modes.A parametric investigation is then repeated in the static and kinematic sense, as shown in Figures 8-11.As visible in these figures, the presence of bending-extension coupling stiffness ( ) 11 i B , for an orthotropic specimen, will always allow to a mixed-mode behaviour both at the interface (see the interface stresses in Figure 8) and within each sublaminate (see the internal forces, displacements and strain in Figures 9-11), also for pure-mode loading conditions (namely for 1  Similar considerations can be also repeated in terms of critical energy which dictates the starting point for the delamination process.Also in this case, the EBT agrees very well with the LM, whose estimations are always lower than predictions based on the EBT.This would have a direct consequence on the global load-displacement response, as will be investigated in a further work. As second example, we consider a more complex orthotropic specimen, with the following sequence of laminae for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness h 1 = 0.75 mm, h 2 = 1.3 mm, h 3 = 0.45 mm (and total thickness Although a similar problem would be cumbersome to be solved in a closed form, the numerical solution provides an easy and efficient strategy to check for the local response of the specimen, for all mixed-modes.A parametric investigation is then repeated in the static and kinematic sense, as shown in Figures 8-11.As visible in these figures, the presence of bending-extension coupling stiffness B (i) 11 , for an orthotropic specimen, will always allow to a mixed-mode behaviour both at the interface (see the interface stresses in Figure 8) and within each sublaminate (see the internal forces, displacements and strain in Figures 9-11), also for pure-mode loading conditions (namely for α N = ±1).Similar considerations can be also repeated in terms of critical energy which dictates the starting point for the delamination process.Also in this case, the EBT agrees very well with the LM, whose estimations are always lower than predictions based on the EBT.This would have a direct consequence on the global load-displacement response, as will be investigated in a further work.

mm). The Graphite-Epoxy material features a Young's modulus E
As second example, we consider a more complex orthotropic specimen, with the following sequence of laminae for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness 1 0.75 mm h = , 2 1.3 mm h = , 3 0.45 mm h = (and total thickness 1 2 3
The Graphite-Epoxy material features a Young's modulus 137.9 GPa , whereas the Glass-Epoxy material has a Young's modulus .
Although a similar problem would be cumbersome to be solved in a closed form, the numerical solution provides an easy and efficient strategy to check for the local response of the specimen, for all mixed-modes.A parametric investigation is then repeated in the static and kinematic sense, as shown in Figures 8-11.As visible in these figures, the presence of bending-extension coupling stiffness ( ) 11 i B , for an orthotropic specimen, will always allow to a mixed-mode behaviour both at the interface (see the interface stresses in Figure 8) and within each sublaminate (see the internal forces, displacements and strain in Figures 9-11), also for pure-mode loading conditions (namely for 1

Geometry Effect
Another possible way of studying the mixed-mode delamination is related to the geometry of the specimen, as considered in the foregoing by varying the thickness ratio η = H 1 /H 2 while assuming a mode-II loading condition with N 2 = 500 N and α N = −1.The asymmetric DCB has length L = 150 mm, precrack length a = 35 mm, width t i = 20 mm and has two sublaminates made of 16-ply unidirectional E-glass/epoxy isotropic material with elastic moduli E i = 25.7 GPa and G i = 2.5 GPa.In Figure 12, the numerical solutions in terms of interface stresses σ and τ show the high sensitivity to the mixed-mode thickness ratio.More specifically, increased values of the thickness ratio η < 1 make the oscillating response of σ and τ less pronounced.This is due to the increase in stiffness of sublaminate 1 with respect to sublaminate 2. The pure mode-II condition is obtained for a perfectly symmetric specimen, that is, for η = 1, which corresponds to the presence of the only tangential stress τ and a complete absence of the normal component σ.An increased value of η > 1 yields to an inversion in sign of σ (Figure 12a) and to a negligible reduction of τ(Figure 12b).The consistency of the formulation is also investigated in terms of internal forces N, T, M and displacements of the sample, as represented in Figures 13 and 14, respectively.As clearly shown in Figure 13a,b, the axial and shear forces N and T keep always symmetric, independently of the thickness ratio η.Bending moments, instead, exhibit always an asymmetric behaviour, except for η = 1 when a pure mode-II condition is reached (see Figure 13c).A similar consideration can be repeated by analysing the kinematic and strain distributions within sublaminates (see Figures 14 and 15).Moreover, all the displacement components in Figure 14, decay to negligible values within a short distance from the crack tip, in agreement with findings by Bennati et al. [13] and by Dimitri et al. [17].A further consistency check aims at computing the mixed-mode angle ψ and the critical energy G c for varying thickness ratios (Figure 16).
the specimen, as considered in the foregoing by varying the thickness ratio .In Figure 12, the numerical solutions in terms of interface stresses σ and τ show the high sensitivity to the mixed-mode thickness ratio.More specifically, increased values of the thickness ratio 1 η < make the oscillating response of σ and τ less pronounced.This is due to the increase in stiffness of sublaminate 1 with respect to sublaminate 2. The pure mode-II condition is obtained for a perfectly symmetric specimen, that is, for 1 η = , which corresponds to the presence of the only tangential stress τ and a complete absence of the normal component σ .An increased value of 1 η > yields to an inversion in sign of σ (Figure 12a) and to a negligible reduction of τ (Figure 12b).The consistency of the formulation is also investigated in terms of internal forces , , N T M and displacements of the sample, as represented in Figures 13 and 14, respectively.As clearly shown in Figure 13a,b, the axial and shear forces N and T keep always symmetric, independently of the thickness ratio η .Bending moments, instead, exhibit always an asymmetric behaviour, except for 1 η = when a pure mode-II condition is reached (see Figure 13c).A similar consideration can be repeated by analysing the kinematic and strain distributions within sublaminates (see Figures 14 and 15).Moreover, all the displacement components in Figure 14, decay to negligible values within a short distance from the crack tip, in agreement with findings by Bennati et al. [13] and by Dimitri et al. [17].A further consistency check aims at computing the mixed-mode angle ψ and the critical energy c G for varying thickness ratios (Figure 16).
The EBT-based numerical results are compared to the theoretical predictions according the LM, with a satisfactory correspondence between them, at least within a certain range of thickness ratios η .A non-monotonic variation of both ψ and cr G is observed when applying the EBT approach, instead of an increasing and monotonic trend, as expected by the LM.A perfect correspondence between the two approaches is also reached in pure mode-II condition, when ψ becomes equal to 90° (Figure 16a) and cr G reaches the limit value 1.68 (Figure 16b).(Isotropic specimen).

Mechanical Effect
In this section, we analyse the possible effect of the mechanical properties of the specimen on its mixed-mode behaviour, for a fixed pure mode-II loading condition.The elastic properties of sublaminate 1 are related to those ones of sublaminate 2, here set to 2 25.7 GPa E = and 2 2.5 GPa G = .
Thus, a dimensionless mechanical parameter is introduced , which is varied from 0. β < corresponds to an increasing stiffness of sublaminate 1, which tends to the one of sublaminate 2 in the limit case given by the condition 1 β = .
Based on the plot in Figure 17, the oscillating behaviour of σ and τ features a reduced magnitude for increasing values of β .This reverts to a pure condition II when  The EBT-based numerical results are compared to the theoretical predictions according the LM, with a satisfactory correspondence between them, at least within a certain range of thickness ratios η.A non-monotonic variation of both ψ and G cr is observed when applying the EBT approach, instead of an increasing and monotonic trend, as expected by the LM.A perfect correspondence between the two approaches is also reached in pure mode-II condition, when ψ becomes equal to 90 • (Figure 16a) and G cr reaches the limit value 1.68 (Figure 16b).

Mechanical Effect
In this section, we analyse the possible effect of the mechanical properties of the specimen on its mixed-mode behaviour, for a fixed pure mode-II loading condition.The elastic properties of sublaminate 1 are related to those ones of sublaminate 2, here set to E 2 = 25.7 GPa and G 2 = 2.5 GPa.Thus, a dimensionless mechanical parameter is introduced β = E 1 /E 2 = G 1 /G 2 , which is varied from 0.2 to 1, to embrace different possible mixed-modes.A pure mode-II axial loading condition is assumed once again, such that N 2 = −N 1 = 500 N.The local response is plotted in Figure 17 in terms of normal and tangential adhesive stresses σ and τ versus the distance s from the crack tip.It is worth noticing that an increasing value of β < 1 corresponds to an increasing stiffness of sublaminate 1, which tends to the one of sublaminate 2 in the limit case given by the condition β = 1.Based on the plot in Figure 17, the oscillating behaviour of σ and τ features a reduced magnitude for increasing values of β.This reverts to a pure condition II when β = 1, for which the normal stress σ becomes zero and only the tangential stress survives along the specimen.Moreover, the outcomes of this parametric study show the variation of the mechanical response of sublaminates for varying mechanical ratios β = 1, as clearly represented in Figure 18 in a static sense, as well as in Figures 19 and 20 in a kinematic and strain sense, or, finally in Figure 21 from an energy standpoint (namely G cr and ψ).More specifically, the plots in Figure 18 show a pronounced sensitivity of the response in terms of shear forces and bending moments and a meaningless variation of the axial forces for varying stiffnesses.G and ψ ).More specifically, the plots in Figure 18 show a pronounced sensitivity of the response in terms of shear forces and bending moments and a meaningless variation of the axial forces for varying stiffnesses.Also, the kinematic and strain response of the specimen is significantly affected by its mechanical properties (see Figures 19 and 20), with a higher deformability of sublaminate 1 compared to sublaminate 2 (Figure 20a-c).Finally, as far as the energy is concerned, the predictions based on the EBT are always more conservative than those ones given by the LM and tend to them when approaching the pure mode condition (see Figure 21a,b).

Composite Specimen
As last example, we study the mechanical response of the specimen in Section 4, with the same loading conditions but for a more complex composite material.The following sequence of laminae is assumed for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness   Also, the kinematic and strain response of the specimen is significantly affected by its mechanical properties (see Figures 19 and 20), with a higher deformability of sublaminate 1 compared to sublaminate 2 (Figure 20a-c).Finally, as far as the energy is concerned, the predictions based on the EBT are always more conservative than those ones given by the LM and tend to them when approaching the pure mode condition (see Figure 21a,b).

Composite Specimen
As last example, we study the mechanical response of the specimen in Section 4, with the same loading conditions but for a more complex composite material.The following sequence of laminae is assumed for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness h 1 = 0.75 mm, h 2 = 1.3 mm, h 3 = 0.45 mm and η = H 1 /H 2 = 1.
Graphite-Epoxy material features a Young's modulus E = 8.96 GPa.This specimen is here selected with the only purpose of demonstrating the great potential and flexibility of the GDQ to solve numerically more complex problems, otherwise cumbersome to be solved in closed form (look at the plots in Figures 22-25 for the mixed-mode responses of the interface and sublaminates).Also, the kinematic and strain response of the specimen is significantly affected by its mechanical properties (see Figures 19 and 20), with a higher deformability of sublaminate 1 compared to sublaminate 2 (Figure 20a-c).Finally, as far as the energy is concerned, the predictions based on the EBT are always more conservative than those ones given by the LM and tend to them when approaching the pure mode condition (see Figure 21a,b).

Composite Specimen
As last example, we study the mechanical response of the specimen in Section 4, with the same loading conditions but for a more complex composite material.The following sequence of laminae is assumed for each sublaminate: Graphite-Epoxy/Glass-Epoxy/Graphite-Epoxy, of thickness

Conclusions
In this work, we have proposed the EBT to deal with the delamination behaviour of isotropic and composite specimens under different loading, geometry and mechanical mixed-mode conditions.The specimen is modelled as an assemblage of two sublaminates partly precracked and partly bonded by an elastic interface, here defined through a continuous distribution of varied from −1 to 1 in order to embrace the whole range of mixed-modes between pure-modes I and II.

Figure 1 .
Figure 1.Mechanical scheme of the laminated composite specimen.

Figure 1 .
Figure 1.Mechanical scheme of the laminated composite specimen.

A
first parametric investigation aims at analysing the mixed-mode response induced by a varying loading condition for a symmetric specimen with length 150 is made of a 16-ply unidirectional E-glass/epoxy material with elastic moduli 25

Figure 2 .
Figure 2. L 2 error norm of the traction vector at the interface for varying grid distributions.

Figure 3 .
Figure 3.Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 3 .
Figure 3.Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 4 .
Figure 4. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).Figure 4. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

Figure 4 .
Figure 4. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).Figure 4. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

Figure 7 .
Figure 7. Effect of the axial loading ratio on the (a) Mode-mixity angle ψ ; (b) Critical energy cr G .

Figure 8 .
Figure 8.Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Orthotropic specimen).

JFigure 7 .
Figure 7. Effect of the axial loading ratio on the (a) Mode-mixity angle ψ ; (b) Critical energy cr G .

Figure 8 .
Figure 8.Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Orthotropic specimen).Figure 8. Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Orthotropic specimen).

Figure 8 .
Figure 8.Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Orthotropic specimen).Figure 8. Effect of the axial loading on the interface local response: (a) Normal stress; (b) Tangential stress.(Orthotropic specimen).

Figure 9 .
Figure 9.Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Orthotropic specimen).Figure 9. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Orthotropic specimen).

Figure 9 .
Figure 9.Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Orthotropic specimen).Figure 9. Effect of the axial loading on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Orthotropic specimen).
ply unidirectional E-glass/epoxy isotropic material with elastic moduli 25

Figure 12 .
Figure 12.Effect of the thickness ratio η on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 12 .
Figure 12.Effect of the thickness ratio η on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 13 .
Figure 13.Effect of the thickness ratio η on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

Figure 13 .
Figure 13.Effect of the thickness ratio η on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

Figure 16 .
Figure 16.Effect of the thickness ratio η on the (a) Mode-mixity angle ψ ; (b) Critical energy cr G .

2 to 1 ,
to embrace different possible mixed-modes.A pure mode-II axial loading condition is assumed once again, such that .The local response is plotted in Figure17in terms of normal and tangential adhesive stresses σ and τ versus the distance s from the crack tip.It is worth noticing that an increasing value of 1

1 β=Figure 17 .
Figures 19 and 20 in a kinematic and strain sense, or, finally in Figure21from an energy standpoint (namely cr G and ψ ).More specifically, the plots in Figure18show a pronounced sensitivity of the response in terms of shear forces and bending moments and a meaningless variation of the axial forces for varying stiffnesses.

Figures 19 and 20
Figures 19 and 20 in a kinematic and strain sense, or, finally in Figure 21 from an energy standpoint (namely cr

Figure 17 .
Figure 17.Effect of the mechanical ratio β on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 17 .
Figure 17.Effect of the mechanical ratio β on the interface local response: (a) Normal stress; (b) Tangential stress.(Isotropic specimen).

Figure 18 .
Figure 18.Effect of the mechanical ratio β on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

Figure 18 .
Figure 18.Effect of the mechanical ratio β on the static response of the sublaminates: (a) Axial force; (b) Shear force; (c) Bending moment.(Isotropic specimen).

.
This specimen is here selected with the only purpose of demonstrating the great potential and flexibility of the GDQ to solve numerically more complex problems, otherwise cumbersome to be solved in closed form (look at the plots in Figures22-25for the mixed-mode responses of the interface and sublaminates).

.
This specimen is here selected with the only purpose of demonstrating the great potential and flexibility of the GDQ to solve numerically more complex problems, otherwise cumbersome to be solved in closed form (look at the plots in Figures22-25for the mixed-mode responses of the interface and sublaminates).