Analysis of bonded joints with laminated adherends by a variable kinematics layerwise model

International Journal of Adhesion & Adhesives 50 (2014) 244–254

Contents lists available at ScienceDirect

International Journal of Adhesion & Adhesives journal homepage: www.elsevier.com/locate/ijadhadh

Analysis of bonded joints with laminated adherends by a variable kinematics layerwise model U. Icardi a,n, F. Sola b a b

Dipartimento di Ingegneria Meccanica e Aerospaziale, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy Dipartimento di Ingegneria Meccanica e Aerospaziale, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy

art ic l e i nf o

a b s t r a c t

Article history: Accepted 6 February 2014 Available online 14 February 2014

A displacement-based zig-zag plate model with variable in-plane and through-the-thickness representation and fixed degrees of freedom is developed and applied to the analysis of bonded joints. Just the in-plane displacements and the shear rotations of the middle plane are used as primary variables. The continuity functions enable an a priori fulfilment of the out-of-plane stress contact conditions at the interfaces between adjacent layers. High-order contributions to displacements are included to meet the stress boundary conditions at the upper and lower faces and to allow the in-plane representation to be varied, e.g. from the adherends to the overlap. In this way, it is possible to better simulate the variation of solutions and the stress boundary conditions at the ends of the overlap. Closed form expressions of these quantities are obtained using symbolic calculus, which, once incorporated in the model, simplify the obtainment of governing equations and speed-up computations. As the representation can vary from point to point, the present model permits an accurate analysis of laminates with general boundary conditions and of bonded joints under a unified approach. Applications are presented to sample cases of bonded joints with laminated adherends using appropriate series expansions of displacements. Linear and nonlinear benchmark test cases for single- and double-lap joints taken from the literature are considered. The results show that accurate stress predictions are computed with a low computational effort in all the cases considered. A good accuracy is achieved just using one component in the series expansion, which implies solving a 3
3 system. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Zig-zag model Hierarchic representation Fixed d.o.f. Bonded joints Deformable adherends

1. Introduction Bonded joints offer a combination of improved performance and manufacturing advantages with respect to traditional mechanical fastening, like better vibration isolation, accommodation of thermal expansion mismatch and hygrothermal swelling, better aesthetic appearance and better aerodynamic capabilities. In particular, adhesive bonding gives a more gradual transfer of load between the structural elements, thus achieving a more uniform stress distribution. As a direct consequence, adhesive bonding is the most suited technique for joining components when stress concentrations due to mechanical fastening have potentially catastrophic effects on strength and fatigue life, like for composite laminates as discussed by Her [1]. Simulating bonded joints is a rather complex matter that requires in many cases to account for geometric and material nonlinearity, to take into consideration the stress boundary

n

Corresponding author. Tel.: þ 39 110906872; fax: þ39 110906899. E-mail addresses: [email protected] (U. Icardi), [email protected] (F. Sola).

http://dx.doi.org/10.1016/j.ijadhadh.2014.02.003 0143-7496 & 2014 Elsevier Ltd. All rights reserved.

conditions in a point form and to carry out the analysis with a sophisticate modelling of out-of-plane stress and strain fields. As shown by the comprehensive literature reviews by Kuno [2], Vinson [3] and by da Silva et al. [4], many analytical, numerical and experimental studies have been published over years. A serious design concern of adhesively bonded joints is the evaluation of the stress and strain fields across the joint and how they are influenced by geometry, materials, loading conditions, temperature and moisture effects. A considerable effort has been put forward by many researchers to fully understand the intricacies in the behaviour of bonded joints and the degree to which the stress states and the failure modes are influenced. The analysis of bonded joints can be carried out via finite element analysis (FEA) without introducing any simplifying assumption that could limit accuracy even in the case of complex joint geometries and complex material models including nonlinear effects. However FEA requires a quite long preparation and a large computing time. To overcome this problem, analytical models (AM) based on simplifying assumptions and finite-difference schemes are extensively used for studying joints. Analytical [5–10], finite element [11–14] and finite-difference [15] solutions show that even when

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

the adherends and the adhesive materials are assimilable as homogeneous and isotropic, stresses and deformations are complex due to the different elastic moduli of adhesive and adherends and the enforcement of the stress boundary conditions. Results of the analyses show non-uniform shear strain and stress distributions in the adhesive layer caused by the progressive reduction of the strain in the adherends along the overlap and the continuity of the adhesive/adherends interface. The out-of plane transverse shear and normal stresses in single- (SLJ) and double-lap (DLJ) bonded joints have a peak close to the edges of the bonding layer, which can lead to the premature failure of the joint in service. The peeling stress can be larger than the shearing stress, thus becoming the dominant effect at the edges, as shown, e.g., by Nemes and Lachaud [16]. The readers can find a comprehensive discussion and extensive assessments in the papers by da Silva et al. [4,17]. Most AM are stress-based models that still maintain some of the simplifying assumptions of the pioneering models (Volkersen [18], Goland and Reissner [19] and Hart-Smith [20]) in order to allow a fast, low cost solution in closed form. Often shear stress is assumed constant across the thickness of the adhesive, or shear and peel stresses in the adhesive layer are calculated solving a plane strain problem, because the adhesive is assumed to deform only in shear and the adherends are assumed to be rigid, though these hypotheses could be not realistic in all cases. With the advent of bonded laminated composites as primary structures in the 80’s, more complex AM that consider the effects due to the deformation of the adherends and of the adhesive have been developed in order to adequately treat joints with laminated adherends. In this case, it is of primary importance to accurately account for the large variation of the transverse shear and normal stresses across the thickness that rise in consequence of the different properties of the layers constituting the joint. However, trying to keep these improved models as simple as possible, often the stress-free boundary conditions at the ends of the overlap are still disregarded, along with the bending effect due to the eccentric load path of SLJ. Moreover, the deformability of adherends is still accounted for with simplified techniques, as in [19,20]. Refined AM with improved predictive capability are reviewed by da Silva et al. [4], Gustafson et al. [5], Diaz Diaz et al. [6] and Radice and Vinson [7]. These models, depending on the complexity of their governing differential equations, can be solved either in closed or numerical form. For examples, the models by Renton and Vinson [21], Srinivas [22] and Allman [23] are cited as example that account for the transverse shear and normal deformations of adhesive and adherends, satisfy the stress-free boundary conditions ([21,22]) and consider the effects of bending, stretching and shearing in the adherends and the tearing actions in the adhesive ([23]). Bonded joints with adherends made of composite materials have been recently studied by Diaz Diaz et al. [6] using a stress based layerwise model obtained by stacking Reissner–Mindlin plates, which requires to solve a system of 29 equations in 29 unknown parameters. Because the assumptions of AM can affect the accuracy of results, an increasing number of 3-D FEA and of finite-difference 3-D schemes results for bonded joints have been published in the last years. The papers by Khalili et al. [12] and Diaz et al. [13] dealing with 3-D FEA of SLJ with CFRP adherends and epoxy adhesive and the paper by Andruet et al. [11] dealing with geometric nonlinear effects are cited as examples. Finite-difference 3-D schemes as an alternative to FEA and AM are overviewed by Xu and Li [15]. As discussed by Adams and Mallick [24], whenever a complex representation aimed at realistically simulating the stress fields is employed, a numerical scheme should be used for solving, even in the case of AM. However, this does not result disadvantageous and unpractical compared to simpler AM allowing for closed-form solutions, because costs are still saved with respect to 3-D FEA without significant accuracy loss. Besides, AM are not affected by stress singularities at the edge interfaces like 3-D FEA.

