An Elastic Interface Model for the Delamination of Bending-Extension Coupled Laminates

: The paper addresses the problem of an interfacial crack in a multi-directional laminated beam with possible bending-extension coupling. A crack-tip element is considered as an assemblage of two sublaminates connected by an elastic-brittle interface of negligible thickness. Each sublaminate is modeled as an extensible, ﬂexible, and shear-deformable laminated beam. The mathematical problem is reduced to a set of two differential equations in the interfacial stresses. Explicit expressions are derived for the internal forces, strain measures, and generalized displacements in the sublaminates. Then, the energy release rate and its Mode I and Mode II contributions are evaluated. As an example, the model is applied to the analysis of the double cantilever beam test with both symmetric and asymmetric laminated specimens.


Introduction
Delamination, or interlaminar fracture, is the main life-limiting failure mode for fiber-reinforced composite laminates [1]. Delamination cracks may originate from localized defects and propagate due to peak interlaminar stresses. A huge number of analytical, numerical, and experimental studies have been devoted to this problem during the last few decades [2][3][4][5][6]. The phenomenon is also relevant for many similar layered structures, such as glued laminated timber beams [7] and multilayered ceramic composites [8].
Within the framework of linear elastic fracture mechanics (LEFM), a delamination crack is expected to propagate when the energy release rate, G, attains a critical value, or fracture toughness, G c . Since delamination cracks preferentially propagate along the weak interfaces between adjoining plies (laminae), propagation generally involves a mix of the three basic fracture modes (Mode I or opening, Mode II or sliding, and Mode III or tearing). Each fracture mode furnishes an additive contribution to the total energy release rate, G I , G II , and G III [9], and corresponds to a different value of interlaminar fracture toughness, G Ic , G IIc , and G IIIc [10]. Thus, to predict delamination crack propagation, it is necessary (i) to assess by experiments the interlaminar fracture toughness (in pure and mixed fracture modes), (ii) to define a suitable mixed-mode fracture criterion, and (iii) to use a theoretical model to evaluate the energy release rate and its modal contributions [11].
Bending-extension coupling and other elastic couplings are a typical feature of composite laminated beams and plates [49]. Nevertheless, only a few theoretical models for the study of delamination take into account elastic couplings. Among these, it is worth citing the pioneering work by Schapery and Davidson on the prediction of the energy release rate in mixed-mode fracture conditions [50]. In recent years, Xie et al. considered bending-extension coupling in the analysis of delamination toughness specimens [36] and composite laminated plates subjected to flexural loading [25]. Dimitri et al. presented a general formulation of the elastic interface model including bending-extension and shear deformability, but limited their analytical solution to the case with no elastic coupling [26]. Valvo analyzed the delamination of shear-deformable laminated beams with bending-extension coupling based on a rigid interface model [51]. Tsokanas and Loutas extended the above-mentioned analysis to include the effects of crack tip rotations and hygrothermal stresses [52]. To the best of our knowledge, a complete analytical solution for the elastic interface model of delaminated beams with bending-extension coupling and shear deformability has not yet been presented in the literature.
This paper analyses the problem of an interfacial crack in a multi-directional laminated beam with possible bending-extension coupling [53]. A crack-tip element is defined as a laminate segment extending ahead and behind the crack tip cross-section, such that the fracture process zone is fully included within [54,55]. In line with the elastic interface modeling approach, the crack-tip element is considered as an assemblage of two sublaminates connected by an elastic-brittle interface of negligible thickness. Each sublaminate is modeled as an extensible, flexible, and shear-deformable laminated beam [56]. The mathematical problem is reduced to a set of two differential equations in the interfacial stresses. A complete analytical solution is derived including explicit expressions for the internal forces, strain measures, and generalized displacements in the sublaminates. Then, the energy release rate and its Mode I and Mode II contributions are evaluated based on Rice's J -integral [57]. By way of illustration, the model is applied to the analysis of the DCB test. Three laminated specimens are analyzed, made up of uni-directional, cross-ply, and multi-directional sublaminates.

