Influence of boundary conditions on the response of multilayered plates with cohesive interfaces and delaminations using a homogenized approach

Stress and displacement fields in multilayered composites with interfacial imperfections, such as imperfect bonding of the layers or delaminations, or where the plies are separated by thin interlayers allowing relative motion, have large variations in the thickness, with characteristic zigzag patterns and jumps at the layer interfaces. These effects are well captured by a model recently formulated by the author for multilayered plates with imperfect interfaces and affine interfacial traction laws (Massabò & Campi, Meccanica, 2014, in press; Compos Struct, 2014, 116, 311-324). The model defines a homogenized displacement field, which satisfies interfacial continuity, and uses a variational technique to derive equilibrium equations depending on only six generalized displacement functions, for any arbitrary numbers of layers and interfaces. The model accurately predicts stresses and displacements in simply supported, highly anisotropic, thick plates with continuous, sliding interfaces. In this paper the model is applied to wide plates with clamped edges and some inconsistencies, which have been noted in the literature for models based on similar approaches and have limited their utilization, are explained. A generalized transverse shear force is introduced as the gross stress resultant which is directly related to the bending moment in the equilibrium equations of multilayered structures with imperfect interfaces and substitutes for the shear force of single-layer theory. An application to a delaminated wide plate highlights the potential and limitations of the proposed model for the solution of fracture mechanics problems.


INTRODUCTION
tress and displacement fields in multilayered composites with interfacial imperfections, such as imperfect bonding of the layers or delaminations, or where the plies are separated by thin interlayers allowing relative motion, have large variations in the thickness, with characteristic zigzag patterns and jumps at the interfaces.These effects cannot be captured using classical first-or higher-order single-layer theories and require models based on a discrete-layer approach, where the number of unknowns is typically large and depends on the number of layers/interfaces and on the layer kinematic fields.This limits the range of problems which can be solved analytically and makes computational solutions necessary for most cases, in particular when the status of the interfaces evolves during loading due to delamination fracture [1][2][3][4][5][6][7].Zigzag theories [8][9][10] were originally proposed for multilayered systems with perfectly bonded interfaces in order to overcome the limitations of discrete-layer approaches and satisfy continuity of transverse and normal stresses at the S interfaces.The theories, through the imposition of interfacial continuity conditions, define a homogenized displacement field which depends on a limited number of unknowns and is able to reproduce through-thickness zigzag patterns due to the material inhomogeneities.Later, the zigzag theories were extended to describe plates and shells with imperfectly bonded, purely elastic interfaces in [11][12][13][14].The extended theories however manifest a number of inconsistencies: (i) they are unable to reproduce the expected transitions in the internal gross stress resultants on varying the stiffness of the imperfect interfaces (first noted in [11,15,16]); (ii) they give rise to unrealistic effects in the transverse displacements, which are larger than those of fully debonded plates in partially bonded plates (this effect was defined shear-locking in [15]); (iii) they show some inconsistencies in plates with clamped edges (noted in [17]).Recently the author formulated in [18,19] a theory which, starting from those in [11][12][13][14], extends the formulation to plates and beams with interfaces characterized by affine interfacial traction laws (to describe piecewise linear cohesive functions) and accounts for the energy contribution of the imperfect interfaces in the derivation of the equilibrium equations.This contribution was erroneously omitted in the original theories and in all theories derived later from the original models (see [18] for a list).Corrected formulations of the models [11][12][13][14] are presented in the Appendices in [18].The accuracy of the theory proposed in [18,19] has been verified against exact 2D elasticity solutions in highly anisotropic, simply supported, multilayered plates, with continuous sliding interfaces and deforming in cylindrical bending.The model accurately describes stress and displacement fields in plates with different numbers of interfaces, equally and unequally spaced in the thickness, over the whole range of interfacial stiffnesses, from fully bonded to fully debonded.A limitation of the theory has been observed when dealing with very thick, highly anisotropic plates with compliant interfaces, where the shear deformations are underestimated in the derivation of the transverse displacements as a consequence of the assumed continuity between interfacial tractions and shear stresses.It is expected that this problem could be solved using a shear correction factor which depends on the interfacial properties (work in progress).In this paper, the issues noted in [17] in beams with clamped ends, where the zigzag theory proposed in [9,10] for fully bonded systems shows some apparent inconsistencies in the transverse shear force, are discussed; and the equilibrium equations derived in [19] for plates in cylindrical bending are restated to clarify the problem.The new formulation introduces a generalized transverse shear force, which is the gross stress resultant directly related to the bending moment in the equilibrium equations of multilayered plates with imperfect interfaces and substitutes for the shear force of singlelayer theory.Finally, a delaminated cantilevered wide plate with a clamped edge is studied as a preliminary investigation of the applicability of the model to fracture mechanics problems.