245

Geometric and material nonlinear effects and adherends with dissimilar thickness were recently considered by Mortensen and Thomsen [25] and Smeltzer [26], using refined AM. Damage of joints under impact loading was considered by Vaidya et al. [27]. The classical laminated plate theory (CLPT) and the first-order shear deformation plate theory (FSDPT) were used for analysis of bonded joints with laminated composite adherends by Mortensen and Thomsen [25], Wah [28] and Yang and Pang [29]. These applications shed light on whether the analysis of bonded joints can be successfully carried out with ordinary displacement-based models which are the tools customarily employed for the analysis of laminated structures. Nevertheless the adhesive peel stress is disregarded in the constitutive equations and the transverse shear stress is neglected or assumed constant across the thickness (as a consequence, equilibrium is not satisfied at the interfaces and the stress free conditions at the end of the overlap are not met), these studies have shown that the results of displacement-based models can be in a good agreement with those of finite element models, as shown e.g. in [29]. Accordingly further research is motivated in this field because even better results could be obtained without the mentioned infringements. Moreover, displacement-based AM are of interest for analysis of bonded joints because the same approach can be used across the overlap and far from it, so designers can carry out a realistic analysis of the structures in the bonded region and outside contemporaneously and with the same tool. As a contribution towards refinement of studies with displacement-based approaches, in this paper the analysis of bonded joints with laminated adherends is carried out assuming as structural model a zig-zag model with variable kinematics and a fixed number of functional degrees of freedom (d.o.f.). It constitutes a refined version of the layerwise model recently developed by Icardi and Sola [30] and successfully applied to the analysis of laminated and sandwich composites in Ref. [31], which is here reformulated in order to treat bonded joints through inclusion of a variable in-plane representation. The paper is structured as follows. First, the modelling assumptions are discussed and justified, along with loading and boundary conditions requirements. Then, numerical results for benchmark test cases are presented and discussed.

2. Structural model As a necessary premise to the discussion of the structural model, the simulation of joints is briefly overviewed at the light of the layerwise effects that rise in the overlap and in the multilayered structures constituting the adherends, along with the loading and boundary conditions that should be fulfilled. 2.1. Layerwise effects A variety of models with a progressively refined predictive capability have been developed in the past years using different techniques, as discussed above. Bonded joints with laminated adherends need appropriate displacement-based models whose representation satisfies the kinematic and stresses boundary conditions and the local indefinite equilibrium equations at any point. Due to this latter condition, the out-of-plane stresses should be continuous at the material interfaces. As a consequence, the displacements should be represented as piecewise continuous functions having appropriate discontinuous derivatives at the layer interfaces. These models are largely adopted for analysis of composites, while they are only seldom employed for studying joints, since the stress boundary conditions cannot be easily enforced. On the contrary, stress based and mixed models that assume the displacements separately from the stresses are widespread because they facilitate the enforcement of the boundary

246

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

conditions. However, irrespectively of the model used, a crucial feature is whether the number of the functional d.o.f. is fixed, as in the majority of displacement-based models, or whether it varies with the number of layers of the adherends and with the geometry of the joint. From a practical viewpoint it is important to evaluate the performances of structural models at the light of accuracy and costs. Usually, the accuracy of the models with variable d.o.f. is generally better, since their representation can be refined across the thickness in the regions with step gradients. But, unfortunately, their memory storage dimension and processing time can be too large for an extensive application in the industrial environment. On the other hand, the models with fixed d.o.f. are efficient but not always accurate. An alternative way for refining the representation across the thickness without consistently increasing the computational effort is to use a displacement-based zigzag model with adaptive representation and fixed d.o.f., like the model of Ref. [30]. This model was proven capable of capturing the stress fields from constitutive equations even when the constituent layers have distinctly different properties (see, e.g. [30,31]). Here a refined version of the model [30] that enables an in-plane variable representation of displacements is developed and applied to the analysis of bonded joints with laminated adherends. The representation can be different in the adherends and in the overlap, with the aim of accurately simulating the variation of solutions across the joint and to predict the effects of such variation on global and local behaviour, its relevance having been shown by Gaudiello and Kolpakov [32] and Kolpakov [33]. With the present model, a refinement of the representation across the thickness does not mean an increased number of d.o.f. because the expressions of the higher-order coefficients incorporated in the displacement fields are calculated once for all in closed form using symbolic calculus (MATLABs) through the enforcement of all physical constraints. The cost to pay is the presence of derivatives of the d.o.f. in the model, since the enforced constraints usually turn into differential equations. The presence of these derivatives, which is unwise for the development of finite elements since they should appear as nodal d.o.f., thus requiring using computationally inefficient interpolation function, is overcome employing e.g. the strain energy updating technique [34].

expressions of the higher-order contributions to displacements. The relations being rather lengthy and tedious are obtained using symbolic calculus. A different representation for the adherends and for the overlap area can be used, because geometry and material properties change across the joint. In this way always appropriate expressions of the coefficients of displacements are obtained together with the possibility of refining the model in the in-plane direction. If one assumes that no variation occurs in the transverse direction y as customarily, the transverse shear and the normal stresses should identically vanish over the free surface Ω1 and Ω2 of the overlap (see Fig. 1). Therefore it is possible to enforce the following relations to hold in integral or point form: R R sxz ¼ 0; Ωi sxz dz ¼ 0; ∑ Ωk sxz dz ¼ 0; Z Zk sxx ¼ 0; sxx dz ¼ 0; ∑ sxx dz ¼ 0; ð1Þ Ωi