Mechanical Model
Let us consider a composite laminated beam of length L, width B, and thickness H = 2h, with a through-the-width delamination crack of length a. The delamination plane ideally subdivides the laminate into two sublaminates (labeled 1 and 2) of thicknesses H 1 = 2h 1 and H 2 = 2h 2 . We fix a right-handed global Cartesian reference system, Oxyz, with the origin O at the center of the crack tip cross-section and the x-, y-, and z-axes aligned with the laminate longitudinal, width, and thickness directions, respectively ( Figure 1). The present mechanical model focuses on a crack-tip element (CTE), here defined as a laminate segment of length , extending ahead and behind the crack tip cross-section, which fully includes the fracture process zone [54,55]. In practice, the end sections of the CTE, A and B, should be chosen such that A is immediately ahead of the delamination front and B is sufficiently far from it to make sure that the laminate behind behaves as a monolithic beam (Figure 2a). We model the sublaminates as two plane beams connected by an elastic-brittle interface of negligible thickness. The generalized displacements for the sublaminates are defined as follows: u α and w α denote the sublaminate mid-plane displacements along the xand z-directions, respectively; ϕ α denote their cross-section rotations about the y-axis (positive if counter-clockwise). Henceforth, the subscript α is used to refer to the upper (α = 1) and lower (α = 2) sublaminates ( Figure 2b). In line with the kinematics of Timoshenko's beam theory [56], we define the strain measures for the sublaminates as the mid-plane axial strain, α , average shear strain, γ α , and pseudo-curvature, κ α . These are related to the generalized displacements as follows: Next, we introduce the conjugated internal forces, namely the axial force, N α , shear force, Q α , and bending moment, M α . For a laminated beam with bending-extension coupling, the internal forces are related to the strain measures as follows: where A α , B α , C α , and D α respectively denote the sublaminate extensional stiffness, bending-extension coupling stiffness, shear stiffness, and bending stiffness. In general, such stiffnesses should be calculated according to classical laminated plate theory [49], as detailed in Appendix A.
For what follows, it is convenient to introduce also the sublaminate compliances, and write the constitutive relations, Equation (2), in inverse form: The elastic interface between the upper and lower sublaminates transfers both normal and shear interfacial stresses, respectively given by: where k z and k x are the interface elastic constants and: are the transverse and longitudinal relative displacements at the interface, respectively. By substituting Equation (6) into (5), we obtain:

Differential Problem
The equilibrium equations for the upper and lower sublaminates are: where n α , q α , and m α respectively are the distributed axial load, transverse load, and bending couple, due to the normal and shear stresses transferred by the interface: We substitute Equation (9) into (8) and differentiate the latter with respect to x. Then, by using Equations (1)-(4), we obtain the set of governing differential equations for the upper sublaminate, and for the lower sublaminate,

Boundary Conditions
The crack-tip element is subjected to known internal forces on the upper and lower sublaminates at cross-section A. Hence, the following six static boundary conditions can be written: At cross-section B, three kinematic boundary conditions can be obtained from the assumption that the laminate behind behaves as a monolithic beam: To sum up, we have nine boundary conditions to complete the formulation of the differential problem. We note that, in general, no restraints are given to prevent rigid-body motions of the CTE.