MODEL
onsider a rectangular multilayered plate of thickness h and in-plane dimensions 1 L and 2  L L , with x , is introduced with the axis 3 x normal to the reference surface of the plate, which is arbitrarily chosen, and measured from it (Fig. 1).The plate consists of n layers exhibiting different mechanical properties and joined by 1  n interfaces, which are described as mathematical surfaces where the material properties and the displacements may change discontinuously while the interfacial tractions are continuous.The layer k, where the index 1  ,.., k n is numbered from bottom to top, is defined by the coordinates , and has thickness ( ) k h , Fig. 1 (the k superscript in brackets identifies affiliation with layer k).Each layer is linearly elastic, homogeneous and orthotropic with material axes parallel to the geometrical axes.The displacement vector of an arbitrary point of the plate at the coordinate The plate is subjected to distributed loads acting on the upper and lower surfaces, x and deforms in cylindrical bending.In addition, the plate is assumed to be incompressible in the thickness direction and the interfaces to be rigid against mode I (opening) relative displacements.This latter assumption, which is often used in the literature, is rigorously correct only in problems where the conditions along the interfaces are purely mode II.The assumption, however, is acceptable in the presence of continuous interfaces, when the interfacial normal tractions are small compared to the tangential tractions and interfacial opening is prevented, e.g. by a through-thickness reinforcement or other means.To describe cohesive delaminations under general mixed mode conditions, the general treatment, which has been proposed in [18], for plates, and in [19], for wide plates in cylindrical bending, and accounts for interfacial opening, must be applied.Based on the C assumptions above, the displacement components then simplify in Following the classical assumption of lower-order plate theories, in the constitutive relationships for the generic layer k the normal stresses, 33  ( ) k for k =1..n-1, are assumed to be negligibly small compared to the other components.This yields: , where the ( ) A are the coefficients of the 6×6 stiffness and compliance matrices (engineering notation).Transverse normal tractions will then be derived a posteriori from local equilibrium.The interfaces are described by interfacial traction laws which relate the interfacial shear tractions, acting along the surface of the layer k at the interface with unit positive normal vector, to the interface relative sliding displacement: The interfacial traction law is generally nonlinear to represent different physical mechanisms, which may include the elastic response of thin interfacial layers, cohesive/bridging mechanisms developed by trans-laminar reinforcements or other means, material rupture, elastic contact along the delamination surfaces [3][4][5][6][7]20,21].The law can be approximated as a piecewise linear function so that the arbitrary branch i is described by an affine function of the relative displacement: where i k S K and i k S B are the interface tangential stiffness and compliance and i k S t is a constant interfacial traction which is assumed to act when 2 0 ˆk v  , and is typically positive/negative for positive/negative 2 ˆk v , i.e.
, the affine law of Eq. ( 4) could describe the bridging mechanisms developed by a through-thickness reinforcement, e.g.stitching, applied to a laminated composite [21].If the relationship ( 4) is used to represent branches of cohesive traction laws, different linear functions may be needed to represent processes inducing loading and unloading of the delaminations and the associated sliding and reverse-sliding mechanisms.In this paper, equilibrium equations will be derived for plates with interfaces described by the arbitrary branch i of Eq. ( 4) and, for the sake of simplicity, the superscript i will be removed.

Two length-scales displacement field
The displacement field is assumed to be given by the superposition of a global field and local perturbation terms (or enrichments).The nonzero components of the displacement vector at an arbitrary point   where C , and, when the reference surface coincides with the mid-surface of the bottom layer, define the generalized displacement components of its points.The third term in Eq. (5a), with summations on the n-1 of interfaces, supply the zig-zag contributions, 2  k , [9,10], which are continuous in 3 x but with jumps in the first derivatives at the interfaces, C , and are necessary to satisfy continuity of the shear tractions at the interfaces in plates with arbitrary stacking sequences; the fourth term, with summations on the n-1 interfaces, define a discontinuous field and supply the contribution of the relative displacements (jumps) at the cohesive interfaces.Eq. ( 5) define a first order model, since the displacements are piecewise linear functions of 3 x .Higher order models have been proposed in the theories in [11,13], which however have no advantages over I order models in the presence of imperfect interfaces, as it was proven in [18].The linear infinitesimal nonzero strain components at the coordinate 3 x within layer k are derived using Eq. ( 5): where the comma followed by a subscript denotes a derivative with respect to the corresponding coordinate.