k

Ωk

The stresses are consistent with the applied loads as R R R R N ¼ Ωi sxx dz; M ¼ Ωi zsxx dz; Q ¼ Ωi sxz dz; N ¼  Ωi sxz dx; ð2Þ N, M, Q being the in-plane, bending and shear resultants. All these conditions are used for computing an appropriate number of coefficients of displacements. Non-vanishing stresses can be enforced, because in real adhesive joints a fillet of surplus adhesive, the so-called spew-fillet can be formed at the end of overlap zone allowing transferring the shear stress. Any other condition used in the reference papers considered for comparisons can be enforced; even a non-vanishing transverse shear can be imposed at clamped edges, in order to avoid the poor behaviour of the models with mid-plane displacements and shear rotations as functional d.o.f due to the simultaneous vanishing of displacements and shear rotations. Moreover, the model can be made consistent with a state of zero transverse shear stress in presence of nonzero bending strain in the thin regime and with a state of nonzero transverse normal stress with a nonzero bending strain in the thick regime. 2.3. Preliminaries and notations

2.2. Boundary conditions When adhesively bonded joints with laminated adherends are considered, the problem becomes a multiple-point boundary value problem. The boundary conditions at the edges of the adherends and of the overlap regions are satisfied finding appropriate

Since the adherends may have a different number of layers, a laminated model is developed whose number of layers is different across the joint, as shown in Fig. 1. The middle surface Ω of the overlap is chosen as the plate reference surface. A rectangular Cartesian reference frame (x, y, z)

Fig. 1. Geometry and reference frame for analysis of: (a) single lap joints; (b) double lap joints.

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

with (x, y) in Ω and z normal to it is assumed as reference system. The effects of non-linear strains, which are important for SLJ and DLJ, are accounted for using the updated Lagrangian methodology, which offers an improved numerical efficiency than the standard Lagrangian approach. Material non-linearity can be accounted for by implementing stiffness coefficients that varies with the load magnitude, as successfully done in Ref. [35].

syz ju ¼ 0 sxz jl ¼ 0

ð13Þ

szz ju ¼ p0 ju szz jl ¼ p0 jl

ð14Þ

szz;z ju ¼ 0 szz;z jl ¼ 0

ð15Þ

Also the equilibrium condition is imposed at discrete points across the thickness. ðaÞ sxx;x þ sxy;y þ sxz;y ¼ 0

ðbÞ sxy;x þ syy;y þ syz;z ¼ 0

2.4. Kinematics

uðx; y; zÞ ¼ U 0 ðx; y; zÞ þ U i ðx; y; zÞ þ U c ðx; y; zÞ þ U c_ip ðx; y; zÞ i

c

vðx; y; zÞ ¼ V ðx; y; zÞ þ V ðx; y; zÞ þV ðx; y; zÞ þ V

c_ip

ðx; y; zÞ

wðx; y; zÞ ¼ W 0 ðx; y; zÞ þ W i ðx; y; zÞ þW c ðx; y; zÞ þ W c_ip ðx; y; zÞ

ð3Þ ð4Þ ð5Þ

Please note that all these physical constraints are enforced using a symbolic calculus tool, thus avoiding the drawback of cumbersome algebraic manipulations. The expressions of the terms Δc are chosen as follows: ni

ni

k¼1

k¼1

ni

ni

k¼1

k¼1

U c ðx; yÞ ¼ ∑ Φkx ðx; yÞðz  zk ÞH k þ ∑ C ku ðx; yÞH k V c ðx; yÞ ¼ ∑ Φky ðx; yÞðz  zk ÞH k þ ∑ C kv ðx; yÞH k ni

0

i

c

The first three contributions, ( ) , ( ) and ( ) , are the same as in the model of Ref. [30]. The terms with superscript 0, which are here indicated as Δ0, have a polynomial representation with a fixed expansion order across the thickness. The terms Δi with the superscript i can vary from point to point across the thickness, thus enabling the model to obtain a variable kinematic across the thickness and to fulfil the boundary conditions. Actually, any possible set of boundary conditions can be satisfied, because these terms provide a variable order of representation across the thickness. The terms Δc with the superscript c are aimed at a priori fulfilling the continuity conditions as prescribed by the elasticity theory for keeping equilibrium at the interfaces. As shown in Ref. [30], their expressions are determined once for all for any lay-up. The last contributions Δc_ip enable the fulfilment of the displacement and continuity conditions at the interfaces of adjacent regions across which the in-plane representation is varied. Details are given hereafter. The Δ0 terms have the following expressions that repeat the kinematics of FSDPT model: U 0 ðx; y; zÞ ¼ u0 ðx; yÞ þ zðγ 0x ðx; yÞ  w0;x Þ 0

V ðx; y; zÞ ¼ v 0

0

ðx; yÞ þ zðγ 0y ðx; yÞ w0;y Þ 0

W ðx; y; zÞ ¼ w ðx; yÞ

ð6Þ ð7Þ ð8Þ

They contain just the primary, starting contributions expressed in terms of the functional d.o.f. The terms Δi are postulated to vary across the thickness as follows: U i ðx; y; zÞ ¼ Ax2 ðx; yÞz2 þ ⋯ þ Axn ðx; yÞzo

ð9Þ

V i ðx; y; zÞ ¼ Ay2 ðx; yÞz2 þ ⋯ þ Ayn ðx; yÞzo

ð10Þ

W i ðx; y; zÞ ¼ Az ðx; yÞz þ ⋯ þ Azn ðx; yÞzo

ð11Þ

so to adapt to the variation of solutions, in order to obtain accurate stress predictions directly from constitutive equations. Their expressions are determined imposing the conditions of Eqs. (1) and (2) and, the boundary conditions prescribed by the elasticity theory:

sxz ju ¼ 0 sxz jl ¼ 0

ð16Þ

ðcÞ sxz;x þ syz;y þ szz;z ¼ 0

The through-the-thickness variation of displacements across the thickness is postulated in the following general, piecewise form that considers four separated contributions, as indicated by the superscripts:

0

247

ð12Þ

ð17Þ

ð18Þ

ni

W c ðx; yÞ ¼ ∑ Ψ k ðx; yÞðz  zk ÞH k þ ∑ Ωk ðx; yÞðz  zk Þ2 H k k¼1 ni

þ ∑

k¼1

k¼1

C kw ðx; yÞH k

ð19Þ