Solution Strategy
Following a solution strategy already adopted to solve some specific problems [22,43,44,53], the interfacial stresses are here assumed as the main unknowns. To this aim, we differentiate the normal and shear interfacial stresses, Equation (7), four and three times, respectively, with respect to x.
Then, we substitute Equations (10) and (11) into the result. After some simplifications, we obtain the following differential equation set: where: are constant coefficients and: is a coupling parameter, here called the unbalance parameter by analogy with the literature on adhesive joints [39]. Two cases have to be considered in solving the differential problem: (1) the balanced case, when β 0 = 0, i.e., This condition is trivially fulfilled, for instance, by homogeneous laminates with mid-plane delaminations (for which b 1 = b 2 = 0, d 1 = d 2 , and h 1 = h 2 ) or symmetrically-stacked laminates with mid-plane delaminations (for which b 1 = −b 2 , d 1 = d 2 , and h 1 = h 2 ). Besides, the balance condition is satisfied by homogeneous sublaminates with bending stiffness ratio D 1 D 2 = h 1 h 2 [38,58]. More generally, it is possible to conceive of non-trivial stacking sequences, where the upper and lower sublaminates have unequal thicknesses, but fulfil the balance condition [59,60]. In the balanced case, Equations (14) are uncoupled and can be solved separately to obtain the normal and shear interfacial stresses; (2) the unbalanced case, when β 0 = 0, i.e., This condition corresponds to laminates with generic stacking sequences and arbitrarily-located delaminations (with respect to the mid-plane). In the unbalanced case, Equations (14) are coupled and must be solved simultaneously.
In the following subsections, we will show how to solve the stated differential problem and obtain the expressions for the interfacial stresses. Then, by using Equations (1), (4) and (8), it is straightforward to deduce also the expressions for the internal forces, strain measures, and generalized displacements.

Interfacial Stresses
In the balanced case, the differential problem, Equation (14), reduces to the following two uncoupled differential equations: whose coefficients are defined in Equation (15). The general solution for the interfacial stresses in the balanced case can be written as follows: where f 1 , f 2 , ..., f 7 are integration constants and: are the non-zero roots of the characteristic equation of the differential problem (17): Strictly speaking, the solution expressed by Equation (18) holds provided that the interface constant k z is not equal to the following "critical" value: for which Equation (20) has coincident roots and the solution requires a different expression. However, for material properties corresponding to common fiber-reinforced laminates, the numerical values of k z are usually greater thank z , and the roots of the characteristic equation are real-valued (see Appendix B.1). For very compliant interfaces, e.g., if the present model is applied to the analysis of adhesively-bonded joints with weak adhesive, it may happen that k z <k z , and the characteristic equation has complex roots. In such cases, the present solution is still valid, provided that the complex exponential functions are considered in Equation (18) and the following. Similar considerations apply also to the unbalanced case.

Internal Forces
By substituting Equation (18) into (9), the latter into (8), and integrating with respect to x, we obtain the expressions for the internal forces in the upper sublaminate, and in the lower sublaminate, where f 8 , f 9 , ..., f 13 are integration constants.

Strain Measures
By substituting Equations (22) and (23) into (4), we obtain the expressions for the strain measures in the upper sublaminate, and in the lower sublaminate,

Generalized Displacements
Lastly, by substituting Equations (24) and (25) into (1) and integrating the latter with respect to x, we obtain the expressions for the generalized displacements in the upper sublaminate, and in the lower sublaminate, where f 14 , f 15 , ..., f 19 are further integration constants.

Supernumerary Integration Constants
In the obtained analytical solution, nineteen integration constants appear. However, seven of them are supernumerary as a consequence of the adopted solution strategy (see Section 3.1).
To determine the supernumerary integration constants, we substitute the solution for the interfacial stresses, Equation (18), and generalized displacements, Equations (26) and (27), into Equation (7). In this way, two identities are obtained, which need to be fulfilled for all values of x. The following seven relations among the integration constants are determined: where: The values of the remaining twelve independent integration constants ( f 1 , f 2 , ..., f 9 and f 17 , f 18 , and f 19 ) have to be determined by imposing the boundary conditions of the particular problem being solved. It can be noted that the last three independent integration constants describe a rigid motion of the specimen and can be assessed only when appropriate kinematic restraints are prescribed. The adopted solution strategy enables the determination of the internal forces, strain measures, and relative displacements, even though the values of the last three constants remain undetermined.