Homogenized displacement field
The n-1 unknown zigzag functions 2 2  ( ) k x for k = 1..n-1 in Eq. (5a) are determined as functions of the global displacement variables and displacement jumps by imposing continuity of the shear tractions, Eq. ( 2), across the layer interfaces, which yields: Through Eq. ( 7), ( 6) and ( 1) the kth zig-zag function is defined in terms of the global displacement variables: where: Once the functions 2 2  ( ) k x , for k =1...n-1, have been defined, the relative displacement at each cohesive interface, 2 ˆk v , is defined through the constitutive law of the interface, Eq. (4b), using Eq.(1a), ( 2), (6b) and (8).The expression for the displacement jump at the  S ( ) k interface at the coordinate 3 k x in terms of the global displacement variables is: Eq. ( 8) and ( 10) can then be inserted into Eq.( 5) to obtain the homogenized displacement field.The displacement components within layer k are: where: ( ) Eq. ( 12) highlight that the displacement field is fully defined by the 3 displacement variables, 02 0 2  , , v w , which describe the global part of the field, and are underlined in the equations, and by parameters which depend on the elastic constants of the material, the layup and the geometry (no line) and parameters depending on the properties of the interfaces through the assumed interfacial traction laws (curved line on top).For perfectly bonded layers, when  k S B 0 for k = 1..n-1, all terms with the curved line on top vanish and the equations are those of first order zig-zag theory [9,10].The strain components in the layer k in terms of the homogenized displacement variables are: The interfacial tractions at the x , in terms of the homogenized displacements are:

Equilibrium equations
Equilibrium equations and boundary conditions are derived in weak form through the Principle of Virtual Works: where V is the volume of the plate and i = 2,3; the virtual displacements are assumed to be independent and arbitrary and to satisfy compatibility conditions.The first term on the left hand side defines the strain energy in the volume of the body; the second term, with the flat line on top, the energy contributions due to the interfacial tractions on the n-1 interfaces.The last terms define the work done by the external forces, with S have been assumed to be zero (refer to [18,19] for more general loading conditions).
The term related to the interfacial tractions was not present in the models where this approach was first proposed for plates with linear-elastic interfaces [11][12][13][14].The terms were also missing in all subsequent models which extended the theories to different problems (see list in [18]).It has been recently proved in [18] that, in the absence of these terms, the solutions are accurate only in the limiting cases of fully bonded and fully debonded interfaces.
Virtual strains and displacements in Eq. ( 16) are defined as functions of the global displacement variables using Eq. ( 10), ( 12) and ( 14).Then, by applying Green's theorem wherever possible and after lengthy calculations, Eq. ( 16) yields the equilibrium equations and boundary conditions.Dynamic equilibrium equations and boundary conditions for multilayered plates with arbitrary stacking sequences, mixed mode interfaces and under arbitrary loadings have been derived in [18]; a particularization of the equations to plates deforming in cylindrical bending have been presented in [19].
The equilibrium equations for wide plates with sliding only interfaces and quasi-static loading with 2 0    , S S F are presented here in a form that highlights similarities and differences with respect to single-layer theory: T n , to the vertical equilibrant of the external forces acting on the portion of the plate to the right of the sections (Fig. 1a):  normal force:  bending moment:  generalized shear force: where  transverse shear force:  gross resultants and couples associated to the multilayered structure:  gross resultants and couples associated to the cohesive interfaces: The terms in the Eq.(18e) vanish in unidirectionally reinforced systems and those in Eq. (18f) vanish in fully bonded systems.The generalized shear force then coincides with the resultant of the transverse shear stresses, , when the layers have the same elastic constants, namely when    , : where Terms with the tilde define prescribed values of generalized displacements and gross forces and couples applied to B .
Equilibrium and boundary conditions can be expressed in terms of the homogenized displacement variables using the constitutive and compatibility Eq. ( 1), ( 12), ( 14) and ( 15).The equations are presented in [19].
Eq. ( 14a) and (9) show that the transverse shear stress, 23  , obtained from the shear strains, 23  , through the constitutive Eq. ( 1), is constant in the thickness and related to the transverse shear force, Eq. (18d), through 23 2   / b Q h .This stress does not describe the effective status of the material, but for the limit case of a system with perfectly bonded interfaces and layers with the same elastic constants, where In the presence of imperfect interfaces, 23  follows the dependence of the interfacial tractions on the stiffness of the interfaces, due to the imposed continuity, Eq. ( 7)- (10), and progressively goes to zero when the stiffness of the interfaces decreases; in fully bonded systems with a multilayered structure, where 2 .The generalized transverse shear stress, 23 g  , averages the actual shear stress distribution which can be obtained a posteriori from the bending stresses by satisfying local equilibrium, , so that .In order to account for the correct shear deformations in the solution of the differential equations, a shear correction factor, 2 K , can be introduced such that K is equal to 5/6 in fully bonded unidirectionally reinforced plates, to account for the approximated constant distribution of 23  in the thickness, and it becomes a problem dependent parameter in multilayered plates, e.g.[22]; in [9] it was shown that, for simply supported plates with common layups and geometrical/loading conditions, the homogenized zigzag theory with K =5/6 (see also [18][19]).