here Hk is the Heaviside unit step function that activates the contribution of the continuity functions from the pertinent interface k. The contributions of Eqs. (17)–(19) are aimed at fulfilling the continuity constraints at the material interfaces in the thickness direction:

sxz jðkÞz þ ¼ sxz jðkÞz  syz jðkÞz þ ¼ syz jðkÞz  sz jðkÞz þ ¼ sz jðkÞz  sz;z jðkÞz þ ¼ sz;z jðkÞz  ujðkÞz þ ¼ ujðkÞz  vjðkÞz þ ¼ vjðkÞz  wjðkÞz þ ¼ wjðkÞz 

ð20Þ c_ip

aimed at fulfilling the continuity conditions in The terms Δ case of variation of the in-plane representation have the following expressions: U c_ip ðx; yÞ ¼ Θkx ðx; yÞðx  xk ÞH k þ Γ kx ðx; yÞðx  xk Þ2 H k

ð21Þ

V c_ip ðx; yÞ ¼ Θky ðx; yÞðx  xk ÞH k þ Γ ky ðx; yÞðx  xk Þ2 H k

ð22Þ

W c_ip ðx; yÞ ¼ ðΘkw þ Θkw* þ χ kw þ χ kw* Þðx; yÞðx  xk ÞH k þ ðΓ kw þΓ kw* þ Λkw þ Λkw* Þðx; yÞðx  xk ÞH k

ð23Þ

As for Eqs. (17)–(19), Hk activates the contribution of the continuity functions starting from the variation of the in-plane representation. It could be noticed that this variation can happen only in one direction (x or y) and not in both at the same time. The piecewise contributions of Eqs. (21)–(23) are incorporated in order to make continuous the stresses and their gradients up the desired order at the interfaces of the regions where the in-plane representation is changed, as prescribed by the elasticity theory. Generally speaking, the summations of Eqs. (21)–(23) can have infinite terms, since from equilibrium it is possible to deduce that the continuity should hold for the stresses and their in-plane gradient of any order. By the practical viewpoint, the summations can be truncated, without significant loss of accuracy. Accordingly, in the present paper the summations of Eqs. (21)–(23) are considered, without affecting accuracy as shown by the results

248

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

Fig. 2. Stages of the procedure employed to compute the stress field of the joint.

of Section 3. In details, the terms of first order in the in-plane coordinate are aimed at fulfilling the continuity of the stresses; those of second order are aimed at satisfying the continuity of the stress gradients. Like for the former terms, the explicit expressions of the continuity functions are evaluated in closed form using a symbolic calculus tool. 2.5. Solution methodology The partial differential governing equations are obtained in a straightforward way from the present model using standard techniques (Principle of Virtual work, Rayleigh–Ritz method, etc.) as shown in Fig. 2. Their expressions are here omitted being rather lengthy. The solution to these equations is found in closed form by the Fourier series method assuming the displacement fields as:     M N mπ nπ u0 ðx; yÞ ¼ ∑ ∑ Amn cos x sin y ; Lx Ly m ¼ 1n ¼ 1     M N mπ nπ x cos y ; v0 ðx; yÞ ¼ ∑ ∑ Bmn sin Lx Ly m ¼ 1n ¼ 1     M N mπ nπ x sin y ; w0 ðx; yÞ ¼ ∑ ∑ C mn sin Lx Ly m ¼ 1n ¼ 1     M N mπ nπ x sin y ; γ x 0 ðx; yÞ ¼ ∑ ∑ Dmn cos Lx Ly m ¼ 1n ¼ 1     M N mπ nπ x cos y ; ð24Þ γ y 0 ðx; yÞ ¼ ∑ ∑ Emn sin Lx Ly m ¼ 1n ¼ 1 Substituting the former representation into the governing differential equations, a set of algebraic equations is obtained that gives the unknown amplitudes Amn to Emn. Then, the stress field can be obtained from the stress–strain and strain–displacement relations. Eq. (24), which represent a general solution for general loading and boundary conditions and either linear or nonlinear strains, were used by Yang and Pang [29] for solving SLJ subdividing into three zones two of which being outside the overlap and the other being the overlap itself. Many other researchers used this

method for solving problems with non-classical, intricate boundary conditions, like those considered hereafter to assess the model. As shown in the literature, convergent results are obtained for these cases using at least a hundred of terms, so processing time and memory occupation dimension are much larger than those required for solving problems with conventional boundary conditions and simple loading schemes. Please note that in twodimensional problems of SLJ v0 and γ0y are neglected since the solution does not vary in the transverse direction y, in accordance with the reference articles used for comparisons. In such cases, the following fast alternative representation of the displacements is also used:     πx πx u0 ðx; yÞ ¼ AU ½ coshðxÞ  1 cos ; w0 ðxÞ ¼ AW sin ; Lx Lx   πx γ x 0 ðx; yÞ ¼ Aγ x ½ coshðxÞ cos ; ð25Þ Lx that just requires the solution of three algebraic equations in the three unknowns Au , Aw , Aγ x . The comparison with analytical and finite element solutions by other researchers will show that quite accurate results can be obtained with such single component expansions at least for the sample cases available in the literature, with a self-evident advantage by the viewpoint of computational costs.

3. Numerical applications and discussion Sample cases of SLJ and DLJ published in the literature are retaken and analysed by the present model, in order to assess whether it can be successfully employed for studying adhesively bonded joints, and thus whether a unified, efficient approach can be used for treating laminates, as well as the joints used to connect them. In order to preliminary assess the capability of the present analytical model to predict the stress fields across bonded components and to treat general boundary conditions, which is a requirement of primary importance at the edges of the overlap, the results for a piezoactuated beam made of an aluminium

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

249