Interfacial Stresses
In the unbalanced case, Equations (14) are coupled. To separate the unknowns, we solve the first of such equations w.r.t. to dτ/ dx and substitute the result into the second equation. After some simplifications, we obtain: where: with α 1 , α 2 , α 3 , and β 0 already defined in Equations (15) and (16), respectively. The general solution for the interfacial stresses in the unbalanced case can be written as follows: where g 1 , g 2 , ..., g 7 are integration constants and µ 1 , µ 2 , ..., µ 6 are the non-zero roots of the characteristic equation of the differential problem (30): By substituting Equation (31) into (33), the latter can be rearranged as: where λ 1 , λ 3 , and λ 5 are given by Equation (19). From Equation (34), it can be easily concluded that the roots of the characteristic equation for the unbalanced case are different from those for the balanced case. Conversely, by assuming β 0 = 0, Equation (34) turns out to be equivalent to Equation (20). In general, Equation (33) can be solved numerically. In Appendix B.2, the conditions are given for which µ 1 , µ 2 , ..., µ 6 are real, along with their analytical expressions in this case. It should be noted that, strictly speaking, the solution expressed by Equation (32) is valid only when the characteristic equation has no coincident roots, which is the case for common geometric and material properties.

Internal Forces
By substituting Equation (32) into (9), the latter into (8), and integrating with respect to x, we obtain the expressions for the internal forces in the upper sublaminate, and in the lower sublaminate, where g 8 , g 9 , ..., g 13 are integration constants.

Strain Measures
By substituting Equations (35) and (36) into (4), we obtain the expressions for the strain measures in the upper sublaminate, and in the lower sublaminate,

Generalized Displacements
Lastly, by substituting Equations (37) and (38) into (1) and integrating the latter with respect to x, we obtain the expressions for the generalized displacements in the upper sublaminate, and in the lower sublaminate, where g 14 , g 15 , ..., g 19 are further integration constants.

Supernumerary Integration Constants
As for the balanced case (see Section 3.2.5), nineteen integration constants appear in the analytical solution for the unbalanced case. Again, seven of them are supernumerary and can be determined by substituting the solution for the interfacial stresses, Equation (32), and generalized displacements, Equations (39) and (40), into Equation (7). In this way, two identities are obtained, which need to be fulfilled for all values of x. Consequently, the following seven relations among the integration constants are determined: where α 4 is furnished by Equation (29). The values of the remaining twelve independent integration constants (g 1 , g 2 , ..., g 9 and g 17 , g 18 , and g 19 ) have to be determined by imposing the boundary conditions of the particular problem being solved. Again, the last three independent integration constants describe a rigid motion of the specimen and can be assessed only when appropriate kinematic restraints are prescribed. The adopted solution strategy enables the determination of internal forces, strain measures, and relative displacements, even though the values of the last three constants remain undetermined.

Energy Release Rate
Under I/II mixed-mode fracture conditions, the energy release rate can be written as G = G I + G II , where G I and G II are the contributions related to Fracture Modes I and II, respectively. In LEFM, the energy release rate identifies the path-independent J -integral introduced by Rice [57]. Choosing an integration path encircling the interface between sublaminates, as detailed in Appendix C, we obtain: or, by recalling Equation (5), where σ 0 = σ(0) and τ 0 = τ(0) are the peak values of the interfacial stresses at the delamination front. The two addends in Equation (43) correspond to the Mode I and Mode II contributions to the energy release rate, respectively [41]. Since compressive normal stresses at the crack tip (σ 0 < 0) do not promote crack opening, the modal contributions to the energy release rate can be written as where H(·) denotes the Heaviside step function. By substituting Equation (18) into (44), we obtain the following expressions for the balanced case: (45) Likewise, by substituting Equation (32) into (44), we obtain the following expressions for the unbalanced case: Lastly, to characterize the ratio between the modal contributions to the energy release rate, we introduce the mode-mixity angle [11], ranging from 0 • (pure Mode I) to 90 • (pure Mode II).