APPLICATIONS
Highly anisotropic, simply supported, multilayered wide plates with imperfect interfaces he homogenized model for multilayered plates with imperfect interfaces has been verified against exact 2D elasticity solutions in [18,19].Simply supported, highly anisotropic multilayered plates, loaded quasi-statically, with one or more weak layers have been examined on varying the properties of both layers and interfaces.The whole transition between fully bonded and fully debonded plates has been considered.The diagrams in Fig. 2 show exemplary results taken from [19] and refer to a highly anisotropic, thick plate, with two imperfect interfaces having different stiffnesses.

Wide plates with clamped edges
The zigzag theory for fully bonded plates [9,10] was applied in [17] to a multilayered cantilever beam subjected to a concentrated force F at the free end.The authors noted that the shear force was not constant along the beam length, as they would have expected given the linear distribution of the bending moment.In addition, the shear force increased from zero (at the clamped edge) to an asymptotic value higher than F.These apparent inconsistencies are explained by Eq. (17b), (18c) and (19c), which show that the internal gross stress resultant related to the bending moment, through Eq. (17b), is the generalized transverse shear force, Eq. (18c), which is statically equivalent to the vertical equilibrant of the external forces acting on the portion of the plate to the right of each arbitrary cross section, 2 g Q F  .
The diagrams in Fig. 3  Q , correctly coincides with the external applied force F at all coordinates and since the plate is multilayered and fully bonded).Diagram (c) compares the values of the interfacial tractions along the plate length, obtained a posteriori from local equilibrium, with those obtained with a classical discrete layer approach [6][7][8] and shows the existence of a boundary region near the clamped edge where the interfacial tractions predicted through the homogenized approach are not correctly described (bending stresses, not shown, are accurately T predicted in all cases).This is a consequence of the geometric boundary conditions imposed at the clamped boundary.The size of this boundary region depends on the interfacial stiffness and is negligible for very stiff and very compliant interfaces.In a multilayered plate the size of the boundary region is nonzero even when the layers are fully bonded, since

Plates with delaminations
As a preliminary investigation of the applicability of the theory to fracture problems, the interface in the cantilevered plate studied before has been assumed to be fully bonded, for x  ; this is due to the imposition of the continuity conditions on the homogenized displacement variables only (see Eq. (12a) and ( 19)).The incompatibility produces unreliable predictions of the stresses in a very small region localized at the crack tip, of size 50 / L  .Accurate predictions of energy release rate and stress intensity factors are obtained using expressions derived for orthotropic layers in [23], which depend on stress resultants and couples at the crack tip.