substrate, a bonding film and a piezoactuator and for a plate with clamped edges are premised to those for SLJ and DLJ. All the computational times provided in the following sections are obtained carrying out the analyses on a laptop computer with dual-core CPU 2.20 GHz, 64 bit operating system and 4 GB RAM. 3.1. Piezoelectrically actuated beam Consider the sample case studied by Robbins and Reddy [36] using 3-D finite elements derived from a displacement-based layerwise theory, which consists of a three layered cantilever beam 152 mm long and made of an underlying aluminium beam substructure, an adhesive film and a piezoactuator bonded on the upper face. Piezoactuator, adhesive and aluminium substructure form an unsymmetrical three material laminate that exhibits bending/extension coupling. For this sample case, whose stress fields are similar to those of bonded joints, the only acting loads are the self-equilibrating loads induced by the piezoactuator. A bending deformation is provided by applying an actuation strain of 0.001 to the piezoelectric layer via an applied electric field. Of course, the main interest here is to emphasize the implications of the localized effects induced by the piezoactuator control layers on through-the-thickness stress distribution across the adhesive layer. The comparison with the reference FEA 3-D results will be used to assess the effectiveness of the present model in treating the intricate stress fields in the most critical regions typical of adhesively bonded structures and their non-classical boundary conditions. The material properties and the thickness (t) of the three constituent layers used in the computations are the following: aluminium substrate: E¼69 GPa, G¼27.5 GPa, υ¼0.25, t¼ 15.2 mm; adhesive: E¼6.9 GPa, G¼2.5 GPa, υ¼0.4, t¼0.254 mm; piezoactuator: E1 ¼ 69 GPa, E3 ¼48 GPa, G¼21 GPa, υ13 ¼0.25, υ31 ¼0.175, t¼1.52 mm. As shown by Robbins and Reddy [36], unwanted dangerous stress concentrations take place at the free edge. As for bonded joints, close to the free edge the transverse normal and shear stresses reach a very large peak nearly in the adhesive layer. Owing to these interlaminar stress concentrations, the adhesive layer may progressively fail till to complete debonding of the piezoactuator. Fig. 3a shows the spanwise variation of membrane, transverse shear and transverse normal stresses near the top of the aluminium substrate (z1 ¼6.61 mm) and at the centre of the adhesive layer (z2 ¼6.84 mm), while Fig. 3b represents the stress fields across the thickness close to the free edge. The stress distributions are presented in the following normalized form, according to [36]: Atot
sx
103 sx ¼ ðE1 AÞalum þ ðE1 AÞad þ ðE1 AÞpiezo

sz ¼

Atot
sz
103 ðE3 AÞalum þ ðE3 AÞad þ ðE3 AÞpiezo

sxz ¼

Atot
sxz
103 ðG13 AÞalum þ ðG13 AÞad þ ðG13 AÞpiezo

ð26Þ

the subscripts alum, ad, piezo and tot being used to indicate the cross sectional area and the elastic moduli of the substrate structure, of the adhesive and of the piezoactuator layer, respectively. In this case, 6 subdivisions are considered across the thickness of the aluminium substrate for computing the coefficients of displacements and the continuity functions of the present model, while 5 subdivisions are used across the adhesive layer and 9 in the piezoactuator layer. These computational subdivisions are refined near the adhesive layer and gradually enlarged as the distance to this layer increases.

Fig. 3. Piezoactuated beam. Comparison beetween: (a) the span-wise distributions and (b) the through-the-thickness stress distribution close to the free edge, as predicted by the present analytical model and by Robbins and Reddy [36] via 3D FEA.

The trial functions at the tip and the root of the beam are selected as follows in order to satisfy the boundary conditions, R i.e. w(0)¼ w,x(0)¼0 and zsxx dz ¼ 0:   x  x i þ 3 w A ð27Þ w0 ðxÞ ¼ 1  ði þ 3Þ  1  L L retaining only the first component, i.e. i ¼1, while the in-plane displacement and the shear rotation are obtained from (27) deriving in x and considering Au and Aγ x as amplitudes. Hermite's cubic polynomials are used elsewhere for all the components: Aa H a ðxÞ þ Ab H b ðxÞ þ Ac H c ðxÞ þ Ad H d ðxÞ

ð28Þ

The predictions being in good agreement with those of Ref. [36] and the boundary conditions being correctly met in a point form, as shown in Fig. 3, the present model is proven to accurately carry out the analysis of adhesively bonded structures. Nevertheless many subdivisions are used across the thickness and in the spanwise direction, just a couple of minutes were required for computing the coefficients and solving the problem (i. e. less than for a FEA 3-D analysis carried out by a commercial code), thus the present model is also efficient. As a further test, a plate with clamped edges is considered, since the models based on the classical middle-plane d.o.f. were proven to give poor results in the literature for this kind of structures. 3.2. Plate with clamped edges The case of a square plate that is simply supported on two opposite edges and clamped on the other two is considered, whose top surface is subjected to a bi-sinusoidal normal load, whereas the bottom surface is traction free. Vel and Batra [37]

250

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

exactly solved the governing equations using a three-dimensional generalization of the Eshelby–Stroh formalism that determines an infinite system of equations in infinite unknowns. The boundary conditions at the edges and the continuity conditions at the interfaces were used to define constants in the general solution, which were determined by the Fourier series method. Truncations of the set of infinite equations introduce errors that are minimized by increasing the number of the terms in the series. The solution to this case is obtained by the present model assuming the following in-plane representation of the displacements:     mπ nπ x cos y ; Lx Ly m ¼ 2;4n ¼ 1;3     M N mπ nπ x sin y ; v0 ðx; yÞ ¼ ∑ ∑ Bmn cos Lx Ly m ¼ 1;3n ¼ 2;4        M N mπ nπ ∑ C mn cos x  ð 1Þm=2 cos y  ð  1Þn=2 ; w0 ðx; yÞ ¼ ∑ Lx Ly m ¼ 2;4n ¼ 2;4 u0 ðx; yÞ ¼

γ 0x ðx; yÞ ¼ γ 0y ðx; yÞ ¼

M

∑

M

∑

N

∑ Amn sin

N

m ¼ 2;4n ¼ 1;3 M

∑

    mπ nπ x cos y ; Lx Ly     mπ nπ cos x sin y ; Lx Ly

∑ Dmn sin N

∑ Emn

m ¼ 1;3n ¼ 2;4

ð29Þ