Examples
By way of illustration, we applied the developed solution to the analysis of the DCB test ( Figure 3). Laminated specimens with the following geometric sizes were considered: L = 100 mm, B = 25 mm, and H = 3.2 mm. A delamination length a = 25 mm and a reference load P = 100 N were used in all of the following calculations. Boundary conditions (12) and (13) were specialized to the present case by taking the crack tip internal forces as N A 1 = N A 2 = 0, Q A 1 = −Q A 2 = P, and M A 1 = −M A 2 = Pa and assuming a CTE length = L − a. Three laminated specimens with different stacking sequences were ideally built as follows. We started from a typical carbon-epoxy fiber-reinforced lamina of thickness t p = 0.2 mm and elastic moduli shown in Table 1 [61]. Based on this, three stacking sequences were designed corresponding to uni-directional (UD), cross-ply (CP), and multi-directional (MD) sublaminates (Table 2). Lastly, the DCB test specimens were obtained by variously assembling the sublaminates one over the other: The sublaminate plate stiffness matrices, calculated as explained in Appendix A, are reported in Table 3. It can be noted that no bending-extension coupling was exhibited by the UD sublaminate (all terms B ij = 0), while bending-extension coupling was shown by both the CP (terms B xx = −B yy = 0) and MD sublaminates (all terms B ij = 0, except B xs = B ys = 0). The sublaminate equivalent beam stiffnesses, calculated via Equations (A10) and (A19), are given in Table 4. It can be noted that the UD sublaminate had the highest extensional and bending stiffnesses, A and D, while the MD sublaminate had the highest bending-extension coupling stiffness, B. The shear stiffness, C, showed little variations among the chosen stacking sequences. A crucial point for the successful implementation of the present model is the appropriate choice of the values of the elastic interface constants. These can be assigned by conducting a compliance calibration with experimental tests or numerical simulations, if available [62]. In a previous study on asymmetric DCB test specimens [22], compliance calibration with finite element simulations gave: which, in the present case, yielded k x = 21,845 N/mm 3 and k z = 19,016 N/mm 3 . Such values were adopted for the present examples. Figure 4 shows the plots of the interfacial stresses and internal forces as functions of the x-coordinate in the neighborhood of the crack tip for the UD//UD specimen. This specimen was symmetric; hence, the unbalance parameter was β 0 = 0, and the analytical solution for the balanced case applied (see Section 3.2). For the DCB test, the specimen response also corresponded to pure Mode I fracture conditions: actually, Figure 4a shows that the shear stresses were identically null, while the normal stresses attained a peak value at the crack tip cross-section. Figure 4b shows that also the sublaminate axial forces were null. Instead, Figure 4c,d shows the trends of the sublaminate shear forces and bending moments, which started from the assigned values at the crack tip cross-section and then approached zero asymptotically. Shear forces and bending moments had symmetric trends in the upper and lower sublaminates because of the specimen symmetry. Figures 5 and 6 show the plots of the interfacial stresses and internal forces as functions of the x-coordinate in the neighborhood of the crack tip for the UD//CP and UD//MD specimens, respectively. Both specimens were asymmetric with unbalance parameter β 0 = 0. Hence, the analytical solution for the unbalanced case applied (see Section 3.3). For the DCB test, the specimen response corresponded to I/II mixed-mode fracture conditions: actually, both Figures 5a and 6a show non-zero normal and shear stresses, which attained peak values at the crack tip cross-section. All of the internal forces started from the assigned values at the crack tip cross-section and then approached zero asymptotically.   Table 5 reports the Mode I and II contributions to the energy release rate computed for the three analyzed specimens, along with the total energy release rate and mode-mixity angle. It is apparent that the relative amount of the Mode II contribution increased with the specimen asymmetry. Table 5. Energy release rate and mode mixity. (c) (d) Figure 6. UD//MD specimen: (a) interfacial stresses, (b) axial forces, (c) shear forces, and (d) bending moments in the crack tip neighborhood. Figure 7 shows a schematic of the double cantilever beam specimen loaded by uneven bending moments (DCB-UBM) test [63]. Here, this particular test configuration was used to illustrate the effects of interface stiffness, i.e., of the values of the interface elastic constants, k x and k z , on the solution according to the proposed model. The same laminated specimens considered in the examples of Section 5.1 were again analyzed, but subjected to different loading conditions. The upper sublaminate was loaded by a unit bending couple, M 1 = 1 N m, while the lower sublaminate was loaded by a variable bending couple, M 2 . Figure 8 shows the dimensionless Mode I and Mode II contributions to the energy release rate, γ I = G I /G and γ II = G II /G, as functions of the bending moment ratio, M 2 /M 1 .