CONCLUSIONS
quilibrium equations were derived in [18,19] for multilayered composite plates with cohesive interfaces and delaminations which depend on only six unknown displacement functions (three for wide plates) for any arbitrary numbers of layers and interfaces.The equations have been particularized here to plates deforming in cylindrical bending and restated in a form similar to that of single-layer theory.This introduces a generalized transverse shear force which is directly related to the bending moment, as the shear force is in single-layer theory, and depends on the multilayered structure of the material and the status of the interfaces.The new equations explain the apparent inconsistencies which have been observed in the shear force when using zigzag theories to model plates with clamped edges.Applications to cantilevered wide plates with imperfect interfaces and delaminations confirm the accuracy of the proposed model in predicting stress and displacement fields in thick, highly anisotropic, multilayered plates.They also highlight the existence of boundary regions, near the clamped edges and at the delamination tips, where gross stress resultants and couples are accurately predicted, while stresses and displacements in the layers are not, as a consequence of the imposition of boundary/continuity conditions on the global displacement variables only.The size of the boundary region at the delamination tip in a unidirectionally reinforced laminate is very small, 50 / L  with L the characteristic inplane dimension.The presence of the boundary regions must be accounted for in the solution of the problems and fracture mechanics predictions must rely on expressions depending on gross stress resultants and couples, which should be calculated at the boundary of the region surrounding the crack tip.It is expected that improvements in the prediction of stresses and displacements in the boundary regions may be obtained through the introduction of a shear factor which depends on the status of the interfaces.

S
, and on the lateral bounding surface, B , is in plane strain conditions parallel to the plane 2 3  x


and the terms on the right hand side of Eq. (5a) denote different contributions in the displacement representation: 02 02 2 define standard first order shear deformation theory terms, which are continuous with continuous derivatives in the thickness direction, 1 3

Figure 1 :
Figure 1: (a) Composite plate showing discretization into layers, imperfect/cohesive interfaces and delaminations.(b) infinitesimal element of layer k showing stress resultants and couples and interfacial tractions.

3 SF  (top), 3 SFF
 (bottom) and B i (lateral) the components of the forces acting along the bounding surfaces of the plate,

are normal force and bending moment in the 2 x
direction and 2 g Q is a generalized transverse shear force which is statically equivalent, at any arbitrary sections of the plate with outward normal n =  

2 z Q and 22 zMand 1 2 
in this case the equilibrium equations coincide with those of single layer theory.In all other cases, the generalized shear force depends on the multilayered structure, through , and on the status of the cohesive interfaces, through 22 S M ˆ, and the classical relation between bending moment and shear force of single-layer theory is modified as in Eq. (17b).The mechanical/geometrical boundary conditions on C , at 2 0  , x L, with n =   normal, are: 02  / b Q h again does not describe the actual stress distribution but for the special case of layers with the same elastic constants where observations above, a generalized transverse shear stress can be introduced, which is the relevant internal stress for strength predictions, Similarly, the transverse shear strain, which is related to the transverse shear stress through the constitutive Eq. (1), describes the shear deformations of the plate whose correct measure within this model is given by a generalized shear strain

2 K 2 K 2 K
= 1 leads to accurate predictions of the displacement field.In plates with imperfect interfaces, must depend on the stiffness of the interfaces.Work in in progress on the derivation of and results are presented here for2

Figure 2 .
Figure 2. Simply supported wide plate with L/h = 4 subjected to a sinusoidal transverse load,

. 2 bQ
refer to a cantilevered wide plate of length L/h = 10 made with two unidirectionally reinforced layers with elastic constants, Diagram (a) depicts the transverse shear force, Eq. (18d), along the plate length for different values of the elastic interfacial stiffness and highlights the apparent inconsistency noted in[17] for a fully bonded multilayered plate.Note that 2 b Q is forced to be zero at the boundary by the geometric boundary conditions at 2 0x  , (b) highlights that the generalized transverse shear force in Eq. (18c),2 g . (18c).The diagrams (d) and (e) show bending and transverse shear stresses at the mid-span for different values of the interfacial stiffness.Predictions are accurate in all cases.

Figure 3 .;
Figure 3. Wide plate clamped at 2 0  x and subjected to a concentrated force F at 2 x L  .Two layers, unidirectionally reinforced

)
Homogenized equilibrium equations have been derived for the two regions in terms of the homogenized displacement variables, , and continuity conditions applied at the delamination tip, at 2x L  .The model accurately predicts gross stress resultants/couples and stress components.Bending and transverse shear stresses are depicted by the solid curves (thick lines for 2 in Figs.3d,e.Incompatible displacements are predicted in the layers to the immediate right and left of the cross section at 2