Owing to the adaptive terms incorporated in the present model, a non-vanishing transverse shear can be enforced with displacements and shear rotations that vanish at clamped edges, along with the stress-free boundary conditions at the free edges. Because the stress boundary conditions at the edges of the overlap of bonded joints are enforced in a similar way, the comparison with the results by Vel and Batra [37] will serve as a test for assessing whether the present model can treat general boundary conditions. As in [37], a length-to-thickness ratio of 5 is considered, assuming that each lamina has the following mechanical properties: EL/ET ¼25; GLT/ET ¼ 0.5; GTT/ET ¼ 0.2; υLT ¼0.25 and the stacking sequence is assumed as [01/901/01]. The through-the-thickness variation of the transverse displacement component at the centre of the plate w and the transverse shear stress sxz at the clamped edge are presented in normalized form as follows, according to [37]:  3 s  Lx þ 0:05Lx ; 0; z 100ET
h sxz ¼ xz 2 w¼ wð0; 0; zÞ ð290 Þ sxz max p0 L4x Lx being the side length, h the thickness of the plate and p0 the amplitude of the bi-sinusidal load. Fig. 4a reports the variation of the transverse shear stress across the thickness predicted by the present model and by Vel and Batra [37] with 250 terms in the series. The present results with dashed lines are computed considering 200 terms, while the ones with a solid line are obtained considering just a single component in the series expansion, i.e. M¼ N ¼1. Fig. 4b shows how the transverse displacement varies across the thickness. In this case, the comparison results are obtained from those presented in [37] in tabular form. It could be seen that the present model accurately predicts the through-the-thickness variation of solutions even when just a single term representation is adopted. Thus it provides accurate results with a low computational effort, with self-evident advantages. This capability is a direct consequence of the a priori fulfilment of the stress and displacement boundary conditions and of the interfacial stress contact conditions obtained determining expressions in closed form for the coefficients of displacements. The computation of these coefficients, which is required just once at a time, takes about a minute, while the solution of

Fig. 4. Plate with two simply supported edges and two clamped edges. Comparison beetween the through-the-thickness distribution of: (a) transverse shear stress and (b) transverse displacement by the present and Vel and Batra [37] analytical models.

governing equations takes four minutes with M ¼N ¼ 200 and just 1.3 s with M¼N ¼ 1. The capability of the model of being accurate with M ¼N ¼ 1 obviously gives a considerable practical advantage when repeated computations are required, as in non-linear problems. 3.3. Adhesively bonded joints The sample cases by Diaz Diaz et al. [6], Radice and Vinson [7], Andruet et al. [11] Vaidya et al. [27] and da Silva and Adams [38] are considered in order to assess the capability of the present displacement-based multilayered plate model of describing the stress fields in adhesively bonded joints. 3.3.1. Case A: Single-lap joint with aluminium adherends As first application, the sample case by Radice and Vinson [7] is considered. In this study, a single-lap joint clamped at the left edge and simply supported at the right one undergoing an axial load with intensity 40 kN/mm is analysed. The adherends are made of an aluminium alloy (E ¼70 GPa, υ ¼0.33), while the adhesive is made of an epoxy resin (E¼ 4.82 GPa, υ ¼0.4). The adherends are 50.8 mm long and 1.62 mm thick, the overlap length is 12.7 mm, while the adhesive thickness is 0.25 mm. It could be remarked that the length- to-thickness ratio of the adherends is more than 30 and, in addition, they are made of isotropic material, therefore the present model does not require the computation of the continuity functions in the adherends, nor an extremely high order expansion across the thickness. Accordingly, a third order expansion is adopted for the in-plane displacement and a fourth order one for the transverse displacement. Fig. 5a shows how the shear stress sxz varies along the bond-line,

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

Fig. 5. Single-lap joint with aluminium adherends. Span-wise distribution of: (a) shear stress and (b) peel stress at the mid-plane of the adhesive layer by Radice and Vinson's elasticity model [7] and by the present model using different shape functions.

whereas Fig. 5b gives the variation of the peeling stress sz along the same line. According to [7], the variation of solutions in the transverse direction y is neglected. The results by the present model reported as a solid line are obtained using as trial functions those of Eq. (25) in the overlap section. In this way, the problem turns into a system of three equations in the three unknown amplitude of the in-plane and transverse displacement components and of the transverse shear rotations. The solution represented as a dashed line is obtained using a Fourier series expansion of Eq. (24) with 30 components in the overlap section. In both cases, the trial functions of the lower adherend are those of Eq. (24), while those of Eq. (24) are chosen into the upper adherend, in order to respect the boundary conditions. In the former case, the problem is solved in 1.5 s, while 80 s are required in the latter one. It could be noticed that also in this benchmark case the present model is proven to be able of capturing the stress fields with just one component in the series, so it appears suitable to efficiently solve the problem examined.

3.3.2. Case B: Single-lap joint with CFRP adherends As next case, the single-lap joint with CFRP adherends analysed by Vaidya et al. [27] and retaken by Khalili et al. [12] is considered. The mechanical properties of the adherends are: E1 ¼ E2 ¼50 GPa, E3 ¼7.2 GPa, G12 ¼5 GPa, G13 ¼ G23 ¼3 GPa, υ12 ¼0.3, υ13 ¼ υ23 ¼0.25. Those of the adhesive are: E¼ 2.8 GPa, υ ¼0.35. The adherends are 87.5 mm long and 3 mm thick, the adhesive thickness is 0.33 mm, while the overlap length is 25.4 mm. The lower adherend is restrained against translation but free to rotate, while the other one is simply-supported at the edge, where a uniformly

251

Fig. 6. Single-lap joint with CFRP adherends. Span-wise distribution of: (a) shear stress and (b) peel stress by Vaidya et al. [27] (FEA) and by the present model enforcing or not the free-edge stress boundary conditions.

distributed traction load with intensity 7500 N is applied. Fig. 6a and b show the span-wise distributions of transverse shear stress and peel stress, as predicted by the present model and in [27]. Similarly to what done in the previous case, the variation of solutions in the transverse direction y is still neglected. Please note that the curves named BC1 are computed imposing the stress-free boundary condition at the edges of the adhesive, those referred as BC2 are instead achieved imposing the same value of the shear as Vaidya et al. [27] at the edge of the adhesive. In both cases, the trial functions for the adherends are those of Eq. (24) while in the overlap the trial functions are those of Eq. (25). As far as the through-the thickness representation is concerned, in both cases a third-order representation of the in-plane displacement and a fourth-order one for the transverse displacement is adopted. Differently to Case A, two computational layers are considered across the thickness of the adherends, in order to accurately represent their stress variation. Nevertheless more computational layers are considered than in the former case, the processing time is still comparable, since just 1.7 s are required for computing the coefficients. Also in this case the representation with a single component gives quite accurate result with the lowest computational effort.

3.3.3. Case C: Double-lap heterogeneous joint We now consider the case analysed by Diaz Diaz et al. [6] of a double-lap joint with adherends made of different materials. The outer adherends are assumed to be made of steel (E¼ 210 GPa, υ¼ 0.3), while the inner adherend is assumed to be orthotropic (E1 ¼120 GPa, E3 ¼ 10 GPa, G13 ¼ 5 GPa, υ 13 ¼0.3). The adhesive is an epoxy resin (E¼ 4 GPa, υ ¼0.3). The thickness of the adherends