Effects of Interface Stiffness
For the UD//UD specimen (Figure 8a), the modal contributions to G turned out to be independent of the interface stiffness. This result was due to the mid-plane symmetry of the specimen and the independence of the crack tip bending moments on the delamination length. In this case, also the mode mixity depended uniquely on the bending moment ratio, M 2 /M 1 . In particular, it can be noted how pure Mode I fracture conditions (γ II = 0) were obtained for M 2 = −M 1 , while pure Mode II fracture conditions (γ I = 0) were obtained for M 2 = M 1 . Values of M 2 /M 1 greater than one also corresponded to pure Mode II because of the incompressibility constraint introduced for G I in Equation (44). Instead, for both the UD//CP (Figure 8b) and UD//MD specimens (Figure 8c), the modal contributions to G depended on the values of the interface elastic constants, k x and k z . For each modal contribution, the plots show five curves with different shades of blue: from lighter to darker, and these correspond to k x = k z = 10 2 , 10 3 , 10 4 , and 10 5 N/mm 3 ; the darkest curves represent the limit case k x = k z − → +∞, which corresponds to a rigid interface model [51]. For the specimens considered, also the mode mixity and the conditions for pure fracture modes depended on the interface stiffness.
The range of values adopted for k x and k z in the previous plots may look quite wide. Indeed, this range corresponds to real values that may be obtained when performing a compliance calibration with experimental results or numerical simulations [62]. The results obtained suggest that the elastic interface model is potentially capable of highly-accurate prediction and interpretation of experimental results, provided that the appropriate values of the interface elastic constants are selected. To this aim, it would be useful to validate the theoretical predictions obtained with experimental results, such as those obtained by Davidson et al. [64,65]. Such a comparison, similar to the one carried out by Harvey et al. to validate different mixed-mode partition theories [66,67], will be the subject of future work.

Conclusions
An elastic interface model was introduced for the analysis of delamination in bending-extension coupled laminated beams with through-the-width delamination cracks. The analytical solution obtained applies to a wide class of multi-directional fiber-reinforced composite laminates, including cross-ply laminates, provided that their stacking sequences allow them to be modeled as plane beams. Examples illustrated the application of the model to the analysis of the double cantilever beam test in I/II mixed-mode fracture conditions. Furthermore, the effects of the interface stiffness on the model predictions were discussed with reference to the DCB specimen loaded with uneven bending moments.
By suitably adapting the boundary conditions, the outlined solution can be easily used to analyze different delamination toughness test configurations, as well as delamination behavior in real structural elements and debonding of adhesively-bonded joints.
The proposed model is potentially capable of highly-accurate prediction and interpretation of experimental results, provided that the appropriate interface stiffness is defined. In this respect, it will be useful to compare the model predictions with the results of experimental tests available in the literature. Future developments also include the extension of the model to fully-coupled laminates, for which plate theory must be used, and I/II/III mixed-mode fracture conditions are expected.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Laminated Beam Stiffnesses
For a homogenous and special orthotropic beam of width B and thickness H, the extensional, bending-coupling, shear, and bending stiffnesses can be calculated respectively as: where E x and G zx respectively are the longitudinal Young's modulus and the transverse shear modulus of the material. For laminated beams with generic stacking sequences, the above quantities have to be calculated as explained in the following sections.