252

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

Fig. 7. Double-lap joint with heterogeneous adherends. Span-wise distribution of shear stress between (a) lower adherend and adhesive and (b) inner adherend and adhesive and of peel stress between (c) lower adherend and adhesive and (d) inner adherend and adhesive by Diaz-Diaz et al. [6] (analytical and FEA) and by the present model using different trial functions.

is 3 mm, while that of the adhesive is 0.5 mm. The overlap length is 40 mm. The joint is undergoing a traction load of 5000 N. A comparison is presented with the results by the structural model of Ref. [6] and by a finite element model therein considered for comparisons. The analysis with the present model is carried out considering a second order expansion of the in-plane displacement and a third order expansion of the transverse displacement across the thickness of the outer adherends and of the adhesive. Instead, a third order expansion of the in-plane displacement and a fourth order expansion of the transverse displacement across the thickness are adopted within the inner orthotropic adherend. As no information about the total length of the adherends and about their boundary conditions can be gathered from Diaz Diaz et al. [6], only the overlap section is simulated by the present model. Fig. 7 shows the span-wise distribution of transverse shear and transverse normal stresses between lower outer adherend and adhesive and between inner adherend and adhesive, as predicted by the present model using the single component representation of Eq. (25), i.e. M,N ¼1, and the series expansion of Eq. (24) with M, N ¼30. As it takes only 1.3 s, obviously the first approach is dramatically faster, since 65 s are required to perform the analysis with 30 components. The comparison with the reference results shows that the fastest approach is also quite accurate and thus efficient, considering that the use of many components in the series does not produce a consistent accuracy improvement, while it considerably increases the processing time and the memory occupation.

3.3.4. Case D: Single-lap joint taking into account non-linearity In order to verify whether the method employed to model nonlinearity provides accurate results, we now consider the single-lap joint studied by Andruet et al. [11]. The joint has aluminium adherends (E ¼68.3 GPa, υ ¼0.3) and epoxy resin adhesive (E ¼2.5 GPa, υ¼0.3); it is clamped at left edge and simply supported at the right one, where an average tensile stress with intensity 200 MPa is applied. The thickness of the adherends is 1.6 mm, that of the adhesive is 0.1 mm. The overlap length is 20 mm, while the adherends are 60 mm long. Fig. 8 shows the results obtained by the present variable kinematics model either imposing the boundary constraints prescribed by the theory of elasticity, which are named as BC1, or imposing the same value of the shear stress computed by Andruet et al. [11] at the free edge of the adhesive, here named as BC2. In both cases, a third order piecewise variation of the in-plane displacement and a fourth order variation of the transverse displacement is chosen, using as trial functions those of Eq. (25). Irrespectively of the constraints imposed, the computational time is 1.5 s. As far as the in-plane representations of displacements in the adherends are concerned, the same trial functions already adopted for Case A are employed. The results show that also in this case the present model provides predictions in a very good agreement with the reference ones. This confirms the correct implementation and the effectiveness of the updated Lagrangian approach in taking into account the effects of geometric nonlinearity, as already pointed out in the literature.

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

253

experiments [38] confirm that the present model accurately predicts the stresses and, besides, that stress-based criteria like [39,40] can be successfully used for determining the strength of bonded joints.

4. Concluding remarks

Fig. 8. Non-linear analysis of single-lap joint. Span-wise distribution of: (a) shear stress and (b) peel stress by Andruet et al. [11] (FEA) and by the present model enforcing or not the free-edge stress boundary conditions.

Table 1 Comparison between failure loads [kN] predicted by the present model and those experimentally measured in Ref. [38] for titanium/composites double lap joints. (FC ¼Failure criterion). Composite inner adherend Supreme 10HT

AV119

Composite outer adherend

Da Silva and Adams 24.40 [38] Present (FC [39]) 26.20 Present (FC [40]) 26.35

38.54 ( 71.21)

Da Silva and Adams 21.89 (7 1.11) [38] Present (FC [39]) 23.48 Present (FC [40]) 23.15

37.33 ( 7 3.93)

40.50 40.30

39.19 39.51

3.3.5. Case E: Failure load prediction for a double lap joint In order to compare the numerical results of the present model to those obtained through experimental tests, the failure loads of titanium/composite double lap joints presented by da Silva and Adams [38] are considered. The adherends are made of titanium (Ti–6Al–4V) and of a 01/901 laminate obtained by staking two plies of carbon fabric HTM552, impregnated with a BMI resin composite. According to [38], different adhesives are considered, namely Supreme 10 HT or AV119. Table 1 reports the comparisons between experimental and numerical results obtained assuming for the d.o. f. the in-plane representation reported in Eq. (24) considering M ¼N ¼35. For what concerns the representation across the thickness, a third order piecewise variation of the in-plane displacement and a fourth order variation of the transverse displacement is chosen. The failure loads reported in Table 1 are computed considering the stress-based failure criteria of Refs. [39] and [40]. The good agreement shown with respect to the results of

A refined zig-zag model with variable kinematics was developed to study adhesively bonded joints with laminated adherends, thus enabling the analysis of laminates and of bonded joints under a unified approach. Its piecewise, three-dimensional displacement field a priori satisfies the out-of-plane stress contact conditions at the interfaces of adjacent layers and the stress boundary conditions at the upper and lower faces of the structure. The representation can vary from point-to point across the thickness, in order to capture step gradients, while the in-plane representation can vary from the adherends to the overlap in order to better simulate the variation of solutions and to satisfy the stress boundary conditions of the joint. However no increase of the memory storage dimension and no consistent increase of the processing time are involved varying the representation, because just five functional d.o.f. are used, like with classical plate model. The relations that follow from the enforcement of the continuity of out-of-plane stresses at the interfaces of adjacent layers and of the boundary conditions are obtained in closed form as expressions of the d.o.f. and of their derivatives through the use of symbolic calculus. As a preliminary assessment, a piezoactuated beam made of an aluminium substrate, a bonding film and a piezoactuator was considered in order to test the capability of the present model to capture the out-of-plane stress fields across the thickness of the adhesive. In order to assess the capability of the present analytical model to treat general boundary conditions, like at the edges of the overlap, a plate with clamped edges was analysed. For both benchmark cases accurate results were obtained, as shown by the comparisons with the reference three-dimensional results in the literature, in few seconds. Applications to linear and nonlinear benchmark test cases of joints taken from the literature were presented. Numerical results show that the model here presented can accurately and efficiently treat bonded joints with laminated adherends. In fact, for all the examined cases, accurate results were obtained in few seconds carrying out the analyses on a laptop computer. References [1] Her S-C. Stress analysis of adhesively-bonded lap joints. Comput Struct 1999;47:673–8. [2] Kuno, JK., Structural adhesives and bonding. In: Proceedings of the structural adhesives and bonding conference, El Segundo; 1979. [3] Vinson JR. Adhesive bonding of polymer composites. Polym Eng Sci 1989;29:1325–31. [4] da Silva LFM, das Neves PJC, Adams RD, Spelt JK. Analytical models of adhesively bonded joints-Part I: Literature survey. Int J Adhes Adhes 2009;29: 319–30. [5] Gustafson PA, Bizard A, Waas AM. Dimensionless parameters in symmetric double lap joints: an orthotropic solution for thermo-mechanical loading. Int J Solids Struct 2007;44:5774–95. [6] Diaz Diaz A, Hadj-Ahmed R, Foret G, Ehrlacher A. Stress analysis in a classical double-lap, adhesively bonded joint with a layerwise model. Int J Adhes Adhes 2009;29:67–76. [7] Radice JJ, Vinson JR. On the analysis of adhesively bonded structures: a higher order semi-elastic adhesive layer model. Compos Sci Technol 2008;68:376–86. [8] Wu XF, Jenson RA. Stress-function variational method for stress analysis of bonded joints under mechanical and thermal loads. Int J Eng Sci 2011;49: 279–94. [9] Duong CN. A unified approach to geometrically nonlinear analysis of tapered bonded joints and doublers. Int J Solids Struct 2006;43:3498–526. [10] Schmidt P. Modelling of adhesively bonded joints by an asymptotic method. Int J Eng Sci 2008;46:1291–324.

254

U. Icardi, F. Sola / International Journal of Adhesion & Adhesives 50 (2014) 244–254

[11] Andruet RH, Dillard DA, Holzer SM. Two- and three-dimensional geometrical nonlinear finite elements for analysis of adhesive joints. Int J Adhes Adhes 2001;21:17–34. [12] Khalili SMR, Khalili S, Pirouzhashemi MR, Shokuhfar A, Mittal RK. Numerical study of lap joints with composite adhesives and composite adherents subjected to in-plane and transverse loads. Int J Adhes Adhes 2008;28:411–8. [13] Diaz J, Romera L, Hernàndez S, Baldomir A. Benchmarking of threedimensional finite element models of CFRP single-lap bonded joints. Int J Adhes Adhes 2010;30:178–89. [14] Kilic B, Madenci E, Ambur DR. Influence of adhesive spew in bonded single lap joints. Eng Fract Mech 2006;73:1472–90. [15] Xu W, Li G. Finite difference three-dimensional solution of stresses in adhesively bonded composite tubular joints subjected to torsion. Int J Adhes Adhes 2010;30:191–9. [16] Nemes O, Lachaud F. Double-lap adhesive bonded-joints assemblies modelling. Int J Adhes Adhes 2010;30:288–97. [17] da Silva LFM, das Neves PJC, Adams RD, Spelt JK. Analytical models of adhesively bonded joints-Part II: Comparative study. Int J Adhes Adhes 2009;29:331–41. [18] Volkersen O. Die nietkraftverteilung in zugbeanspruchten nietverbindungen mit konstanten laschenquerschnitten. Luftfahrtforschung 1938;35:4–47. [19] Goland M, Reissner E. The stresses in cemented joints. J Appl Mech 1944; A1:17–27. [20] Hart-Smith, LJ. Adhesive-bonded single-lap joints, technical report NASA CR 112236; 1973. [21] Renton WJ, Vinson JR. The efficient design of adhesive bonded joints. Int J Adhes 1975;7:175–93. [22] Srinivas, S. Analysis of bonded joints, NASA technical note TN D-7855; 1975. [23] Allman DJ. A theory for elastic stresses in adhesive bonded lap joints. Q J Mech Appl Math 1977;30:415–36. [24] Adams RD, Mallick V. A method for the stress analysis of lap joints. Int J Adhes 1992;38:199–217. [25] Mortensen F, Thomsen OT. Analysis of adhesive bonded joints: a unified approach. Compos Sci Technol 2002;62:1011–31.

[26] Smeltzer, SS. An inelastic analysis methodology for bonded joints with shear deformable, anisotropic adherents. PhD dissertation. North Carolina State University, Raleigh NC; 2003. [27] Vaidya UK, Gautam ARS, Hosur M, Dutta P. Experimental-numerical studies of transverse impact response of adhesively bonded lap joints in composite structures. Int J Adhes Adhes 2006;26:184–98. [28] Wah T. Stress distribution in a bonded anisotropic lap joint. ASME J Eng Mater Technol 1973;95:174–81. [29] Yang C, Pang SS. Stress analysis of single-lap composite joints under tension. ASME J Eng Mater Technol 1996;118:247–55. [30] Icardi U, Sola F. Development of an efficient zig-zag model with variable representation of displacements across the thickness. J Eng Mech 2014;140:531–41. [31] Icardi U, Sola F. Recovering critical stresses in sandwiches using through-thethickness reinforcement. Composites Part B 2013;54:269–77. [32] Gaudiello A, Kolpakov AG. Influence of non-degenerated joint on the global and local behaviour of joined rods. Int J Eng Sci 2011;49:295–309. [33] Kolpakov AG. Influence of non-degenerated joint on the global and local behaviour of joined plates. Int J Eng Sci 2011;49:1216–31. [34] Icardi U. C0 plate element based on strain energy updating and spline interpolation, for analysis of impact damage in laminated composites. Int J Impact Eng 2007;34:1835–68. [35] Icardi U, Sola F. Indentation of sandwiches using a refined zig-zag model and a mesoscale damage model. Univers J Mech Eng Educ 2014;2:6–19. [36] Robbins DH, Reddy JN. Analysis of piezoelectrically actuated beams using a layerwise displacement theory. Comput Struct 1991;11:265–79. [37] Vel SS, Batra RC. Analytical solution for rectangular thick laminated plates subjected to arbitrary boundary conditions. AIAA J 1999;11:1464–73. [38] da Silva LFM, Adams RD. Techniques to reduce the peel stresses in adhesive joints with composites. Int J Adhes Adhes 2007;27:227–35. [39] Rotem A, Hashin Z. Failure modes of angle-ply laminates. J Compos Mater 1975;9:191–206. [40] Puck A, Shurmann H. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 1998;58:1045–67.