Appendix A.1. Extension and Bending Stiffnesses
For a laminated plate with a generic stacking sequence, the constitutive laws are: where N x , N y , and N xy are the in-plane forces and M x , M y , and M xy are the bending moments (per unit plate width); correspondingly, x , y , and γ xy are the in-plane strains and κ x , κ y , and κ xy are the bending curvatures; A ij , B ij , and D ij , with i, j ∈ {x, y, s}, are the coefficients of the extensional, bending-extension coupling, and bending stiffness matrices, respectively, to be computed following classical laminated plate theory [49]: where Q (k) ij , with i, j ∈ {x, y, s}, are the in-plane elastic moduli of the kth lamina, included between the ordinates z = z k−1 and z = z k , in the global reference system, and n is the total number of laminae. For an orthotropic lamina, whose principal material axes, 1 and 2, are rotated by an angle θ k with respect to the global axes x and y, where the in-plane elastic moduli of the lamina in its principal material reference system are: , and Q (k) 21 , and G (k) 12 being the in-plane elastic moduli of the lamina in engineering notation. For a plate to be modeled as a plane beam in the zx-plane, the out-of-plane internal forces and moments must be N y = N xy = 0 and M y = M xy = 0. Besides, the shear-extension and bend-twist coupling terms must vanish: A xs = A ys = 0 and D xs = D ys = 0. The bending-shear coupling terms must also be null: B xs = B ys = 0. Such conditions are fulfilled by laminates with single or multiple special orthotropic layers and antisymmetric cross-ply laminates. Furthermore, antisymmetric angle-ply and generic multi-directional laminates may satisfy the above-mentioned conditions for particular stacking sequences. As a consequence, Equations (A2) and (A3) reduce to: and: For a laminated beam of width B, the axial force and bending moment respectively are N = BN x and M = BM x ; the axial strain and pseudo-curvature are = x and κ = κ x . Thus, from Equations (A7) and (A8), the constitutive laws are derived as follows: where the laminated beam stiffnesses are evaluated as:

. Shear Stiffness
The shear stiffness of a laminated plate can be evaluated by imposing a strain energy equivalence, as usually done for homogeneous beams [56] and recently extended to laminated plates by Xie et al. [25]. Here, we detail such a procedure for a laminated beam [68][69][70]. To this aim, first we recall the expression of the normal stress in the x-direction within the kth lamina [49]: xs γ xy + zκ xy . (A12) Next, we substitute the expressions of the plate strain measures obtained in the previous section into Equation (A12). Thus, we obtain: where: (A14) Following the classical derivation of Jourawski's shear stress formula [56], we impose the static equilibrium in the x-direction of a short portion of beam defined by a given value of the z-coordinate. Thus, by recalling the inverse constitutive Equation (4) and static equilibrium Equation (8), with no distributed loads, the following expression is obtained for the shear stress in the kth lamina: where Q is the transverse shear force (acting in the z-direction) and: R k = S k z k + T k z 2 k + n ∑ j=k+1 S j z j − z j−1 + T j z 2 j − z 2 κ0 , and T k = where b and d respectively are the bending-extension and bending compliances of the laminated beam, as given by Equation (3). The shear compliance, c, can now be calculated by imposing the following energy equivalence for a short beam segment: where: G where Γ is an arbitrary integration path surrounding the crack tip, ω is the strain energy density (per unit volume), t x and t z are the stress vector components (referred to the outer normal to Γ), and u and w are the displacement vector components along the coordinates x and z.
For the CTE, we choose an integration path, Γ, encircling the interface between sublaminates as depicted in Figure A1. The path is subdivided into three segments, Γ 12 , Γ 23 , and Γ 34 , along which quantities entering the J -integral have the expressions listed in Table A1.
Accordingly, the following expression results: where t is the interface thickness and Φ = Φ(∆w, ∆u) is the interface potential energy [27], such that: By recalling the definition of relative displacements, Equation (6), and assuming a vanishing interface thickness (t − → 0), we obtain: Next, by substituting Equation (A28) into (A29) and recalling that dΦ is an exact differential, we deduce: which holds for a general cohesive interface with potential-based traction-separation laws. For a brittle-elastic interface, the potential energy function specializes to: Φ(∆w, ∆u) = 1 2 k z ∆w 2 + 1 2 k x ∆u 2 . (A32) By substituting Equation (A32) into (A31), we obtain: which is identical to Equation (42).
For the present elastic interface model, the two addends in Equation (A33) naturally identify the Mode I and II contributions to the energy release rate (see Section 4.1). This result holds for both the balanced and unbalanced cases, for which the analytical solution has been derived in Section 3. For a general cohesive interface model, Equation (A31) yields the value of the J -integral, but its decomposition into Mode I and II contributions is not straightforward. The interested reader can find further hints about this issue in the papers by Wu et al. [37,38] and in the references recalled therein.