Assessment of recent zig-zag theories for laminated and sandwich structures

Accepted Manuscript Assessment of recent zig-zag theories for laminated and sandwich structures Ugo Icardi, Federico Sola PII:

S1359-8368(16)30477-2

DOI:

10.1016/j.compositesb.2016.04.058

Reference:

JCOMB 4257

To appear in:

Composites Part B

Received Date: 12 October 2015 Revised Date:

15 March 2016

Accepted Date: 21 April 2016

Please cite this article as: Icardi U, Sola F, Assessment of recent zig-zag theories for laminated and sandwich structures, Composites Part B (2016), doi: 10.1016/j.compositesb.2016.04.058. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

ASSESSMENT OF RECENT ZIG-ZAG THEORIES FOR LAMINATED AND SANDWICH STRUCTURES Ugo Icardi, Federico Sola Dipartimento di Ingegneria Meccanica e Aerospaziale Politecnico di Torino - Corso Duca degli Abruzzi 24, 10129 Torino, Italy

RI PT

E-mail address: [email protected]; [email protected]

Abstract

The flexural response of laminated composite and sandwich beams/plates under static distributed

loading, classical and non-classical boundary conditions (simply-supported, cantilever and proppedcantilever) and geometric/ constitutive heterogeneity of layers is analysed. As structural models, a

SC

recently developed equivalent single-layer zig-zag model by the authors, a discrete-layer model and a sublaminate model developed from it in this paper are used. Their contribution is to consider the

M AN U

continuity of the transverse normal stress and its through-thickness gradient at layer interfaces, as prescribed by the elasticity theory in addition to kinematic and transverse shear stress interlayer continuity customary considered in the literature. To this purpose, a piecewise variation of the three displacement components is adopted. The zig-zag amplitude expressions are obtained in closed-form from the enforcement of stress continuity conditions. To be refined without affecting costs, the equivalent single-layer model has variable kinematics and just five unknowns. Instead, sublaminate and discrete-layer models are refined by increasing the number of computational layers and variables. The

TE D

aim of this paper is to assess whether the equivalent single-layer model having a computational cost comparable to that of classical models can be as accurate ad discrete-layer and sublaminate models. Benchmarks are presented, for which exact elasticity and approximate solutions are available for comparisons in literature. It is illustrated the utility of considering a variable kinematics for obtaining accurate stress predictions from constitutive equations and the transverse normal deformability for

EP

keeping equilibrium. The equivalent single-layer model is shown as accurate as discrete-layer and

AC C

sublaminate models in all cases examined.

Keywords: C Computational modelling; C Numerical analysis; A Layered structures.

1

ACCEPTED MANUSCRIPT Acronyms

3D FEA

Finite element analysis using solid elements

FSDT

First-order shear deformation theory [29]

ADZZ

Present adaptive zig-zag model

HSDT

High-order shear deformation theory [30]

ADL

Present discrete-layer model

KC

ASB

Present sublaminate model

LIN-ADZZ

CUB-ADZZ

Present cubic mixed model

RZT

DL

Discrete-layer plate theories

SC

ESL

Equivalent single-layer plate theory

ZZ-F

F k (z )

Murakami’s zig-zag function [28]

RI PT

SC

M AN U

A2x … J A3z

Models with Murakami’s zig-zag function Di Sciuva’s zig-zag function [27]

Higher-order coefficients of ADZZ

Uc,Vc,Wc

Zig-zag terms of ADZZ

Higher-order coefficients of ADL, ASB

(x,y,z)

Cartesian reference frame

TE D

J

ZZ-K

Models with Di Sciuva’s zig-zag function

( z − zk ) H k

Nomenclature

Ax 2 ... Azn

Kinematic continuity cond. Present linear mixed model Refined zig-zag theory [48] Stress continuity cond.

Upper+ and lower -

Zig-zag amplitudes that restore the continuity of displacements

C11 , C12 ...C66

Elastic coefficients of orthotropic media

zk

Position of the k-th layer

E11, E22 , E33 , G12 , G23 , G31

Elastic moduli of orthotropic media

ε xx , ε xy ...ε zz

Strain components

P0

Distributed transverse load

Φxk, Φyk, Ψk, Ωk

Zig-zag amplitudes of ADZZ

AC C

EP

Cuk, Cvk, Cwk

u,v,w

Elastic displacement components

u0, v0, w0,γx0, γy0

J

0 u 0 , J v 0 , J w 0 J γ x0 , J γ y

U°,V°,W° i

i

U ,V ,W

i

k) + (k) -

z,

z

J

k Φ kx , J Φ y , J Ψ k ,

J

Ω

k

surfaces of k--th layer

Zig-zag amplitudes of ADL, ASB

Functional d.o.f. of ADZZ

υ12 , υ 21 , υ23 , υ32 , υ13 , υ31

Poisson’s ratios

Functional d.o.f. of ADL, ASB

Ω

Reference mid-plane

Basic terms of ADZZ

σ xx , σ xy ...σ zz

Stress components

Higher-order terms of ADZZ

2

ACCEPTED MANUSCRIPT 1 Introduction Laminated and sandwich composites have become increasingly important as primary structures components in many key engineering areas. In aeronautics, where fibre-reinforced and sandwich materials are used since the late 70’s, higher speed, longer range, larger pay loads, less engine power and better operating economy are achieved owing to their excellent structural performances and

RI PT

low weight. This paper deals with the ply-level stress analysis of laminated structures and sandwiches. Three-

dimensional stress states rise in these materials that generate local damaging at layer interfaces as

adhesive bond separation and at the fibre-matrix level as fibre rupture and matrix cracking, as discussed by Reddy[1] and Reddy and Robbins [2]. The damage size may grow under the loadings in service

SC

causing a relevant loss of strength and stiffness. The damage mechanism and how they reflect into a loss of strength and stiffness are extensively discussed by Liu and Islam [3] and Vachon et al. [4]. Modeling of fibre-reinforced laminated and sandwich composites is challenging because the through-

M AN U

thickness displacement field is no longer C1-continuous as in homogeneous non-layered materials, but instead is C°-continuous. The slope of in-plane displacements must change at layer interfaces because strains should be discontinuous for obtaining continuous transverse shear and normal stresses, as required to satisfy local equilibrium equations. The difference in magnitude of transverse shear strains drives the slope discontinuity of the displacement field, which becomes much severe if the difference in transverse shear and normal moduli of layers increase (transverse anisotropy), and links kinematics to stresses. In the current terminology, the slope change is termed as the zig-zag effect.

TE D

The difficulty in accurately modelling the above mentioned layerwise effects is that displacement and through-thickness stress fields are interdependent. As shown in [8], [9], formulations based on Reissnermixed variational theorem with independent assumptions for transverse shear stresses may be highly inaccurate when the number of layers exceeds three. The importance of satisfying the kinematic and

EP

stress continuity requirements at layer interfaces has been extensively discussed in literature and their practical importance has been proven for obtaining accurate stress predictions since Lekhnitskii [5] (1935) and Ambartsumian [6], [7] (1958) up to the recent study by Groh and Weaver [8] (2015).

AC C

The ply level stress analysis of laminated structures is carried out generally describing them as a stack-up of linear elastic, homogeneous orthotropic layers perfectly bonded together and arbitrarily oriented from each other. Inelasticity and imperfect bonding are only considered in the progressive failure analyses. The interface, i.e. the adhesive thin film bonding two adjacent layers, is usually disregarded even in the exact three-dimensional solutions ( see, e.g., Pagano [10]), but lately its effects have been shown to be relevant (see, Faraboschi [11]-[14]).

Sandwich structures are often described as three-layered plates (e.g., Faraboschi [11] and Mattei and Bardella [15]) or they are homogenized in smeared form as laminates with a thick intermediate relatively weak and deformable layer constituting the core. An accurate modelling of the cell walls behaviour of honeycomb core is only required for simulating the crushing mechanism [16]. In this case, it is mandatory to consider the transverse normal deformability.

3

ACCEPTED MANUSCRIPT Many laminate theories based on different assumptions have been developed for ply-level stress analysis of laminated and sandwich structures.Hereafter these theories are just classified on the basis of their assumptions, accuracy and costs. Their unknowns can either increase with the number of layers, or be fixed irrespectively of the number of layers, hence analysis and accuracy costs can considerably vary with the number of variables employed. As shown by the review papers on laminate theories by Noor et al. [17], [18], Carrera [19], Ghugal et al.

RI PT

[20], Reddy [21], Kreja [22] and Khandan et al. [23], a large number of research findings are related to discrete-layer models (DL) wherein unknowns are assumed layer by layer, a separate representation being adopted across each computational/physical layer. As relevant examples of DL, the papers by Reddy [1], Yaqoob Yasin and Kapuria [24] and Plagianakos and Saravanos [25], [26] are cited. DL

subdivide each constituent layer into one or an appropriate number of mathematical layers in order to

SC

capture step stress gradients within critical regions. As a consequence, their memory storage dimension and processing costs are dependent upon the number of computational layers adopted. Owing to the

M AN U

possibility they offer of improving accuracy by increasing the number of computational layers, DL are always accurate irrespectively of lay-up, loading and boundary conditions. However, due to their large number of variables, DL could overwhelm the computational capacity when structures of industrial interest made of many layers (even up to a hundred) have to be analysed.

Zig-zag plate theories based on a combination of global higher-order and local layerwise functions have been developed to achieve the accuracy of DL with a lower number of variables. Pioneered by Di Sciuva [27] (1984) and few years later with a different approach by Murakami [28] (1986), these layerwise

TE D

theories assume a coarse representation for describing the behaviour on the total laminate thickness scale, like equivalent single layer models (ESL), e.g.[29], [30], and a refined representation for describing the behaviour at the layer thickness scale. The latter contribution is constructed as the product of zig-zag amplitudes and linear zig-zag functions. In Di Sciuva’s like models the expressions of zig-zag amplitudes are derived in analytic form by the enforcement of the interfacial stress continuity conditions. In

EP

Murakami’s like models, the zig-zag function is constructed by a priori assuming a periodic change of the slope at layer interfaces, as in the exact elasticity solution of laminates with periodical stack-up [31]. Then, stress and displacement fields are assumed apart within the framework of Reissner’s mixed

AC C

variational theorem.

A comprehensive discussion of models based on Murakami’s zig-zag function for analysis of laminated structures is given by Carrera and co-workers in [19] and [32] to [36]. Comparisons among models with the same coarse representation at the overall scale and Di Sciuva’s or Murakami’s zig-zag functions are carried out by Groh and Weaver [8] and Gherlone [9]. These latter two studies show that the choice of the layer thickness scale representation is of paramount importance for laminates with general lay-up, having layers with different transverse shear and normal elastic moduli lower than those of the adjacent ones and for sandwiches with large face-to-core stiffness ratios and asymmetric properties of faces. For these cases, the zig-zag models based on Murakami’s approach are less accurate than those based on Di Sciuva’s approach with the same order of representation of displacement and stresses, while the two modelling options are equally accurate for laminates with periodic lay-up. For this reason, just Di

4

ACCEPTED MANUSCRIPT Sciuva’s like models are discussed hereafter along with discrete-layer and sublaminate zig-zag models that constitute an extension of laminate modelling by Pagano and Soni [37] merging the ideas of discretelayer and zig-zag models [38]-[40]. In sublaminate zig-zag models, the laminate is described grouping several layers within a computational layer with a zig-zag representation. The capability of zig-zag models to accurately predict displacement and stress fields has been proven for thin to thick cross-ply plates with simply-supported edges under sinusoidal distributed mechanical,

RI PT

thermal and piezo-actuating loadings in papers [41] to [47].

The reason for this choice of lay-up, loading and boundary conditions is that exact 3D solutions can be obtained , thus these benchmarks are generally assumed as test cases for comparing the accuracy of the structural models. A limited number of research findings concerns cases with more general lay-up, like angle-ply lay, asymmetric staking, distinctly different thickness and material properties of constituent

SC

layers, localized loading and other boundary conditions than simply-supported edges. Clamped edge is an interesting boundary condition by the theoretical standpoint, because the kinematic variables have to

M AN U

vanish without resulting into vanishing transverse shear resultant forces, as erroneously predicted by ESL theories having mid-plane displacements and rotations as variables.

Theories that accurately predict displacement and stress fields of laminates and sandwiches with clamped edges are developed by Tessler et al. [48], Vel and Batra [49], Swaminathan and Ragounadin [50] by adopting equivalent single-layer approaches. Mattei and Bardella [15] study cantilever and proppedcantilever sandwich beams assuming the Timoshenko kinematics across the face layers, while core is described by further considering the transverse normal deformability (Krajcinovic theory and Jourawski

TE D

post-processing). The results show that a negligible transverse shear stress in the core or a large one can rise, the behaviour being largely dependent upon the combination of boundary conditions, geometric and material properties of layers.

The reviews of models for fibre-reinforced laminated and sandwich structures [17]-[23] show that the only alternative to 3D FEA for obtaining always high-fidelity displacements, strains and stress field

EP

predictions are DL. Recently, refined zig-zag models have been developed that have shown an improved predictive capability than models reviewed in [17]-[23]. Due to the lack of extensive assessments for general lay-up, loading and boundary conditions, whether Di Sciuva’s like low cost zig-zag models are

AC C

always capable of accurately predicting displacements, strain and stress fields like discrete-layer models remains a topic still open to discussion. Even recent papers have shown discrepancies between the results of zig-zag models and exact solution for laminates with a pronounced transverse anisotropy[8], [9]. As the analyses are often carried out neglecting the transverse normal stress and strain, many authors have ascribed to this assumption the cause of discrepancies. Paper [41] demonstrated that for static problems and mechanical loading the transverse normal stress can assume a significant bearing for keeping equilibrium for laminates with a marked transverse anisotropy like soft-core sandwiches and a marked asymmetry of geometric and material properties. In this case, the transverse displacement can neither be assumed constant across the thickness, nor in a polynomial form. Instead, a layerwise piecewise variation of the transverse displacement should be assumed. Studies have shown that the same happens in static thermo-elastic

5

ACCEPTED MANUSCRIPT problems when stresses rise owing to temperature variations across the thickness. The importance of accounting for the transverse normal stress and strain in laminates under static mechanical loading is also focused by Iurlaro et al. [51].Many authors have demonstrated that a variable transverse displacement and high order zig-zag approximation of in-plane displacements is unnecessary for predicting global response quantities, in particular the first undamped free vibration frequencies (see, e.g., Iurlaro et al.[52]). However, previous considerations are correct only for the specific cases cited above, namely

RI PT

they are only valid for the cases of quasi-static response, lower vibration modes, distributed global loads and uniform plate geometries. They are not true for high-frequency vibrations and transient responses in laminates and sandwich configurations, as shown by Rekatsinas et al. [53] and they are also not true for piezo-actuating loadings, as shown by Saravanos and Heyliger [46] and Plagianakos and Saravanos [47]. Accordingly, the structural models may be required or not to describe the through-thickness variation of

SC

the transverse displacement, in order to account for the transverse normal stress and strain effects, and to use or not a higher order representation of in-plane displacements depending upon the problem

M AN U

considered.

The senior author developed a zig-zag model with a piecewise cubic representation of in-plane displacements and a constant transverse displacement for analysis of laminates with symmetric or asymmetric stack-up in [54]. A displacement-based 3D zig-zag model incorporating linear zig-zag functions in in the in-plane displacements representation and parabolic a zig-zag function in the transverse displacement was developed in [41], that satisfies the continuity of out-of-plane stresses at layer interfaces and the continuity of the transverse normal stress gradient like in the exact solutions [10].

TE D

This latter model has been refined by Icardi and Sola [43] in a hierarchical formulation that allows the representation of displacement to vary region-by-region across the thickness keeping fixed the unknowns, in order to accurately carry out the stress analysis of structures of industrial interest without needing post-processing operations.

The contribution of the present paper is to assess whether the latter zig-zag model by the authors [43],

EP

developed in the spirit of Di Sciuva’s approach as an equivalent single layer model with few variables and consequently with a low computational cost, can accurately describe displacement, strain and stress fields of laminated and sandwich plates and beams with asymmetric staking and different thickness,

AC C

transverse shear and normal moduli of constituent layers, for which relevant transverse anisotropy effects rise. The aim of the paper is to contribute to understanding whether zig-zag equivalent single-layer models have achieved the accuracy of discrete-layer and sublaminate models. Comparisons are made with recently developed multi-layered plate models for giving a contribution toward understanding the state of the art in modelling displacement and stress fields of laminates and sandwiches. It is also assessed whether accurate predictions can be obtained by the model [43] at clamped edges, other models based on the same d.o.f. being inadequate for this boundary condition as they incorrectly predict vanishing transverse shear force resultants. Understanding whether the authors’ model [43] can be always accurate like discrete-layer and sublaminate models and when its piecewise zig-zag representation of the transverse displacement can be important requires the comparison with discrete-layer and sublaminate zig-zag models for each

6

ACCEPTED MANUSCRIPT benchmark. Because no results are still available by these models for each important benchmark that present strong layerwise effects, in this paper a discrete-layer model and a sublaminate model are derived from the kinematics of model [43] . As very often in the literature, the results are presented for extremely thick laminates that seem of no practical interest for industrial applications. These cases are considered since they constitute very severe tests for the models, due to their very strong layerwise effects, and since

comparisons. The benchmarks considered in this paper are:

c)

comparison with global-local zig-zag theory by

bi-sinusoidal loading

Zhen and Wanji [45]

simply-supported rectangular sandwich plate

comparison with the zig-zag model based on

with asymmetric thickness and material

Murakami’s zig-zag function by Brischetto et

properties of faces

al.[34]

simply supported laminated [0°/-90°/0°/-90°]

comparison with the sublaminate model by Averill

beam

SC

b)

simply supported sandwich plate undergoing

M AN U

a)

RI PT

results by other researchers that are significant for assessing the accuracy of models are available for

and Yip [38] and the global-local zig-zag model by Zhen and Wanji [45]

d)

cantilever sandwich beam

comparison with RZT mixed model by Tessler et al. [48]

e)

simply-supported three-layer model

comparison with the exact solution by Faraboschi

[11]

[0°/-90/-0°] simply supported rectangular plate

g)

TE D

f)

simply supported sandwich beam with

comparison among ZZ-F model and cubic RZT model by Iurlaro et al. [51] comparison of models [41], [54] by Icardi

marked material asymmetry of faces h)

simply-supported, [90/05/90], [0/90/0/90] and

comparison with modified RZT model [8]

l)

propped-cantilever sandwich beam

AC C

i)

EP

[0/90/02/90] laminated beams

simply supported sandwich plate

comparison with model by Mattei and Bardella [15] comparison with discrete-layer model by Plagianakos and Saravanos [26]

The paper is structured as follows. A brief discussion of key features of discrete-layer, sublaminate and multilayered zig-zag models not described into details above is premised to give an overall view of available theories. Then, the constitutive equations assumed, other basic modelling assumptions and notations are summarized. After, through-thickness displacement and stress continuity requirements are discussed and the ZZ-F model [43] is explained into details. Then a discrete-layer and a sublaminate models are derived from model [43]. Finally displacements and stress fields predicted by these models

7

ACCEPTED MANUSCRIPT for the examined benchmarks and those by other authors are discussed. The displacement field of these latter models are reported in Appendix A1.

2. Considerations on discrete-layer, sublaminate and zig-zag theories Discrete-layer models can always provide the appropriate interdependent representation of displacement

RI PT

and stress fields, irrespectively of lay-up, loading and boundary conditions, by choosing an appropriate number of computational layers. They split into partial and full theories if just transverse shear stresses

are accounted for, or if also the transverse normal stress is considered in addition to membrane stresses. The kinematic and stress continuity conditions at the layer interfaces are enforced either in a point form, or in integral form. The penalty method adopted by Ascione and Fraternali [55] and Fraternali and Reddy

SC

[56] is cited as an efficient technique that enforces the constraints in a limiting sense as the penalty

parameter approaches a large value. The analysis can be carried out by discrete-layer models even using a single computational layer, using few layers for representing the whole laminate, or subdividing each

M AN U

constituent layer into one or more mathematical layers. With the first option, the order of representation and the quality of results is poor like that of equivalent single-layer models, because the stress continuity is disregarded at the interfaces. Using a highly refined representation, a high computational effort is required hence there could be the possibility to overwhelm the computational capacity when a complex industrial structure has to be analysed, like with 3D FEA( still widespread, as shown in [17]-[23] and [57]- [60]). The analysis carried out by the plate elements derived from DL is advantageous over 3D FEA because the meshing preparation is simpler and no restrictions about thickness-to-side ratios are

TE D

required, the element aspect ratio being restricted only in the two in-plane dimensions. On the contrary, solid elements should have sides with approximately a length equal to the thickness of constituent layers, so they should be very small. Hierarchical and predictor-corrector models that offer a refined modelling just in the local regions of interest and a partial modelling capability outside have been developed to limit

EP

the computational burden of DL [17],[18].

Zig-zag plate theories as equivalent single-layer theories based on a combination of a coarse global higher-order representation and local layerwise functions have been developed for achieving a

AC C

satisfactory level of accuracy with lower costs than discrete-layer models. These models are developed either on the basis of Di Sciuva’s [27] or Murakami’s [28] approaches. Di Sciuva constructs the layerwise contribution to in-plane displacements as the product of a zig-zag amplitude, which rules the magnitude of the zig-zag effect and an a priori assumed linear zig-zag function “(z-zk) Hk “ (Hk is the Heaviside unit step function: Hk=0 for z < zk; Hk=1 for z ≥zk ; z is the thickness coordinate) that makes discontinuous the slope of in-plane displacements at the layer interfaces. The expressions of zig-zag amplitudes are derived in analytic form once for all for any lay-up by enforcing the continuity of transverse shear stresses at layer interfaces. Because they depend only upon geometric and material properties of layers, the models that incorporate Di Sciuva’s zig-zag function are here termed as physically-based zig-zag models (ZZ-F). The unknown kinematic variables of the coarse representation are the only unknown and thus their number do not depend on the number of constituent layers.

8

ACCEPTED MANUSCRIPT Murakami’s zig-zag function [28] is constructed by a priori assuming a periodic change of the slope at layer interfaces, like exact elasticity solution of laminates with periodical stack-up [31]:

Fk ( z ) = ( −1) 2 k

z − zMk . For this reason, the models that incorporate Murakami’s zig-zag function are here hk

termed as kinematic-based zig-zag models (ZZ-K). Stress and displacement fields are assumed apart within the framework of Reissner’s mixed variational theorem. ZZ-K much easily allow for development

RI PT

of efficient C0 finite elements, stresses being described apart from displacements. Advantages in terms of accuracy are obtained with respect to self-equilibrated polynomial stress assumptions (Auricchio and

Sacco [62]) deriving the expressions of interlaminar stresses by integrating local differential equilibrium equations wherein in-plane stresses are obtained from kinematic relations under simplifying assumptions (Tessler et al. [48]). The accuracy of the strain energy being contemporaneously dependent upon the

SC

accuracy of stresses and strain fields, excessively simplified kinematic relations should be avoided, otherwise the response quantities might be inaccurately predicted nevertheless stress are accurate.

M AN U

Moreover, simplified kinematics could inappropriately describe the interdependence of displacement and stress fields of laminates with transverse anisotropy [8].

In the above mentioned zig-zag models, the transverse displacement is assumed constant across the thickness, only in-plane displacements being enriched with layerwise functions. In the spirit of Di Sciuva’s approach, the senior author developed in [41] (2001) a zig-zag model with a piecewise variable transverse displacement. This model adopts the Di Sciuva’s zig-zag function Φ ( x, y )( z − z k ) H k in k

each in-plane displacement to satisfy the continuity of transverse shear stresses at the interfaces and a parabolic zig-zag function Ψ ( x, y )( z − zk ) H k + Ω ( x, y )( z − z k ) H k in the transverse k

2

TE D

k

displacement that satisfies the continuity of the transverse normal stress and of its gradient across the thickness, as prescribed by the elasticity theory and as observed in the exact solutions. The continuity of the transverse normal stress and gradient directly follows from the continuity of transverse shear stresses

EP

by local equilibrium equations. As in Di Sciuva’s model [27], the expressions of zig-zag amplitudes

Φ k , Ψ k , Ω k are derived in analytic form by enforcing the continuity of transverse shear stresses, of the transverse normal stress and of its through-thickness gradient at layer interfaces.

AC C

.An interesting modelling option that merges the idea of incorporating a physically-based zig-zag function to enrich the displacement field of an equivalent single layer model, together with the idea of assuming displacements separately from out-of-plane stresses within the framework of Reissner’s mixed variational theorem is represented by the RZT zig-zag theory by Tessler et al. [48]. The displacement field is assumed as the superposition of FSDT kinematics and a zig-zag field. The characteristic feature of RZT is that the transverse shear stresses is not enforced to be continuous across the layer interfaces, the continuity being enforced only on the zig-zag dependent contribution. Relaxing the transverse shear stress continuity enforcement, the RZT model allows for a piecewise constant distribution of transverse shear stresses that is accurate in an average sense and that gives rise to a C° zig-zag theory suitable for development of efficient finite elements and for analysis of structures with clamped edges.

9

ACCEPTED MANUSCRIPT Various zig-zag theories with an enriched representation and a larger number of variables than Di Sciuva’s zig-zag models, the so-called global-local models, have been recently developed. A general global-local higher-order theory satisfying the continuity of all displacement and through-thickness stresses was presented by Shariyat [61] in 2010. The global-local pioneering theory is that by Li and Liu [63] assuming a constant transverse displacement across the thickness and a layer-by-layer representation of in-plane displacement components as the sum of a global component for each

RI PT

displacement and two local contributions for each component which satisfy the continuity of out-of-plane stresses at the interfaces of physical layers. Refinements have been brought to this theory by Zhen and Wanji [65]-[68] up to the third-order.

Sublaminate zig-zag models, which merge the ideas of discrete-layer and zig-zag models as an extension of the concept of local and global laminate modelling introduced by Pagano and Soni [37], are developed

SC

by Averill and Yip [38],, Gherlone and Di Sciuva [39] and Icardi [40]. In these models, a zig-zag

representation is assumed within each computational sublaminate to represent the through-thickness

M AN U

variation of displacements that satisfies the stress continuity conditions at layer interfaces. Analyses can be carried out either using a single computational layer, stacking several computational layers or even subdividing a physical layer into one or more computational layers. An important consequence of the possibility of refining the representation where stress gradients rise, sublaminate zigzag models can obtain accurate stress predictions from constitutive equations with a computational burden lower than DL models.

3. Theoretical formulation of the structural models of this study

TE D

Next paragraphs describe the theoretical framework, starting from material equations/ assumptions and arriving to discussion of structural models used in the simulations. Features and kinematics of the ZZ-F zig-zag model [43] are reviewed, then a discrete-layer model and a sublaminate model are derived, which are used in the numerical applications.

EP

3.1 Preliminaries and notations

For generality, a multilayered plate of overall thickness h made of nl arbitrarily oriented, orthotropic layers and ni= nl -1 interfaces is considered. The middle plane Ω is assumed as the reference plane of the

AC C

laminate and each point is referred to a Cartesian coordinate system (x, y, z), with z as the thickness coordinate (z ∈ [-h/2, h/2]). Each layer has its natural coordinates 1,2 arbitrarily oriented with respect to x,y. The position of the upper+ and lower - surfaces of the generic k-th layer are denoted by

(k) +

z and (k)z-,

respectively. All the quantities that belong to a generic layer k are denoted with the superscript(k) (see Figure 1). The strain and stress tensor components are denoted respectively as εxx, εyy, εzz, εxy, εxz and εyz and σxx, σyy, σzz, σxy, σxz and σyz. A linear elastic behaviour is assumed, inelasticity being pointless for dimensioning and assessment. Indeed, it is of interest only in the ultimate limit states, while laminates are dimensioned to satisfy the service-ability deflection that entails linear behaviour, as discussed by Faraboschi [11]-[14]. A relatively compliant connection of the layers through the interlayer gives rise to inelastic behaviour beyond the ultimate limit state, while a relatively stiff one gives rise to delamination. Such inelasticity effects are neglected, neither failure nor delamination analyses being carried out in this

10

ACCEPTED MANUSCRIPT paper. As a consequence, displacement jumps are not allowed at layer interfaces, hence each layer is assumed to be perfectly bonded to adjacent ones. Geometric nonlinearity is also disregarded, buckling problems being not considered in this paper. Just infinitesimal strains being considered, a distinction between various measures of stresses and strains is immaterial. The interlayer is described as a thin elastic constituent layer in the numerical applications when comparisons with reference solutions that accounts for its presence are carried out. Sandwich plates are

RI PT

analysed in homogenized form as laminated plates, wherein core is described as a thick intermediate elastic layer with elastic behaviour, because local deformations and stresses at the cell level of honeycomb core are outside the purpose of the analysis.

According to previous assumptions, the behaviour of each arbitrarily oriented layer is described in terms

orthotropic medium by the generalized Hooke’s law:

C13 C23 C33

0 0 0

0 0 0

0 0

C44 C45

C45 C55

C36

0

0

C16  ε xx    C26  ε yy  C36  ε zz    0  ε xz  0  ε yz    C66  ε xy 

M AN U

σ xx   C11 C12 σ    yy  C12 C22 σ zz  C13 C23  = 0 σ xz   0 σ yz   0 0    σ xy  C16 C26

SC

of averaged apparent elastic properties (macromechanics approach) as that of an homogeneous elastic

(1)

The elastic coefficients C11 to C66 are related to Young’s moduli E11, E22, E33 in the natural coordinate system(1, 2, 3) (3 and z are coinciding), shear moduli G12, G23, G31 in the planes 1-2, 2-3, 3-1,

known relations [69].

υ 21 , υ23 , υ32 , υ13 , υ31 (νij=-εj/εi) and to the orientation angle θ

TE D

Poisson’s ratios υ12 ,

by the well-

3.2 Kinematic and stress interfacial continuity requirements The kinematic and stress continuity requirements already discussed by Lekhnitskii [5] and

EP

Ambartsumian [6], [7] many years ago should be fulfilled when layers made of different materials and/or having their material axes oriented with a different angle each other.

AC C

If one assumes the layers perfectly bonded together, the displacements should be continuous across the layer interfaces:

u v w

(k ) z+

(k ) z+

(k ) z +

=u

=v

=w

(k ) z−

(k ) z−

(2)

(k ) z−

the sign – being used to indicate that the interface is approaching, while + indicate that it has been just crossed. The former relations constitute the kinematic continuity requirements (KC). The following stress continuity conditions (SC) should be contemporaneously fulfilled to satisfy Cauchy local equilibrium equations:

11

ACCEPTED MANUSCRIPT a ) σ xx , x + σ xy , y + σ xz , z = bx ; b) σ xy , x + σ yy , y + σ yz , z = b y ; c ) σ xz , x + σ yz , y + σ zz , z = bz (3) In order to satisfy equations a), b), transverse shear stresses have to be continuous at the interfaces, namely: (k )

z+

σ yz |

(k )

+

z

= σ xz | ( k ) = σ yz | ( k )

z− z

(4)

−

as a smooth in-plane variation of membrane stresses along with their derivatives

σ yy, y ,

and

σ xx, x , σ xy, x , σ xy, y

RI PT

σ xz |

is required to satisfy Cauchy equations. As a consequence of (4), it directly follows that

also the transverse normal stress and stress gradient should be continuous at the layer interfaces: (k )

σ z,z |

z+

(k )

= σ z | (k )

z

+

z−

= σ z,z | (k )

z

SC

σz |

−

in order to satisfy the third equilibrium equation c) [41].

(5)

(

) (

)

(

M AN U

The following stress-free boundary conditions hold at upper and lower bounding faces:

)

σxz z =±h2 =σyz z =±h2 =σzz,z z =±h2 = 0

(

) (

)

σzz z =±h2 = p0 z =±h2

(6)

The kinematic continuity requirements KC are satisfied apart from stress continuity requirements SC at the layer interfaces in kinematic-based zig-zag models that incorporate Murakami’s zig-zag function, in the framework of mixed variational theorems. On the contrary, displacement based zig-zag models that incorporate Di Sciuva’s zig-zag function satisfy KC and SC contemporaneously by assuming

TE D

displacement fields whose slope appropriately varies before and after the interfaces. In the next sections, the zig-zag multi-layered plate theory by the authors [43] here referred ad the adaptive zig-zag model is discussed and recast as a discrete-layer theory (variables at the layer level) and as a sublaminates zig-zag theory (variables at the sublaminate level). The predictions of these three

EP

models are compared in this paper. 3.3 Adaptive zig-zag model

The displacement field is postulated as the sum of three contributions that describe their through-the-

AC C

thickness variation:

u( x, y, z) = U 0 ( x, y, z ) + U i ( x, y, z ) + U c ( x, y, z ) v( x, y, z) = V 0 ( x, y, z ) + V i ( x, y, z ) + V c ( x, y, z )

(7)

w(x, y, z) = W 0 ( x, y, z ) +W i ( x, y, z ) +W c ( x, y, z ) The meaning of symbols and the basic steps required for deriving their expressions are explained forward. 3.3.1 Basic contribution ∆0 It brings a linear contribution in z

U 0 ( x, y, z) = u0 ( x, y) + z γ x0 ( x, y) − w0, x ( x, y) V 0 ( x, y , z ) = v 0 ( x, y ) + z γ y0 ( x, y ) − w0, y ( x, y ) 

(8)

12

ACCEPTED MANUSCRIPT W 0 (x, y, z) = w0 (x, y) which contains the five functional d.o.f. of our zig-zag theory: the displacements u0 (x, y), v0 (x, y), w0 (x, y) and the transverse shear rotations γx0 (x, y), γy0 (x, y) at the middle plane (comma indicates differentiation: [.] ,j=∂[.] ⁄∂xj). Kinematics and d.o.f. are reported in Figure 2. 3.3.2 Variable contribution ∆i

RI PT

A representation across the thickness with the desired expansion order is enabled by terms:

U i ( x, y , z ) = Ax 2 z 2 + Ax 3 z 3 + Ax 4 z 4 + ... + Axn z n

V i ( x, y, z ) = Ay 2 z 2 + Ay 3 z 3 + Ay 4 z 4 + ... + Ayn z n

(9)

W i ( x, y , z ) = Az1 z + Az 2 z 2 + Az 3 z 3 + Az 4 z 4 + ... + Azn z n

(

) (

)

(

)

σxz z =±h2 =σyz z =±h2 =σzz,z z =±h2 = 0

(

SC

Coefficients Ax2, …., Azn are computed imposing the fulfilment of stress boundary conditions:

) (

)

σzz z =±h2 = p0 z =±h2

M AN U

and satisfaction of local equilibrium equations:

σ xx , x + σ xy , y + σ xz , z = bx σ xy , x + σ yy , y + σ yz , z = by σ xz , x + σ yz , y + σ zz , z = bz

(10)

(11)

at selected points across the thickness, as described in Appendix A2. It should be remarked that no additional variables are involved, because the analytical expressions of these coefficients are obtained as function of the d.o.f. u0 , v0 , w0, γx0 , γy0 and their spatial derivatives.

A fourth-order representation for w and a third-order one for u, v is required for satisfying boundary

w = w0 for z = 0 [41]. Hence coefficients up to Ax3 , Ay 3 , Az 4 are

TE D

conditions (10) and to obtain

computed by enforcing such boundary conditions. The other higher-order coefficients are computed by enforcing the fulfilment of local equilibrium equations at a specified number of points that depends upon the expansion order that is adopted and whose position is chosen as outlined in Appendix A2.

EP

An alternative option that limits the expansion order and provides accurate results [43] is to consider

zn

just up to n = 3 for in-plane displacements and n = 4 for the transverse displacement. This means that

z n are incorporated, e.g. ( Ax3 + Aˆ x3 ) z 3 , with the former

AC C

two coefficients multiplying the same power

Ax 3 computed by enforcing the boundary conditions and the latter Aˆ x3 computed by enforcing equilibrium equations. In the numerical applications presented forward a piecewise cubic-quartic representation that is the same everywhere across the thickness is adopted, whose coefficient are recomputed for each fictitious computational layers used to subdivide the laminate thickness. Instead of using a unique representation across the thickness whose expansion order depends upon the number of boundary and equilibrium conditions that are enforced, a lower-order representation is chosen whose degree of expansion is the same elsewhere, but whose coefficients are recomputed across the thickness, this option having beenn shown to be cost-effective in [43]. Sandwiches are described using three computational layers, the lower one that incorporates the lower face and a portion of core, the intermediate layer incorporates the central portion of core across the mid-plane

13

ACCEPTED MANUSCRIPT and the upper one that incorporate the remaining part of core and the upper face. Laminates are also studied grouping the constituent layers into two or three computational layers. It is reminded that irrespectively of the number of computational layers chosen, the d.o.f. remain fixed to five. The satisfaction of local equilibrium equations is enforced at points whose position across the thickness is defined in section 3.4.4 . As shown in [43] and by the present numerical results, the fulfillment of local equilibrium equations enforced for computing the coefficients of displacements considerably improves

RI PT

the accuracy of stress predictions with respect to the model [41]with the same expansion order of

displacements across the thickness and the same zig-zag functions. A very good agreement of the

through-thickness distributions of stresses across the thickness computed from constitutive equations is obtained with respect to exact solutions, without the need for additional stress recovery steps.

Thanks to contributions (9), the theory can adapt to local variation of solutions. Indeed, the expansion

SC

order can be varied from point-to point across the thickness. For this reason it is referred throughout the paper as the adaptive model. The contributions (9) also enable satisfaction of general boundary conditions (e.g., clamped edges), because one or more coefficients Axn, Ayn, Azn can be computed by enforcing the

M AN U

fulfillment of the boundary conditions of interest. The contributions ∆c discussed in the next section 2.2.3 are aimed at making continuous displacements and stresses where coefficients are recomputed. Numerical tests have shown that this strategy preserves accuracy and considerably reduces the computational effort. 3.3.3 Zig-zag piecewise contribution ∆c

Piecewise contributions to displacement are incorporated as: ni

ni

TE D

U c ( x, y ) = ∑ Φ kx ( x, y )( z − z k ) H k + ∑ C uk ( x, y ) H k k =1

k =1

ni

ni

V c ( x, y) = ∑ Φ ky ( x, y)( z − z k ) H k + ∑ Cvk ( x, y ) H k k =1

(12)

k =1

ni

ni

ni

k =1

k =1

k =1

EP

W c ( x, y ) = ∑ Ψ k ( x, y)( z − z k ) H k + ∑ Ω k ( x, y )( z − z k ) 2 H k +∑ C wk ( x, y ) H k An appropriate kink in the slope of displacements that a priori satisfies the stress compatibility equations (4), (5) is obtained by deriving the appropriate expressions of zig-zag amplitudes Φxk, Φyk, Ψk, Ωk. Term

AC C

( z − zk ) H k is Di Sciuva’s zig-zag function, while term ( z − z k ) 2 H k is Icardi’s parabolic zig-zag function. Hk is the Heaviside unit step function (Hk=0 for z < zk, while Hk=1 for z ≥zk) and n j is the number of interfaces crossed for reaching the position z (numbered from the bottom face of the laminate). The zig-zag amplitudes Cuk, Cvk and Cwk restore the continuity of displacements (2) at the mathematical layer interfaces where the representation (7) is varied across the thickness assuming different contributions (9). They also restore the continuity of displacements at the interfaces of computational layers of discretelayer and sublaminate models derived from model (7), whose representation of displacement is assumed at the layer level or at the sublaminate level, respectively, as explained hereafter.

14

ACCEPTED MANUSCRIPT The expressions of Φxk, Φyk are obtained from the enforcement of stress compatibility equations (4), while those of Ψk, Ωk from (5), as shown in Appendix A3. The expressions of continuity amplitudes Cuk, Cvk and Cwk are given in Appendix A4. SEUPT technique [70] can be applied for obtaining a C0 formulation that overcomes the presence of unwise derivatives of the d.o.f. in the expressions of continuity functions, so to develop efficient finite elements.

RI PT

3.4 Discrete-layer and sublaminate theories

In this section, the laminated theory described in the previous paragraphs 3.3.1 to 3.3.3, here referred as the parent theory, is reformulated in order to obtain a discrete-layer theory and a sublaminate theory,

whose predictions are compared to ones of the parent theory in order to assess its accuracy, along with exact results and the results of theories by other authors taken from the literature.

SC

3.4.1 Discrete-layer theory

A discrete-layer theory is obtained assuming a separate representation of displacements within each

M AN U

computational layer. In this case, the reference plane is the mid-plane of the computational layer, while the functional d.o.f. are assumed as the in-plane components and the two shear rotations

J

J

u 0 , J v 0 and transverse

γ x0 , J γ y0

J

w 0 displacement

at the mid-plane of the layer J. Each computational

layer can either represent the kinematics across the thickness of a constituent layer, or the kinematics across a sub-layer obtained subdividing the constituent layer into sub-layers with the same material properties and orientation. In the present paper, the following displacement field is assumed across the thickness of a generic computational layer J:

] ]+ A

J

w( x, y , z ) = J w o + J A1z J z + J A2z J z 2 + J A3z J z 3 +

TE D

J

[ v( x, y, z )= v + z[ γ

J

u( x, y, z )= J u 0 + J z J γ x0 − J w,ox + J A2x J z 2 + J A3x J z 3 + 0

J

J

0 y

− J w,oy

J

y 2 J

Φ kx ( J z − J z k ) H k

z 2 + J A3y J z 3 + J Φ ky ( J z − J z k ) H k J

EP

J

J

Ψ k ( J z − J zk ) H k +

k being the position of the upper face of the computational layer J

J

Ω k ( J z − J z k ) 2 H k (13)

( J zk = J h / 2) , which represent the

AC C

interface with the next computational layer. At this interface, the stress continuity conditions (4) and (5) are enforced for obtaining the expressions of

J

k Φ kx , J Φ y , J Ψ k , J Ω k . Three variables are eliminated

u=J u ,

J +1

In (13), the explicit dependence by z is shown, all the other quantities except

Hk

at the interface by enforcing the displacement continuity conditions

The expressions of the unknown coefficients

J

A2x ,

J

A3x , J A2y ,

J +1

J

A3y , J A1z ,

v= J v ,

J +1

w= J w .

are function of (x,y). J

A2z , J A3z are

determined by enforcing: (i) the satisfaction of the stress boundary conditions (4) and (5) just at the upper or lower faces of the laminate, respectively when the computational layer J is the first or the last layer from the bottom, or (ii) the continuity of displacements across the interface with the next layer, i.e. J

u top = J +1 u bottom , J v top = J +1 v bottom , J w top = J +1 w bottom , (iii) the fulfillment of the three equilibrium

15

ACCEPTED MANUSCRIPT conditions at the point (0,0,

J

h / 2 ) if J is an intermediate layer and, in case J is an intermediate layer

(iv) only the third equilibrium equation at (0,0,0), see Figure 3. In the numerical applications, a single computational layer (13) is used to describe the behavior of each constituent layer of laminates, while three layer are used to describe sandwiches with monolithic faces (a layer for each face and two across the core). Sandwiches with laminated faces are described using a

RI PT

computational layer for each face ply and two layers for core. 3.4.2 Sublaminate theory

A sublaminate theory is obtained assuming the displacement field of the previous discrete-layer theory postulating that several layers are grouped within each computational layer (13), then in order to satisfy the stress continuity conditions at the interfaces of these layer, the zig-zag contributions by (12) are

J

[

]

SC

incorporated:

u( x, y, z )= J u 0 + J z J γ x0 − J w,ox + J A2x J z 2 + J A3x J z 3 + J Φ kx ( J z − J z k ) H k +

ni

ni

k =1

k =1

J

[

]

M AN U

∑ Φ kx ( z − zk ) H k + ∑ Cuk H k + ( Ax 4 J z 4 + ...)

v( x, y, z )= J v 0 + J z J γ y0 − J w,oy + J A2y J z 2 + J A3y J z 3 + J Φ ky ( J z − J zk )H k +

ni

ni

k =1

k =1

∑ Φ ky ( z − zk ) H k + ∑ Cvk H k + ( Ay 4 J z 4 + ...)

w( x, y , z )= J w o + J A1z J z + J A2z J z 2 + J A3z J z 3 + ni

∑ Ψ k ( z − zk ) H k + k =1

ni

TE D

J

J

Ψ k ( J z − J zk ) H k +

J

Ω k ( J z − J zk ) 2 H k +

ni

∑ Ωk ( x, y)( z − zk ) 2 H k +∑ Cwk H k + ( Az 4 J z 4 + ...) k =1

(14)

k =1

Terms incorporated in order to satisfy the stress continuity at the upper interface where the sublaminate is

EP

connected to another computational layer, or to another sublaminate, are written in separate form from those of summations (12) to underline their different meaning. The purpose of terms in the summations is making continuous interlaminar stresses at the interfaces of inner layers of the generic sublaminate J,

J

J

k Φ kx ( x, y )( J z − J z k ) H k , J Φ y ( x, y)( J z − J zk ) H k

AC C

while

,

J

Ψ k ( x, y)(J z − J zk ) H k ,

Ω k ( x, y )( J z − J z k ) 2 H k make continuous these stresses at the interfaces of adjacent layers within the

computational layer J+1, as explained for the discrete-layer model (13). Like in (13) ,all quantities except

H k are function of (x, y),

but only the explicit dependence by z is reported in (14).

The same rule described before for the discrete-layer model represented in Figure 3 holds for connecting the computational layer J to next one J+1 and the same unknown variables are involved. The position of points at which local equilibrium equations are enforced for computing

J

A2x , J A3x ,

is selected in the same way described above for the discrete-layer model. The point where local equilibrium equations are enforced for computing the unknown coefficients Ax 4 J z + ... , 4

16

ACCEPTED MANUSCRIPT Ay 4 J z 4 + ... , Az 4 J z 4 + ...

of higher-order terms have to be chosen at different locations from

previous ones. Accordingly, they are chosen uniformly spaced across the thickness and they are moved to another position when they coincide or are too close to previous ones. In the numerical applications, laminates and sandwiches are described using two computational layers (14) of equal thickness. 3.4.3 Specifications for boundary conditions

RI PT

The governing equations for the structural models of sections 3.3, 3.4.1, 3.4.2, here not reported, have been obtained in a straightforward way using the Total Potential Energy Functional. These governing

equations are solved within the framework of Rayleigh-Ritz method assuming the displacement variables expressed as series expansions of unknown amplitudes and trial functions specific for each boundary condition. for z = 0 , i.e. on the reference mid-surface:

w0 (0,0) = 0

;

w0 (0,0), x = 0

;

;

γ x0 (0,0) = 0

M AN U

u 0 (0,0) = 0

SC

Cantilever beam. The following essential boundary conditions are enforced at the clamped edge x = 0

(15)

In order to simulate identically vanishing d.o.f. across the thickness at the clamped edge, the following further boundary conditions are enforced:

u 0 (0, z), z = 0

;

w0 (0, z), z = 0

w0 (0, z), xz = 0

;

(16)

The natural boundary condition is also enforced at the clamped edge:

∫

h/2

−h / 2

σ xz (0, z )dz = T

(17)

edge x=L

TE D

in order to impose the transverse shear stress resultant force equals the constraint force, while at the free

∫

h/2

−h / 2

σ xz ( L, z ) dz = 0

(18)

representation:

EP

is enforced . An approximate solution is obtained assuming for each displacement d.o.f. the following

x d = ∑∆   L i =1

AC C

8

i d

i

( d ≡ u 0 , w 0 , γ x0 )

(19)

Boundary conditions (16) to (18) are enforced using Lagrange multipliers method. The expressions of nine amplitudes ∆ d , three for each displacement, are obtained from the enforcement i

of boundary conditions (15) to (18), the remaining fifteen amplitudes, five for each displacement, are determined by Rayleigh-Ritz method. Propped-cantilever beam. The same boundary conditions (15) to (18) hold, but the additional one adds at x=L, being assumed that the support condition is realized at

h z = − , i.e. on the lower face opposite to 2

that the load is applied, in accordance with the reference solution [15]. In this case, the expansion order is increased to 9:

17

ACCEPTED MANUSCRIPT 9 x d = ∑ ∆id   L i =1

i

( d ≡ u 0 , w 0 , γ x0 )

(20)

Simply-supported edges. The following essential boundary conditions are enforced at x = 0 ,

y = 0 , y = Ly

w0 (0, y) = 0 ; w 0 ( Lx , y ) = 0 ; w0 (0, y), xx = 0 ; w0 ( Lx , y), xx = 0 w0 ( x,0) = 0 ; w0 (0, Ly ) = 0 ; w0 ( x,0), yy = 0 ; w0 ( x, Ly ), yy = 0

RI PT

and

x = Lx

(21)

according to [10]. The solution for this case is obtained in the following form: M N   nπ   mπ x  sin  y ; v 0 ( x, y ) = ∑∑ Bmn sin  L  m =1 n =1   y   Lx

M N  mπ w 0 ( x, y ) = ∑∑ C mn sin  m =1 n =1  Lx

  nπ  x  sin  y ;     Ly 

N

 mπ  Lx

m =1 n =1

M AN U

M

γ x 0 ( x, y ) = ∑∑ Dmn cos

  nπ  x  cos y ;     Ly 

SC

M N  mπ u 0 ( x, y ) = ∑∑ Amn cos m =1 n =1  Lx

M N   nπ   mπ 0 x  sin  y ; γ y ( x, y ) = ∑∑ E mn sin  L  m =1 n =1   y   Lx

For cylindrical bending, the solution reads:

(22)

  nπ  x  cos y ;     Ly 

TE D

M M M  mπx   mπx   mπx   ; w0 ( x, y ) = ∑ Cm sin  ; γ x0 ( x, y ) = ∑ Dm cos  u 0 ( x, y ) = ∑ Am cos m =1 m =1 m =1  Lx   Lx   Lx 

(23)

having assumed that bending occurs is in the plane (x,z), hence the variable y is neglected.

 mπx 

 nπy 

0  or sinusoidal  sin  If a bi-sinusoidal P = P sin   L  L  x   y 

 mπx   load is applied P = P 0 sin  L  x 

EP

(usually on the upper face in the literature), just the first component m=1, n=1 is required in (22), (23) for solving governing equations.

AC C

3.4.4 Number of computational layers and equilibrium points adopted The computational layers and the positions of points across the thickness where equilibrium equations are enforced for deriving the expressions of higher-order coefficients are chosen as follows. Sandwiches. When the analysis is carried out by the adaptive model [43], three equally spaced fictitious computational layers across the thickness are considered, each with a piecewise third-order representation of in-plane displacements and a fourth-order one for the transverse displacement. A total of np=Nlay · ord_u – 2 points are required for computing the higher-order coefficients, Nlay being the number of computational layers and ord_u the order of the expansion of the in-plane displacement (i.e. ord_u=3), hence seven points should be selected across the thickness. In all benchmarks involving sandwich plates and beams their positions are chosen at: [-0.45 h; -0.42 h; -0.3h; 0.1h; 0.3h ;0.42h; 0.45h]. The analysis with the discrete-layer model of section 3.4.1 is carried out subdividing the sandwich plate or the sandwich beam into four computational layers, two corresponding to the face layers and two

18

ACCEPTED MANUSCRIPT obtained subdividing the core into two computational layers of equal thickness, for a total of 20 unknowns.. The points of Figure 3 are adopted for enforcing equilibrium equations and stress boundary conditions. The analysis with the sublaminate model of section 3.4.2 is carried out representing the sandwich plate or beam with two computational layers of equal thickness, for a total of 10 unknowns. Coefficients

J

A2x , J A3x , J A2y

J

A3y , J A2z , J A3z are computed by enforcing the satisfaction of

The coefficients

Ax 4 , Ay 4 , Az 4

RI PT

equilibrium equations as discussed in section 3.4.2 and shown in Figure 3. These are termed as points P1. of fourth-order powers of the thickness coordinate in Eq. (14) are

recomputed across the thickness of each computational layer also enforcing equilibrium equations at

selected points. These points, here termed P2, are assumed at the same position as those of the adaptive model. Whenever the selected point P1 and P2 are too close each other, either the point P1 or the point

SC

P2 are retained, as both at the same position cannot improve accuracy. The present choice of the number of computational layers will result into a degree of accuracy of discrete-layer and sublaminate models equal or lower than the adaptive model. Of course, increasing the number of computational layers and/or

M AN U

the expansion order a better accuracy could be obtained by sublaminate and discrete-layer models Laminates. The analysis of laminates is carried out by adaptive and the discrete-layer models representing each constituent layer with a computational layer, while the whole laminate is subdivided into two computational layers when the analysis is carried out by the sublaminate model. Due to this choice and the related number of points across the thickness used for computing higher-order coefficients, the accuracy of the sublaminate model is lower or equal to that of the adaptive model, while that a of the

adaptive model.

TE D

discrete-layer model is better or equal. In this way upper and lower error bounds are established for the

The following subdivision into computational layers and positions of equilibrium points are adopted for the specific cases considered forward in the numerical assessments. The [0°/-90°/0°/-90°] laminate is studied with the adaptive and discrete-layer models using four computational layers coinciding with the

EP

constituent layers, like any other laminate with four layers. The points [-0.45h; -0.4h; -0.2h; -0.15h; 0.1h; 0.05h; 0.1h; 0.2h; 0.27h; 0.45h] uniformly distributed across the thickness are selected for the adaptive model. This uniform subdivision can be adopted whenever the layers are made of the same

AC C

material, otherwise the positions of points should be chosen where the largest errors are shown in the fulfilment of equilibrium equations. It is reminded that irrespectively of the number of computational layers considered, the adaptive model has always 5 d.o.f, while the d.o.f of discrete-layer and sublaminate models have 5xN d.o.f., N being the number of layers. The analysis with the sublaminate model is carried out computing two sets of coefficients

Ax 4 , Ay 4 , Az 4 , two computational layers being

adopted, selecting the position of the points at which equilibrium is enforced as indicated in Figure 3 and the additional points at the same positions of those of the adaptive model The [0°/-90°/0°] laminate considered forward is subdivided into 3 computational layers (like any threelayer laminate with the layers of equal thickness), when the analysis is carried out by adaptive and discrete-layer models. The points positioned at [-0.425h; -0.3422h; -0.16h; 0h; 0.15h; 0.3468h; 0.425h]

19

ACCEPTED MANUSCRIPT are selected for the adaptive model. The analysis is carried out by the sublaminate model choosing the points as in Figure 3 and the remaining points at the same positions of the former ones. This strategy is followed also for the other laminates considered in the benchmarks. The equilibrium points of the adaptive model for the laminate [90/05/90] are selected at [-0.475h; 0.425h; -0.365h; -0.345h; -0.315h; -0.25h; -0.2h; -0.15h; -0.08h; 0.02h; 0.08h; 0.15h; 0.2h; 0.25h; 0.315h; 0.345h; 0365h; 0.425h; 0.475h], those for laminate [0/90/0/90] at [-0.45h; -0.4h; -0.2h; -0.15h; -0.1h;

RI PT

0.05h; 0.1h; 0.2h; 0.27h; 0.45h], and for laminate [0/90/02/90] at [-0.45h;-0.25h; -0.15h;-0.1h;--0.07h; 0.05h; 0.1h; 0.125h; 0.19h; 0.25h; 0.36h; 0.42h; 0.47h]).

For these laminates each computational layer of the sublaminate model is expanded again up to the fourth-order and coefficients

Ax 4 , Ay 4 , Az 4

of fourth-order powers are recomputed across the

SC

thickness for obtaining a number of equilibrium points similar to that of the adaptive model.

4. Simply-supported sandwich plate: benchmark S1

M AN U

As first test, a three-layered plate with sides parallel to x-axis of length L and sides parallel to y-axis of length B, made of two external isotropic layers each one of thickness h and an isotropic interlayer of thickness 2 λ , restrained at the boundary and subjected to a single-component lateral bi-sinusoidal load P. The boundary conditions (21) hold for this case and the solution is searched as in (22) assuming m,n,=1. The material properties of the layers are the elastic modulus Poisson’s ratio

E , the shear modulus G and

υ . The elasticity solution for this sample case was obtained by Faraboschi [11] under

assumption of anti-symmetric behavior with respect to the middle plane (middle plane is immovable),

TE D

linear elastic behavior of constituent layers and fulfillment of Kirchhoff-Love hypotheses in each layer and continuity of displacements at the interfaces (2). The transverse shear forces in the interlayer and the bending moments across its cross-section are assumed negligible. To solve this benchmark problem, the structural model of section 3.3 is recast by assuming the

simplified form

EP

displacement field made of the basic contribution (8), the higher order contributions (9) in the following

U i ( x, y , z ) = Ax 2 z 2 + Ax 3 z 3

AC C

V i ( x, y, z ) = Ay 2 z 2 + Ay 3 z 3

(24)

Wi = 0

and also zig-zag contributions (12) is assumed in simplified form as: ni

ni

k =1

k =1

ni

ni

k =1

k =1

U c ( x, y ) = ∑ Φ kx ( x, y )( z − z k ) H k + ∑ C uk ( x, y ) H k

V c ( x, y) = ∑ Φ ky ( x, y)( z − z k ) H k + ∑ Cvk ( x, y ) H k

(25)

W c ( x, y) = 0

20

ACCEPTED MANUSCRIPT The expressions of zig-zag amplitudes Φ x , k

Φ ky

are derived from the enforcement of transverse shear

stress continuity conditions (4) at the layer interfaces. The expressions of coefficient

Ax 2 , Ax3 , Ay 2 ,

Ay 3 are obtained by enforcing the fulfilment of stress boundary conditions at top and bottom faces of the three-layered plate:

σ xz ( z = ± h 2 ) = σ yz ( z = ± h 2 ) = 0 and satisfaction of local equilibrium equations:

σ xx, x+σ xy , y + σ xz , z = bx ; σ xy , x + σ yy , y + σ yx, z = by

RI PT

(25)

(26)

at the top and bottom interfaces of the interlayer. The computational layers used with adaptive, discrete-

SC

layer and sublaminate models of sections 3.3. and 3.4 and the points across the thickness where the satisfaction of local equilibrium equations are enforced are described in the former section 3.4.4.

The present adaptive zig-zag, discrete-layer and sublaminate models cannot provide a linear zig-zag

M AN U

representation as in the reference solution [11] and do not neglect bending and transverse shear stresses across the interlayer. However they give a piecewise-linear-like zig-zag representation of the displacement field for the relatively thin plates with thin interfaces considered in the numerical applications. In this case, bending and transverse shear stresses become negligible across the interlayer, like in the exact solution [11]. The solution is obtained under a double sinusoidal distribution of the transverse load by assuming the trigonometric distribution of the variables (22).

TE D

Various sample cases, the same considered by Faraboschi [11], are studied whose elastic and geometric properties are reported in Table 1, along with the results. The maximum deflection the maximum stress in the plate

wMax across the span,

σ Max , and the maximum values tum , tvm of the shear stress transferred

through the upper interface between the interlayer and the upper layer predicted by the present adaptive

EP

zig-zag, discrete-layer and sublaminate models are compared with those by Faraboschi. The symbols ADZZ, ADL and ASB in Table 1denote the results by the present adaptive zig-zag, discrete-layer and sublaminate models of sections 3.3., 3.4.1 and 3.4.2, respectively. These results are

AC C

shown in a good agreement each other and with the elasticity solution [11], irrespectively of the length-tothickness ratios, the elastic properties and the dimensions considered. This proves the capability of the adaptive zig-zag model to provide results as accurate as discrete-layer and sublaminate models with a lower number of variables, in a well agreement with the exact solution for this benchmark. In the next sections other benchmark are considered to assess whether this capability is always shown and to assess whether and when a higher-order piecewise representation of the transverse displacement provides a superior accuracy.

21

ACCEPTED MANUSCRIPT 5. Cantilever and propped-cantilever sandwich beams: benchmarks S2, S3 In this section, a sandwich beam with a uniformly distributed transverse load acting on the upper face, which is clamped at x=0 and free at the other edge, case S2, and a sandwich beam that is clamped at one edge and restrained at other one at the lower face, case S3, are considered. 5.1 Cantilever sandwich beam S2

RI PT

Now the cantilever sandwich beam formerly considered by Tessler et al. [48] is analysed by the models of sections 3.3 and 3.4. The sandwich beam is subjected to a uniform transverse loading acting on the upper face is considered. Since for this case no exact solution can be obtained, the comparison is made with the 3D finite element results of [48] and the results by the RZT model [48]. In this model, a piecewise linear representation of in-plane displacements, a constant transverse displacement and a transverse shear stress

framework of Reissner’s mixed variational theorem.

SC

obtained from the membrane one by integrating local equilibrium equations are adopted within the

The faces are made of unidirectional Carbon-Epoxy laminates, while the core is made of PVC foam. The

M AN U

material properties for this case are reported in Table 2.

The length-to-thickness ratio considered for this case is 10, as in [48] and the thickness ratios of the constituent layers are [0.1h/0.8h/0.1h]. The results for this case are reported in normalized form as:

σ xz =

 Lx L y  L L  h h2 σ , , z σ = σ xx  x , y , z  xx  xz  0 2 0 Lx p Lx p  5 2   5 2 

(27)

Figure 4 shows that results by the structural models of section 3.3. and 3.4 are in good agreement with those of the RZT model [48], and with 3D FEA [48]. This can be viewed as a confirm of how powerful is

TE D

the mixed approach adopted by the RZT model in improving accuracy of stress predictions, thus it can be worthy to develop a mixed version of the adaptive model. On the other hand, since the results by the present higher-order models agree with those obtained by RZT, it results that for this benchmark just a piecewise linear representation of in-plane displacements and a constant transverse displacement suffice.

EP

The results presented next will show that it is not true in general that a higher-order zig-zag model with a piecewise representation of the transverse displacement has no practical meaning. As shown for the benchmark S3 [15], the transverse normal deformability of soft-core sandwiches is important to obtain

AC C

accurate stress predictions near the supported edge of the propped-cantilever sandwich beam. Also benchmarks S6 and T2, will demonstrate that a superior performance is obtained with a refined representation of the transverse normal deformability. The results of Figure 4 show the suitability of the adaptive model for capturing the transverse shear stress distribution across the thickness near the clamped edge of a cantilever sandwich beams, with no additional variables required. Indeed, all boundary constraint relations like those of 3.4.3 are enforced in the symbolic calculus phase carried out once for all, obtaining closed-form relations that involve only the chosen functional d.o.f. and material/geometric properties of the structure. The same can be done for arbitrary boundary conditions, like at the ends of the overlap region of bonded joints as shown in [71]. The computational efficiency of the adaptive model is demonstrated by the comparison of the processing time required to compute stress and displacement fields, which is of 1.721 s on a laptop computer, and

22

ACCEPTED MANUSCRIPT that of discrete-layer model which is more than twice. However, a symbolic calculus run that takes up to 600 s have to be carried out once for all for obtaining the closed form relations that define all quantities of the adaptive model. 5.2 Propped-cantilever sandwich beam S3 The sandwich beam with faces made of the same material but having a different thickness, formerly studied by Mattei and Bardella [15] is solved using the model of sections 3.3. with the number of

RI PT

unknowns unrelated to the number of layers, and the discrete-layer and sublaminate models of 3.4. The lower face layer has a thickness tl = c/2, while the upper face layer has a thickness tu=c/4, c being the thickness of core (h=tl+c+tc). All constituent layers are isotropic, with the following ratios of elastic moduli Eu/El=1.6, Eu/Ec=1.66.6 and a Poisson’s ratio

υ = 0.3 . A length-to-thickness ratio Lx/h=5.714 is

considered and a uniform transverse load is applied on the upper face.

SC

The beam is clamped at x=0 and restrained on the lower face at x=L. The boundary conditions for the propped-cantilever beam and the in-plane representation (20) of displacement hold for this case. As

M AN U

shown by Mattei and Bardella [15], the core compliance has a larger influence in amplifying the contribution of shear deformation to the deflection than in the cantilever beam. The exact solution for this benchmark being not available, the results by Mattei and Bardella are compared to FEA. However the finite element model predicts a different reaction force on the support than the zig-zag model adopted in [15] of about 4%. This results in a larger relative error in the central region of the beam, where the shear force changes sign. The results of Figure 5a) represents the transverse shear stress field across the thickness of the beam at

x / Lx = 7 / 8 , i.e. near the “severe boundary condition” at the restrained edge

edge

TE D

as termed by Mattei and Bardella , while Figure 5 b) represents the same stress field at the restrained

x / Lx = 1. The results are reported in normalized form as

σ xz =

Aσ xz ( x, z ) P 0 Lx

(28)

EP

A = Al + Au + Ac being the cross section of the beam. The points across the thickness where the satisfaction of equilibrium equations is enforced are defined in section 3.4.4. The results by the present

AC C

adaptive, discrete-layer and sublaminate models being in good agreement each other and with the reference solution by Mattei and Bardella confirm that the discrepancies with FEA shown in [15] should be ascribed to a physiological error of the finite element model. As shown by the inset of Figure 5a, the present adaptive model of section 3.3 provides results in a well agreement with the reference solution [15] also at the clamped edge

x / Lx = 0 . Of primary importance by the viewpoint of modeling the stress

fields at the supported edge is the solution [15] for S=0 reported in Figure 5 a) and b) that has been obtained in obtained [15] neglecting the transverse normal deformation. The comparison with the other solutions reported in Figure 5 a) and b) show that this contribution is essential for predicting the stress fields near the restrained edge of this benchmark. Thus cases/locations exist where an accurate modeling of the transverse deformation improves accuracy.

23

ACCEPTED MANUSCRIPT The deflection along the span of the beam is not reported to contain the number of figures, however the results obtained with the present models are also in a good agreement with the reference solution when computed as the average displacement across the thickness at any x . It was chosen to omit these results because they are less critical to capture than stresses, as shown in [15] the FEA results being in a much better agreement with the model by Mattei and Bardella than for stresses.

RI PT

6. Simply-supported sandwich plates: benchmarks S4 to S6 6.1 Sandwich plate S4

The simply supported sandwich plate undergoing a sinusoidal loading formerly studied by Zhen and

Wanji [45] is analysed, in order to compare the predictions of their global-local zig-zag model and the

present adaptive, discrete-layer and sublaminate models with a piecewise zig-zag representation of the

SC

transverse displacement. A length to thickness ratio of 4 is considered. Faces are 0.1h thick (MAT p) and core is 0.8h thick (MAT c) of Table 2. Figure 6 shows the transverse shear stresses in normalized form:

M AN U

σ xz

 L  L  hσ xz  0, y , z  hσ yz  x ,0, z   2  σ =  2  = yz 0 Lx p Lx p0

(29)

predicted by the global–local approach of Ref. [45], the exact 3D solution, the adaptive zig-zag model of section 3.3. and by the discrete-layer and sublaminate models of section 3.4. The exact solution shows that for this case the load applied on the upper face results into asymmetric stress fields across the thickness, thus a full 3D representation is important.

TE D

The analysis is carried out by the present adaptive model using 5 computational layers, 2 for each face and 1 for the core. The representation (22) is adopted for the in-plane variation of the functional d.o.f, assuming m=1, which is the exact solution under bi-sinusoidal loading. The fulfilment of equilibrium equations is imposed at the selected points indicated in section 3.4.4 for the adaptive, discrete-layer and sublaminate structural models considered in this paper.

EP

The results of Figure 6 show that while the global-local approach [45] only gives a hint of the asymmetry in σ̅yz due to the assumed constant transverse displacement, the present adaptive model is able to precisely describe the stress field as compared to the exact solution. This capability is obtained by the

AC C

present equivalent single layer, adaptive zig-zag models thanks assumption of a piecewise variable transverse displacement, as well as to the possibility of computing higher order coefficients by enforcing the fulfilment of local equilibrium equations. As a result, accurate stress predictions are obtained from constitutive equations, i.e. no post-processing operations are required. The present adaptive zig-zag model computes the solution in 1.752 s while its sublaminate formulation takes 9.01 s. So the sublaminate formulation, even if accurate, claims for extremely higher computational times. The same can be said for the discrete-layer model, a little larger processing time of 10.3 s than the sublaminate model being required for carrying out the analyses on a laptop computer with dual-core CPU 2.20 GHz, 64 bit operating system and 4 GB RAM. To understand what such computational times mean, the FSDT theory [29] claims for computational times ranging from 1.285 s to 1.458 s to analyse the all

24

ACCEPTED MANUSCRIPT the benchmarks considered in this paper. Similar times are required by HSDT [30]. However, the results by FSDT and HSDT theories are not reported in this paper because stresses are incorrectly predicted.

6.2 Sandwich plate S5 The simply supported rectangular sandwich plate under bi-sinusoidal loading analysed by Brischetto et al. [34] is now considered in order to compare the results by the present zig-zag model incorporating zig-zag

RI PT

functions whose amplitudes are computed by enforcing the continuity of stresses at the layer interfaces to those of a higher-order zig-zag model with variable transverse displacement in mixed form

incorporating Murakami’s zig-zag function. Note that the layerwise theory [34], which considers various sub-theories with a different order of representation, describes the through-the-thickness variation of

are reported as comparison results in Figure 7.

SC

displacements using up to a seventh order representation. The results by the most accurate sub-theory [34]

The sandwich plate has a length to thickness ratio of Lx/h=4 and a length-side ratio of Ly/Lx=3. A

M AN U

geometrical asymmetry due to the difference between the thickness of the upper and lower faces characterizes this case, the thickness ratios of the three constituent layers being [2h/10; 7h/10; h/10]. There is also an asymmetry in the material properties, the lay-up being [MAT bl/ MAT bc/ MAT p] (properties in Table 2). Figure 7 reports the comparison between the exact solution and the numerical results by the present adaptive zig-zag, discrete-layer and sublaminate models (equilibrium points as described in 3.4.4) . The quantities reported in Figure 7 are normalized as follows:

 ,zh 2   q 0 Lx

σ xz =

Ly

 L  u  0, y , z  E1MATbc h 2 2  u=  q 0 L3x

(30)

TE D



σ xz  0,

The computational time by the present adaptive model is 1.789 s, while the discrete-layer and the sublaminate models are up to twice. The results of Figure 7 show that all structural models give accurate stress predictions. However Figure 7b shows the limit of using Murakami zig-zag function for a structure

EP

with asymmetric lay-up. As a matter of fact, the physically-based zig-zag theory [43] is more accurate in describing the through-the-thickness distribution of the in-plane displacement and related in-plane stress and strain. This represents the advantage that can be achieved using Di Sciuva’s zig-zag function. The

AC C

theory [34] accurately predicts stresses nevertheless displacements are not predicted very accurately, as stresses are computed separately from displacement within the framework of the Reissner’s mixed variational theorem.

6.3 Sandwich beam S6

A case with strongly asymmetric stress fields, previously considered in [41], having a strong transverse anisotropy is now analysed. The present benchmark is that of a simply supported sandwich beam with material asymmetry of faces and a length- to-thickness ratio of 4 under sinusoidal loading. It is simulated as a multilayered beam with a (MAT 1/2/3/1/3/4)s stacking sequence and the following thickness ratios of the constituent layers (0.010/0.025 /0.015/0.020/0.030/0.4)s., whose properties are reported in Table 2. The results for this case predicted by the adaptive, discrete-layer and sublaminate models are compared

25

ACCEPTED MANUSCRIPT to the exact 3D solution [41] in Figure 8 and to those predicted by the PWM 3 model [41] with a constant transverse displacement and a piecewise cubic zig-zag representation of the in-plane displacement. For this benchmark, the core is subdivided into three computational layers of uniform thickness and each face into two computational layers. Initially the properties of faces are those of the health material reported in Table 2, then the elastic modulus E3 of the core and of the upper face layers are reduced of a factor 10-2 , as the residual properties

RI PT

after failure accordingly with the ply-discount theory This enhances the importance of the transverse

normal deformability. As shown by the results reported in Figure 8, which are normalized as follows:

σ xz ( 0, z ) p

0

σ zz =

u=

E2u ( 0, z ) hp 0

L  E2 w  x , z   2  w= hp 0

(31)

SC

σ xz =

 Lx  ,z  2  p0

σ zz 

assumption of very soft properties in the thickness direction under an applied transverse load implies that

M AN U

the transverse displacement and the transverse normal stress and strain assumes a primary role, differently to the initial case (UNDAM). The numerical results show that the different thickness of face sheets and their rather different elastic properties result into strong layerwise effects, owing to a strong transverse anisotropy, so undamaged and damaged cases exhibit totally different trends. When the material properties are asymmetric, i.e. the damaged case, also displacement and stress fields are strongly asymmetric. It is interesting to notice the importance of the modulus E3, that though considerably smaller than E1 has a significant bearing in determining stress and displacement fields for the damaged case. It is

TE D

underlined that the theory PWM3 with a constant transverse displacement is unable to capture the central role played by the transverse displacements stress and strain of this case . These results confirm that it is mandatory having a variable representation of the transverse displacement when soft materials are interfaced to other materials. The comparison with the exact solution shows that even when very strong layerwise effects rise, zig-zag theories with a high-order representation of displacements can be

EP

successfully used for saving cost.

AC C

7. Benchmarks T1 to T6 7.1 Laminate T1

The simply supported laminated [0°/-90°/0°/-90°] beam undergoing sinusoidal loading formerly analysed by Zhen and Wanji [45] and by Averill and Yip [38] is now studied. Its four layers are each h/4 thick and made of MAT p (properties in Table 2). Also for this case, a length to thickness ratio of 4 is assumed. Figure 4 reports the results by the exact 3D elasticity solution and the ones by the present adaptive, discrete-layer and sublaminate models, by the global-local zig-zag model by Zhen and Wanji [45] and by the sublaminate zig-zag model by Averill and Yip [38]. Displacement and stresses are normalized as:

σ xx =

 Lx  ,z  2  p0

σ xx 

σ xz =

σ xz ( 0, z ) p

0

u=

E2u ( 0, z ) hp 0

(32)

26

ACCEPTED MANUSCRIPT The results for this benchmark re reported in Figure 9, which show that all theoretical models of this paper are everywhere in good agreement with the exact solution. Unexpectedly, the results by the sublaminate theory [38] are less accurate in describing the through the thickness distribution of the transverse shear stress than those by the other theories. For what concerns processing time costs, in this case the adaptive model requires 1.66 s to perform the analysis while a much larger time of 8.62 s is required if cast in sublaminate form and a comparable processing time for the present discrete-later

The effects of an extreme variation of the elastic modulus E3

RI PT

model.

is also reported in Figure 9, the case with

2

E3 of the upper layer reduced of a 10 factor being presented as the damaged case, which is indicated as the damaged case. The results show that this variation of E3 has non-significant effects on stresses and displacements for a laminate with the layer having the same properties and a different orientation, while it

SC

was shown important for a soft-core sandwich in a previous benchmark. In the present case variations lower than 10% are shown with respect to the undamaged case. The importance of considering the

M AN U

transverse normal deformability will become even smaller for laminates with a large length-to-thickness ratio. However, consideration of the transverse normal deformability in the analysis of these cases by the present adaptive model does not result into an increased cost due to a longer processing time because the number of unknown variables is the same of lower-order models and the high-order contributions, which currently have a marginal importance, are fast computed from closed-form expressions obtained once for all by symbolic calculus.

The assumption of a piecewise variable transverse displacement is important for sandwiches, as shown

TE D

for benchmark S6, because the elastic properties of core are extremely lower than those of faces and the core is a very thick layer, thus a strong transverse anisotropy effect rises. 7.2 Laminate T2

Here a [0°/-90°/0°] simply supported rectangular (Ly=3Lx) plate, with a length-to-thickness ratio of 5, subjected to bi-sinusoidal transverse loading is considered. All three constituent layers are made of MAT

EP

g (properties in Table 2) The outer ones have thickness 0.25h, the inner one has thickness 0.5h. The results by the RZT cubic theory by Iurlaro et al. [51] are compared to the exact 3D solution and the results by the present adaptive, discrete-layer and sublaminate models. For this benchmark, are also reported the

AC C

results obtained recasting the adaptive model in mixed form, by assuming the transverse shear and normal stresses obtained from in-plane stresses, the latter being computed by the kinematics of the model, and integrating local equilibrium equations. Two mixed models are derived, one with a piecewise linear variation of displacements (LIN-ADZZ) and one with a piecewise cubic representation (CUB-ADZZ). Stresses and displacements are reported in Figure 10 in the following normalized form according to the reference paper [51]:

27

ACCEPTED MANUSCRIPT

{σ

}

h L p0

xz

; σ yz =

xx

; σ yy ; σ zz =

{u; w} = 10

4

2 x

  Ly   Lx   σ xz  0, , z  ; σ yz  , 0, z    2     2 

}

L L  L L   h 2   Lx Ly  σ , , z  ; σ yy  x , y , z  ; σ zz  x , y , z   2 0  xx  Lx p   2 2   2 2   2 2  

D11 L4x p 0

  Ly   Lx Ly   u  0, , z  ; w  , , z     2   2 2  

(33)

RI PT

{σ

D11 being the bending stiffness. Nevertheless the material properties of constituent layers and the lay-up are not dissimilar to those of previous cases, the results of Figure 10 show that a piecewise linear

SC

representation with a constant w cannot be accurate for this benchmark. Indeed, the results by the present adaptive, discrete-layer and sublaminate models, as well as those of the CUB-ADZZ model are in a well agreement with the exact solution, while those by the LIN-ADZZ model are not. It is also shown that the

M AN U

results by the present adaptive model are much accurate of those by the RZT cubic model [51] that assumes a piecewise cubic/quadratic displacement field for in-plane/transverse components respectively, while transverse shear and normal stresses are obtained integrating local equilibrium equations within the framework of the Reissner’s mixed variational theorem.

A high order piecewise representation of displacements is required in this case, because otherwise little errors in all displacement and stress fields are made, though a correct trend is shown across the thickness. . In terms of computational effort, for this case, the theory [43] carries out the analyses in 1.678 s, the

TE D

LIN-ADZZ theory takes 1.75 s, the CUB-ADZZ theory takes 7.55 s. A larger processing time is required by the LIN-ADZZ theory to carry out post-processing operations (integration of local equilibrium equations) in order to obtain accurate out-of-plane stress predictions.

EP

7.3 Sandwich T3

In this section, the predictions by the present discrete-layer model are compared to those by adaptive and sublaminate models, to the exact solution and to the results by the discrete-layer model by Plagianakos

AC C

and Saravanos [26] for a simply supported sandwich plate undergoing a bi-sinusoidal loading. The sandwich plate is simulated as a four-layer multilayered plate with graphite-epoxy faces and foam core (material properties in Table 2 as Gr-Ep and foam respectively) with the following thickness ratios of the constituent layers (0.1h/0.4h)s. The plate is square with sides of length 1 m and a length-to-thickness ratio of 10. For this case, the stress and displacement fields across the thickness are reported in Figure 11. The variation of the transverse shear stress transverse displacement

σ xz across the thickness is shown in Figure 11a, while the

u is shown in Figure 11b. According to [26], the transverse shear stress and the

in-plane displacement across the thickness under load amplitude 1kPa are reported without any normalization. Five computational layers have been used for solving the problem with the present discrete-layer theory, the core being subdivided into three layers of equal thickness and two computational layers being used to represent the faces. This initially results into 25 unknown variables,

28

ACCEPTED MANUSCRIPT which reduces to 13 after enforcement of the continuity of displacements at the interfaces (the continuity of stresses is achieved by computing the zig-zag amplitudes

J

k Φ kx , J Φ y , J Ψ k , J Ω k of Eq. (13) as

indicated in sections 3.4.1), since at each of the four interfaces three unknowns variables are eliminated. The same number of variables is used by Plagianakos and Saravanos [26]. The results of Figure 11 show that all the theories compared in this paper obtain results in a good agreement with the exact solution and

RI PT

each other. This result shows again that the present adaptive model with just 5 unknown variables is adequate for displacement and stress analysis. However, the comparison of the processing time by adaptive model in its original form [43] with that of the discrete-layer theory obtained from it as

described in the section 3.4.1 reveals the advantage of the former: it requires just 1.778 s to perform the analysis, while the latter requires 113 s.

SC

7.4 Laminates T4 to T6

Simply-supported laminated beams with a length-to-thickness ratio of 8 and subject under sinusoidal loading, with layers made of different materials and with a different thickness are analyzed. Current

M AN U

benchmarks T4 to T6 correspond respectively to laminates G, J and L with a length-to-thickness of 8 and a strong transverse anisotropy considered in [8], [9]. T4 is a [90°/0°5/90°] beam with the following thickness ratio of layers [0.12/0.23/0.12] made of materials [p2/pvc/h/pvc/p2] whose properties are reported in Table 2, that is made of a combination of three distinct materials whose 90° outer are weak layers. T5 is an anti-symmetrically [0°/90°/0°/90°] beam with layers of equal thickness [0.244] made of the same material properties, whose layerwise effects are generated by the lay-up asymmetry. T6 is a [0°/90°/0°2/90°] beam whose layers have the following ratios [0.3/0.2/0.15/0.25/0.1] and are made of

TE D

materials [p3/m/p]. Material p is a carbon-fibre reinforced plastic, m is a material with increased transverse stiffness, pvc is a closed-cell polyvinyl chloride foam, while h is representative of honeycomb core and is modeled as a transversely isotropic material with a significantly lower transverse shear stiffness than other materials.

EP

The RZT model [8] overcomes the problem suffered by RZT [48] of modeling laminates with external weak layers shown in [9], where the same benchmarks of [8] are considered. The RZT model [8] accounts for the geometric and constitutive heterogeneity of plies through a system of springs in series

AC C

combined with a system of springs in parallel that describe the transverse shear mechanics. The results of benchmarks T4 to T6 are reported in Table 3 and Figures 12 to 14 in the following normalized form:

106 h 2 w= 0 4 P Lx

L ∫−h / 2 w( 2x , z)dz ; σ xx = h/ 2

Lx , z) 2 ; P 0 L2x

h 2σ xx (

σ xz =

σ xz (0, z ) P0

;

σ zz =

Lx , z) 2 P0

σ zz (

(34) Figure 12 refers to case T4, while Figures 13 and 14 refer to cases T5 and T6, respectively. Table 3 reports the normalized maximum deflection, that is represented by the normalized average throughthickness deflection for the present models and for the exact solution, the maximum absolute axial stress and maximum absolute transverse shear stress, compared to those predicted by the RZT model [8]. The

29

ACCEPTED MANUSCRIPT results of Table 3 are reported as relative percentage errors with respect to Pagano’s exact 3D elasticity solution. For the present models, also the relative error on the transverse normal stress is reported, a piecewise variable transverse displacement being assumed instead a constant one as in [8]. The results by [8] are indicated HR-RZT, those by the present adaptive model as ADZZ and the ones by the discretelayer and sublaminate models obtained from it as ADL and ASB, respectively. Figures 12 to 14 qualitatively compare the through-thickness variation of the transverse shear stress

RI PT

predicted by the present adaptive zig-zag, sublaminate and discrete-layer models, by the RZT model [8] and the exact solution. The results by RZT are obtained capturing points from figures of paper [8] with a frame grabber, like for all reference models of previous benchmarks. The results of models being closer each other than in former cases, the errors inherent this procedure can affect the accuracy of reproduced curves at a larger extent. No result is presented for the variation of the transverse normal stress, like for

SC

the majority of benchmarks presented in this paper, it being not challenging since all models accurately capture its variation through the entire laminate thickness for all cases considered in [8].

M AN U

The results of Table 3 show that errors less than 3% , that is considered as the acceptable accuracy margin, affect the results of all models considered. Figures 12 to 14 show that also the through-thickness variation of the transverse shear stress is accurately captured by all models. Because the results by the RZT model [8], that are obtained neglecting the transverse normal deformability, are in a good agreement with the exact solution and with the results by the present models that account for it, this deformability is of marginal importance for benchmarks T4 to T6, whose material properties enhance just the effects of transverse shear deformability. On the contrary, for benchmark S6, where the elastic modulus E3 of the

TE D

core and of the upper face layers are reduced of a factor 10-2 and for benchmark S3, where a nonclassical boundary condition is considered, the transverse normal deformability was shown important.

8. Evaluation of errors in computing the strain energy density across the thickness

EP

Previous results have shown the importance of a high-order zig-zag representation of all displacement components across the thickness for obtaining accurate stress predictions or, as an alternative, techniques (e.g., developing models in mixed form) that allows to achieve accurate stresses also from a not extremely

AC C

accurate layerwise variation of displacements, owing to simplifying assumptions introduced. In applications of industrial interest such as optimization of variable stiffness composites [72] and damage evaluation [73], local energy needs to be exactly evaluated in order to obtain effective results. A measure of the errors involved by the theories previously discussed in the evaluation of the strain energy density is given in Table 4 where the following energy indicator:

EE ij

[σ ( x) =

ij

( z)

S αβ P

] 10 2

5

0

(35)

is reported . Sαβ denote the elements of the matrix in the inverse relation (1) that defines strains from stresses ( S11 is used for

σ xx , while S 44 is used for σ xz ). The energy indicator EE ij

is computed at

30

ACCEPTED MANUSCRIPT the points across the thickness where the largest errors are shown with respect to reference solutions. Previous numerical results have demonstrated that various theories accurately predict stresses. However, as shown by Table 4, even small errors in the computation of stresses turn into non-negligible errors on strain energy density, as shown by the index EE ij . A better agreement is shown with respect to exact solutions than to FEA results in the cases where exact results are not available. It is thought that these

RI PT

larger errors are inherent to the nature of finite elements, as their accuracy depends either by the accuracy of the element formulation or by the meshing adopted. This is what suggest the results by the present

adaptive model, because when compared to all other theories whose results are obtained in analytical

form it always gives lower errors than compared to finite element results. The errors of mixed models like RZT appear in some cases larger of those by other theories, because they predict the stresses very

SC

accurately but not always the displacements, as shown in the figures.

These errors show the importance of achieving the maximal accuracy also predicting displacements and not only stresses, for obtaining the best prediction of the energy stored locally, as required, e.g. for failure

M AN U

analysis through energy methods [73]. No result is given in Table 4 for the benchmark S1 because the exact solution [11] just shows the maximum values of displacements and stresses across the thickness, not their distributions. The theories that obtain accurate stresses but not equally accurate displacements are effective using stress-based failure criteria that just require accurate stress predictions. Use of strainbased or energy-based failure criteria needs theories that accurately predict displacements and stresses. Hence, Table 4 indirectly confirms the importance of a high-order representation of displacements like in

energy density.

TE D

[43], better estimates of strains being obtained to which correspond better local estimates of the strain

9 Concluding remarks

EP

The flexural response of laminated composite and sandwich beams and plates under static distributed loading has been examined using an equivalent single-layer zig-zag model (with variable through-

AC C

thickness kinematics, a piecewise variation also of the transverse displacement component and just five variables) that has recently developed by the authors and a discrete-layer model and a sublaminate models developed in this paper from it. These three laminated plate models a priori make continuous displacements and transverse shear stresses at the layer interfaces by deriving the expressions of zig-zag amplitudes in closed-form from the enforcement of the fulfillment of kinematic and transverse shear stress conditions at the layer interfaces. A symbolic calculus tool is used for deriving these closed-form expressions that overcomes the algebraic difficulties. The new contribution by these models is the fulfillment of the continuity of the transverse normal stress and of its gradient across the thickness prescribed by the elasticity theory but usually neglected by the models published in literature. The displacement fields of the present models consider a linear variation across the thickness as basic contribution and higher order terms whose coefficients are derived from the enforcement of stress

31

ACCEPTED MANUSCRIPT boundary conditions at upper and lower faces, like for equivalent single-layer models. Higher-order contributions are incorporated, along with zig-zag layerwise contributions, in order to enforce the fulfillment of local equilibrium equations at selected points (this enables obtainment of accurate stress predictions from constitutive equations) and to change the order of the representation with the purpose of locally improving accuracy and/or saving costs. It is still open to discussion whether an equivalent single-layer zig-zag model with a fixed number of

RI PT

degrees of freedom can predict through-thickness displacements, strain and stress fields as accurately as discrete-layer and sublaminate models for cases with strong transverse anisotropy and layerwise effects. In the present paper, the predictive capability of the equivalent single-layer zig-zag model is assessed

comparing its results to those of its counterpart discrete-layer and zig-zag models for benchmarks with

strong layerwise effects and transverse anisotropy taken from the literature with exact and approximate

SC

solutions by other researchers available for comparisons. Thick, simply-supported, cantilever and

propped- cantilever beams or plates under sinusoidal or uniform loading with geometric and constitutive

M AN U

heterogeneity of layers are analyzed. For these cases, for which discrepancies between the results of models and exact solutions are shown in literature, it is illustrated whether and when the transverse normal deformability is important.

The equivalent single-layer zig-zag model is shown as accurate as discrete-layer and sublaminate models and in a good agreement to exact solutions (beams/plates with simply-supported edges and sinusoidal transverse loading) or to FEA 3D (cantilever and propped-cantilever beams under uniform transverse loading). The advantage of the equivalent single-layer zig-zag model consists in providing very accurate

TE D

results even for benchmarks that are challenging for models always with a memory storage dimension and a processing time comparable to lowest ones of FSDT and HSDT equivalent single layer models. The transverse normal deformability of soft-core sandwich is shown important to obtain accurate stress predictions near the supported edge of the propped-cantilever sandwich beam (benchmark S3). Also benchmarks S6, a sandwich beam with damaged upper face and core (E3 reduced by a factor 10-2) and

EP

benchmark T2, a simply-supported rectangular cross-ply plate, demonstrated that a superior performance is obtained with a refined representation of the transverse normal deformability, nevertheless the material properties of constituent layers and the lay-up are not dissimilar to those of other cases for which

AC C

they are unimportant. Benchmark S6 demonstrates that the different thickness of face sheets and their different elastic properties give rise to strong layerwise and transverse anisotropy effects that make the behaviour of undamaged and damaged sandwiches totally different. The modulus E3 of faces and core is shown to have a significant bearing on stress and displacement fields. Benchmark T2 shows that a high order piecewise representation of displacements is required, because otherwise errors in all displacement and stress fields are made, though a correct trend is shown across the thickness. Laminate T1 and T4 to T6, as well as sandwich T3 are shown benchmarks for which a piecewise variable transverse displacement has a marginal importance. Benchmark S4, a simply-supported thick sandwich plate with thin faces show that the piecewise representation of the transverse displacement offered by the present zig-zag model and its capability to

32

ACCEPTED MANUSCRIPT adapt to the variation of solutions by changing the representation across the thickness improve precision also when transverse anisotropy and Layerwise effects are not extreme. Benchmark S5, a rectangular simply-supported sandwich plate with the faces having different thickness and material properties, shows that kinematic-based Murakami zig-zag function is less accurate than physically-based Di Sciuva’s zig-zag functions for a structure with asymmetric geometric and constitutive properties, confirming other similar results in literature.

RI PT

Even small errors in the computation of stresses are shown to turn into non-negligible errors on strain energy density.

References

SC

[1] Reddy JN, A generalization of two-dimensional theories of laminated composite plates. Commun Appl Numer Meth 1987; 3: 173-180. [2] Reddy JN, Robbins DH Jr, Theories and computational models for composite laminates, Appl Mech Rev 1994; 47 (6): 147-169.

M AN U

[3] Liu PF, Islam MM, A nonlinear cohesive model for mixed-mode delamination of composite laminates. Comp. Struct 2013; 106: 47-56. [4] Vachon PL, Brailovski V, Terriault P, Prediction of the propagation of impact-induced delamination in carbon/epoxy laminates. Comp. Struct 2013; 95: 227–235. [5] Lekhnitskii, S.G, Strength calculation of composite beams. Vestn Inzh Tekh, 1935. [6] Ambartsumyan S.A., On the theory of bending plates. Isz Otd TechNauk AN SSSR 1958; 5: 69-77.

TE D

[7] Ambartsumyan S.A., On a general theory of anisotropic shells, Prikl Mat Mekh 1958; 22: 226-237.

EP

[8] Groh RM Jr, Weaver PM, On displacement-based and mixed-variational equivalent single layer theories for modeling highly heterogeneous laminated beams. Int J Solids and Struct 2015; 59: 147170. [9] Gherlone M, On the use of zigzag functions in equivalent single layer theories for laminated composite and sandwich beams: a comparative study and some observations on external weak layers. ASME Appl. Mech. 2013; 80(6) : 1-19.

AC C

[10] Pagano NJ (1969) Exact solutions for composite laminates in cylindrical bending. J Composite Mater 1969; 3: 398–411. [11] Faraboschi P, Three-layered plate: Elasticity solution. Composites Part B 2014; 60: 764-776 . [12] Faraboschi P, Optimal design of glass plates loaded transversally. Mat and Des 2014; 62: 443-458. [13] Faraboschi P, Layered plate with discontinuous connection: exact mathematical model. Composites Part B 2013; 47,:365-378. [14] Faraboschi P, Hybrid laminated-glass plate: design and assessment. Composite Struct 2013; 106: 250-263. [15] Mattei O, Bardella L, A structural model for plane sandwich beams including transverse core deformability and arbitrary boundary conditions. Europ J o Mech Part A Solids, In press. [16] Icardi U, Sola F, Indentation of sandwiches using a plate model with variable kinematics and fixed degrees of freedom. Thin-Walled Structures 2015; 86: 24-34.

33

ACCEPTED MANUSCRIPT

[17] Noor AK, Burton WS, Bert CW (1996) Computational models for sandwich panels and shells. Appl Mech Rev 1996; 49: 155-199. [18] Noor AK, Malik M An assessment of modelling approaches for thermo-mechanical stress analysis of laminated composite panels. Comput Mech 2000; 25: 43–58.

RI PT

[19] Carrera, E., 2001, “Developments, ideas and evaluations based upon Reissner’s mixed variational theorem in the modeling of multilayered plates and shells. ASME Appl. Mech. Rev. 2001; 54(4): 301–329. [20] Ghugal, YM, Shimpi RP, A review of refined shear deformation theories of isotropic and anisotropic laminated plates. J Reinf Plast Comp 2002; 21:775- 813.

SC

[21] Reddy JN, Mechanics of laminated composite plates and shells: theory and analysis. 2nd Edition, CRC Press, Boca Raton, 2003.

M AN U

[22] Kreja I, A literature review on computational models for laminated composite and sandwich panels. Central European Journal of Engineering 2011; 1: 59 – 80. [23] Khandan R, Noroozi S, Sewell P, Vinney J, The development of laminated composite plate theories: a review. J Mater Sci 2012; 47: 5901-5910. [24] Yaqoob Yasin M, Kapuria S, An efficient layerwise finite element for shallow composite and sandwich shells. Comp. Struct 2013; 98: 202-214. [25] Plagianakos TS, Saravanos DA, Higher-order layerwise mechanics and finite element for the damped characteristics of sandwich composite beams. Int. J. Solids Struct 2004; 41 : 6853-6871.

TE D

[26] Plagianakos TS, Saravanos DA, Higher-order layerwise laminate theory for the prediction of interlaminar shear stresses in thick composite and sandwich composite plates. Comp. Struct 2009; 87: 23-35. [27] Di Sciuva M, A refinement of the transverse shear deformation theory for multilayered orthotropic plates. L’Aerotecnica Missili e Spazio 1984; 62: 84–92.

EP

[28] Murakami H, Laminated composite plate theory with improved in-plane responses. ASME Appl Mech 1986; 53: 661–666.

AC C

[29] Whitney JM, Pagano NJ, Shear deformation in heterogeneous anisotropic plates. ASME J Appl Mech 1970; 37: 1031-1036 [30] Reddy JN, A simple higher-order theory for laminated composite plates. ASME J Appl Mech 1984;51:745-752. [31] Murakami H, Maewal A, Hegemier GA, A mixture theory with a director for linear elastodynamics of periodically laminated media. Int. J. Sol. Str.1981; 17: 155–173. [32] Carrera E, A Priori vs. a posteriori evaluation of transverse stresses in multilayered orthotropic plates, Compos Struct 2000; 48: 245–260. [33] Carrera E, An assessment of mixed and classical theories on global and local response of multilayered orthotropic plates. Compos. Struct 2000; 50: 183–198. [34] Brischetto S, Carrera E, Demasi L, Improved response of asymmetrically laminated sandwich plates by using Zig-Zag functions. J. Sandwich Struct. & Mat 2009; 11: 257- 267.

34

ACCEPTED MANUSCRIPT

[35] Carrera E, Brischetto S, A survey with numerical assessment of classical and refined theories for the analysis of sandwich plates. Appl Mech Rev 2009; 62: 1–17. [36] Ferreira, AJM, Roque CMC, Carrera E, Cinefra M, Polit O, Radial basis functions collocation and a unified formulation for bending, vibration and buckling analysis of laminated plates, according to a variation of Murakami’s zig-zag theory,” Eur. J. Mech. A Sol. 2011; 30: 559–570.

RI PT

[37] Pagano NJ, Soni SR, Global-local laminate variational model. Int J Solids Struct 1983; 19: 207228.

[38] Averill RC, Yip YC, Development of simple, robust finite elements based on refined theories for thick laminated beams. Computers & Struct 1996; 59: 529-546.

SC

[39] Gherlone M, Di Sciuva M, Thermo-mechanics of undamaged and damaged multilayered composite plates: a sub-laminates finite approach. Computers & Struct 2007; 81:125–36.

M AN U

[40] Icardi U, Layerwise mixed element with sublaminates approximation and 3D zig-zag field, for analysis of local effects in laminated and sandwich composites. Int. J. Num. Meth. Eng. 2007; 70: 94-125. [41] Icardi U, Higher-order zig-zag model for analysis of thick composite beams with inclusion of transverse normal stress and sublaminates approximations. Comp. Part B 2001; 32: 343-354. [42] Zhen W, Wanji C, An assessment of several displacement-based theories for the vibration and stability analysis of laminated composites and sandwich beams. Comp. Struct. 2008; 84 : 337–349.

[43] Icardi U, Sola F, Development of an efficient zig-zag model with variable representation of displacements across the thickness. J. of Eng. Mech. 2014; 140: 531-541.

TE D

[44] Icardi U, Layerwise mixed element with sublaminates approximation and 3D zig-zag field, for analysis of local effects in laminated and sandwich composites. Int. J. Num. Meth. Eng. 2007; 70: 94-125

EP

[45] Zhen W, Wanji C, An efficient higher-order theory and finite element for laminated plates subjected to thermal loading. Comp. Struct. 2006; 73 : 99–109 [46] Saravanos DA, Heyliger PR, Mechanics and computational models for laminated piezoelectric beams, plates, and shells. Applied Mechanics Reviews 1999; 52 (10): 305-320.

AC C

[47] Plagianakos TS, Saravanos DA, Coupled high-order shear layerwise analysis of adaptive Sandwich composite beams with piezoelectric actuators and sensors. AIAA Journal 2005; 43 (4): 885-894. [48] Tessler A, Di Sciuva M, Gherlone M, A Refined Zigzag Beam Theory for Composite and Sandwich Beams. J. Compos. Mat. 2009; 43: 1051–1081. [49] Vel SS, Batra RC, Analytical solution for rectangular thick plates subjected to arbitrary boundary conditions. AIAA J. 1999; 37: 1464 – 1473. [50] Swaminathan K, Ragounadin D, Analytical solutions using a higher-order refined theory for the static analysis of antisymmetric angle-ply composite and sandwich plates. Comp. Struct. 2004; 64: 405 – 417. [51] Iurlaro L, Gherlone M, Di Sciuva M. A mixed cubic zigzag model for multilayered composite and sandwich plates including transverse normal deformability. In Proc.: 11th World Conference on Computational Methods, Barcelona, Spain 2014.

35

ACCEPTED MANUSCRIPT [52] Iurlaro L, Gherlone M, Di Sciuva M, Tessler A, Assessment of the Refined Zigzag Theory for bending, vibration, and buckling of sandwich plates: a comparative study of different theories. Composite Structures 2013; 106 777–792 [53] Rekatsinas CS, Nastos CV, Theodosiou TC, Saravanos DA, A time-domain high-order spectral finite element for the simulation of symmetric and anti-simmetric guided waves in laminated composite strips. Wave Motion 53, 2015 1-19

RI PT

[54] Icardi U, Eight-noded zig-zag element for deflection and stress analysis of plates with general layup. Comp. Part B 1998; 29:425-441. [55] Ascione L, Fraternali F, A penalty model for the analysis of composite curved beams. Comput Struct 1992;45:985–99.

SC

[56] Fraternali F, Reddy JN. A penalty model for the analysis of laminated composite shells. Int J Solids Struct 1993;30 (24): 3337–55.

M AN U

[57] Li W, Li QN, Jiang WS, Jiang L, Seismic performance of composite reinforced concrete and steel moment frame structures – state-of-the-art. Compos Part B: Eng 2011;42 (2):190–206. [58] Shanmugam NE, Lakshmi B. State of the art report on steel–concrete composite columns. J Construct Steel Res 2001;57 (10): 1041–80. [59] Spacone E, El-Tawil E. Nonlinear analysis of steel-concrete compos struct: state of the art. ASCE J Struct Eng 2004;130(2). Special issue: Compos Hybrid Struct159–168. [60] Yeoh D, Fragiacomo M, De Franceschi M, Boon KH. State of the art on timber concrete compos struct: literature review. J Struct Eng 2011;137(10):1085–95.

TE D

[61] Shariyat M, A generalized global-local higher order theory for bending and vibration analyses of sandwich plates subjected to thermo-mecanical loads. Int J Mech Sci 2010; 52: 495-514 [62] Auricchio F, Sacco E. Refined first-order shear deformation theory models for composite laminates. ASME Appl. Mech. 2003; 70: 381-390

EP

[63] Li XY, Liu D, Generalized laminate theories based on double superposition hypothesis. Int J Num Meth Eng 1997; 40:1197–212.

AC C

[64] Zhen W, Wanji C, An efficient higher-order theory and finite element for laminated plates subjected to thermal loading. Comp. Struct. 2006; 73 : 99–109

[65] Zhen W, Wanji C, A study of global–local higher-order theories for laminated composite plates. Comp. Struct 2007; 79 : 44–54. [66] Zhen W, Wanji C, A quadrilateral element based on refined global-local higher-order theory for coupling bending and extension thermo-elastic multilayered plates. Int J Solids Struct 2007; 44: 3187–3217. [67] Wanji C, Cheung YK, Zhen W, Augmented higher order global–local theory and refined triangular element for laminated composite plates. Comp Struct 2007; 81 : 341–352

[68] Zhen W, Wanji C, Xiaohui R, Refined global–local higher-order theory for angle-ply laminated plates under thermo-mechanical loads and finite element model. Comp Struct 2009; 88 : 643–658

36

ACCEPTED MANUSCRIPT [69] Jones RM, Mechanics of composite materials 2nd Ed, Taylor and Francis, 1999.

[70] Icardi U, Sola F, C0 Fixed Degrees of Freedom Zigzag Model with Variable In-Plane and Out-ofPlane Kinematics and Quadrilateral Plate Element. Journal Aerosp. Eng 2014; 28: 040141351 0401413514.

RI PT

[71] Icardi U, Sola F Analysis of bonded joints with laminated adherends by a variable kinematics layerwise model. Int. J. of Adhesion & Adhesives 2014; 50: 244-254.

[72] ] Icardi U, Sola F, Response of sandwiches undergoing static and blast pulse loading with tailoring optimization and stitching. Aerospace Sci. & Tech. 2014; 32 : 293-301.

AC C

EP

TE D

M AN U

SC

[73] Ladevèze P, Lubineau G, On a damage mesomodel for laminates: micromechanics basis and improvement. Mech. of Materials 2003; 35: 763-775.

37

ACCEPTED MANUSCRIPT Appendix A1: displacement field of models by other researchers Averill and Yip [38] The displacement field of a laminated beam is assumed as: k −1

m

i =1

j =1

u( x, z) = u 0 ( x) + z Ψ ( x) + z 2 β ( x) + z 3η ( x) + ∑ ξi ( x)( z − zi ) − ∑ ξ j ( x)( z − z j )

RI PT

w( x, z ) = w0 ( x)

(A1.1)

u 0 and w0 denote the in-plane and the transverse displacement components of the reference surface

zm , Ψ is the linear rotation of layer m + 1 (the layer above the reference surface zm ), β and η are

( z − zl ) , which

ξl

(l=i,j) is an additional kinematic variable associated with the zig-zag function

SC

higher-order rotations,

are constructed to vanish in an arbitrarily specified layer z j . The methodology for

zag functions vanish. The functions

M AN U

defining the zig-zag function lacks invariance with respect to the choice of the fixed layer at which zig-

β , η and ξl

are to be eliminated by enforcing the continuity of

transverse shear stresses at the layer interfaces and the transverse shear stress-free boundary conditions at the upper and lower interfaces transverse shear stress. The stress continuity is fulfilled in a limiting sense, as an exterior penalty function formulation is adopted from which the stress continuity is obtained as the

TE D

penalty parameter approaches a large value. With this technique, the expression of

ξl

is obtained in

terms of constants depending just on the shear stiffness and thickness of each layer. After the enforcement of stress boundary and interfacial contact conditions as outlined, the final expression of the displacement field (17) obtained is C°- continuous and has only the three engineering degrees of freedom

u 0 , w0 and

Ψ . The following representation is used in the sublaminate zig-zag model [38] for representing the

EP

displacement field across a generic sublaminate: ni

u( x, y, z ) = ub ( x, y) + zΨ x ( x, y) + ∑ ξk ( x, y )( z − zk )

AC C

k =1

ni

v( x, y, z) = vb ( x, y) + zΨ y ( x, y) + ∑ηk ( x, y)( z − zk )

 z w( x, y , z ) = wb ( x, y )  1 − h sb 

(A1.2)

k =1

  z   + wt ( x, y )     hsb 

where hsb is the thickness of the sublaminate; ub and vb are the axial displacements in x and y directions at the mid-plane of the sublaminate , i.e. z=0, respectively; Ψi is the rotation of the normal at z=0; wt and wb are the transverse deflections of the top and bottom surface of the kth sublaminate. The zig-zag layerwise contributions are ∑ξk (x, y) (z-zk) and ∑ηk (x, y) (z-zk), ξk and ηk being the zig-zag amplitudes computed imposing conditions the continuity of transverse shear stresses at the interfaces.

38

ACCEPTED MANUSCRIPT Brischetto et al. [34]

Each elastic displacement component is expanded as:

ui = ui0 + zui1 + z 2ui2 + ... + z N uiN + (−1) k ξ k uiz (u1 = u ; u 2 = v ; u3 = w)

(A1.3)

ui0 being the displacement ui on the reference plane Ω , i.e. for z = 0 , u i1 being an additional variable

are layer independent, while (−1)

k

Ω and u iN , uiz being additional displacement variables that

RI PT

representing the rotation of the normal to

ξ k is the Murakami’s zig-zag function ( ξ K = 2( z − zMk ) / h k ).

Applications in [34] refer to cases ED4, EDZ1 and EDZ3 that correspond to N=4 and N=3, respectively k

ξ k uiz ). The number of degrees of freedom for such

SC

(only latter two include the zig-zag contribution (−1) theories are: EDZ1=9, ED4=EDZ3=15.

M AN U

Faraboschi [11]

Exact governing equations for a three-layered plate with elastic interlayer within the two-dimensional elasticity framework are obtained under assumption that the behavior is linear-elastic, each layer is made of isotropic material and behaves according to the Kirchhoff-Love plate hypotheses, the displacements are continuous at the interfaces and the interlayer does not transmit any stresses through the cross

TE D

sections:

∂4w ∂4w ∂ 4 w 6ζGt  ∂ut ∂vt  6ζβGt  − + 2 + + + ∂x 4 ∂x 2∂y 2 ∂y 4 λαh3  ∂x ∂y  αh3

 ∂ 2 w ∂ 2 w  6P  2 + 2  = 3 ∂y  αh  ∂x

3 ∂ 2 ut ∂ 2 ut ut Gt ∂w ∂ 2 vt h ∂3w  h  ∂ w + υα + Gh + α + G − G + β + ( G + υα ) =0   t 2 2 ∂x 3  2 h ∂x ∂x∂y λh ∂x 2 ∂y 2  ∂x ∂y

α

3 ∂ 2ut ∂ 2 vt ∂ 2 vt v G ∂w h ∂3w  h  ∂ w + υα + Gh + ( G + υ + α ) + α + G − Gt t + β t =0   3 2 2 2 2 ∂y 2 ∂x∂y h ∂y λh ∂y ∂x   ∂x ∂y

AC C

EP

α

where

(A1.4)

β = 1 if the elastic interlayer undergoes shear distortion (foams, metals, polymers), β = 0

interlayer undergo anti-symmetric bending (periodic cellular cores, compliant isotropic materials);

if the

λ is

half the thickness of the interlayer; h is the thickness of upper and lower layers bonded to the interlayer;

υ

is the Poisson’s ratio of the layers; G ,

respectively;

Gt shear moduli of the layers and of the interlayer,

P surface load; ς = h + 2λ lever arm of the couple of in-plane forces; u,v, w components

39

ACCEPTED MANUSCRIPT of displacements in the x,y,z directions;

u t , v t in-plane components of the relative displacement of the

upper interface with respect to the immovable plate middle surface;

α = E /(1 − υ 2 ) .

Groh and Weaver [8]

RI PT

A refined version of the RZT model [48] is developed that overcomes the problem of modelling

laminates with externally weak layers. The displacement field is assumed like in the RZT model as the superposition of FSDT kinematics and a zig-zag field that enhances the representation:

SC

U α( k ) ( X , z ) = uα ( X ) + zθ α ( X ) + Φ α( k ) ( z )ψ α ( X )

U z ( X , z ) = w( x, y )

α = 1,2; 1 ≡ x, 2 ≡ y; X ≡ x, y ; like in [48] and [51] a piecewise constant variation

M AN U

Like in (A1.8)

(A1.5)

is adopted for enhancing in-plane displacements, the continuity of transverse shear stresses being enforced only on the zig-zag dependent contribution. The effective shear modulus

Gα is found using

the stiffness equation of a set of springs in series (its expression coincides with that of [48] reported above for the model by Iurlaro et al. [51]). The springs-in-series analogy is used to define an effective inplane stiffness E = 1 / t

N

∑t Q k

k

TE D

k =1

(k )

. The change in displacement slope Uα , z at layer interfaces and the

change in curvature of the layerwise transverse shear stress profile are ruled by the layerwise stiffness k

ratio Gα / Gαz and the layerwise in-plane stiffness ratio

EP

Iurlaro et al. [51]

e k = Q k / E , respectively.

A piecewise cubic/quadratic displacement field for in-plane/transverse components respectively is

AC C

postulated as:

U α( k ) ( X , z ) = uα ( X ) + zθ α ( X ) + z 2 χ α ( X ) + z 3ωα ( X ) + Φ α( k ) ( z )ψ α ( X ) U z ( X , z ) = H bw ( z ) wb ( X ) + H tw ( z ) wt ( X ) + H aw ( z ) w ( X )

(A1.6)

across the generic k-th layer, which constitutes a refined version of the RZT model by Tessler et al. [48]. Symbols uα ,

θα , Φ α( k )

and ψ α

(α = 1,2; 1 ≡ x, 2 ≡ y; X ≡ x, y) represent the in-plane

displacement, the rotation along the

β − axis ( β

zig-zag amplitude, respectively, while

χα

and

is perpendicular to

ωα

α

) , the zig-zag function and the

represent higher-order rotations accounting for the

deformation of the normal in thick laminates. Like in [48], as zig-zag function, a piecewise constant variation is adopted for enhancing in-plane displacements. In contrast to Di Sciuva’s zig-zag fucntion, it

40

ACCEPTED MANUSCRIPT does not require the transverse shear stresses to be continuous across the layer interfaces, the continuity being enforced only on the zig-zag dependent contribution. Relaxing the transverse shear stress continuity enforcement, the model allows for a piecewise constant distribution of transverse shear stresses that is accurate in an average sense. Accordingly, the slopes ∂Φ α / ∂z = (k )

difference between the transverse shear rigidity

β α are defined by the

Gαz of a layer, and the effective transverse shear rigidity −1

k αz

quantities that belong to the generic layer k of thickness

t k , while t is the laminate thickness.

w

wb , wt and w , respectively

SC

The top, bottom and average transverse displacements are indicated as w

RI PT

N   k of the laminate, i.e. βα = G / Gα − 1 and Gα = 1 / t ∑ t / Gαz  the superscript k denoting  k =1  k

w

and H b , H t and H a are expressed as

M AN U

1 1 3 1 1 3 3 3 Hbw(z) =− − z + 2 z2 Htw(z) = − + z + 2 z2 Haw(z) = − 2 z2 4 2h 4h 4 2h 4h 2 2h

(A1.7)

The transverse normal stress is assumed in the following form:

wherein

qTz = { qz(b ) , q z(t )

the bottom

(b )

and top

(t )

(A1.8)

TE D

σ za = P ( z ) qv + L ( z ) q z

} is the vector collecting the external transverse normal pressure applied at

plate surfaces, while

L( z ) that rule the shape of σ z across the thickness

EP

independent coefficients P( z ) and

qv is the vector collecting two kinematic variables-

1  z  1  z  L( z ) =   − 1,  + 1  2  h  2  h 

(A1.9)

AC C

  z  2   z  2   P ( z ) = h 2   − 1 , zh 2   − 1  ;   h    h   

The expression of stransverse shear stresses is derived by integration of first two local equilibrium equations under the cylindrical bending assumptions.

Mattei and Bardella [15]

The axial displacement

d x ( x) across the upper face, core and lower face layers of sandwich beams is

assumed as:

c t   d x ( x) = ul ( x) −Φ l ( x) z − − l  ; ( z ∈ lower skin) 2 2 

41

ACCEPTED MANUSCRIPT tu tl   uu ( x ) − Φ u ( x ) 2 + ul ( x ) + Φ l ( x ) 2  − Φ c ( x ) z ; ( z ∈ core)   c t   d x ( x) = uu ( x ) −Φ u ( x ) z + + u  ; ( z ∈ upper skin) 2 2 

(A1.10)

RI PT

ul ( x) , uu ( x) being the displacements at lower and upper faces of each layer of the three-layered beam; Φ l , Φ u , Φ c denote the rotation of lower and upper faces and core, respectively; tl , tu and c denote the c d z ( x) = w( x) + s( x) = vl ( x) ; ( z ∈ lower skin) 2 d z ( x) = w( x) + s ( x ) z ; ( z ∈ core)

c = vu ( x) ; ( z ∈ upper skin) 2

M AN U

d z ( x) = w( x) − s( x)

SC

thickness of faces and core. The transverse displacement across the three layers is assumed as:

(A1.11)

v( x) being the deflection at the mid-plane and vl ( x) , vu ( x) those at the lower and upper faces resulting from rotations

Φ l , Φ u , Φ c d and the transverse normal strain of core s( x) = ε zz .

TE D

Plagianakos and Saravanos [26]

The displacement field across the k-th discrete-layer (k=1,2,…,N) is assumed as:

u k ( x, y , ζ k ) = U k ( x, y ) Ψ1k (ζ k ) + U k +1 ( x, y ) Ψ2k (ζ k ) + α xk ( x, y ) Ψ3k (ζ k ) + λ kx Ψ4k (ζ k )

EP

v k ( x, y, ζ k ) = V k ( x, y)Ψ1k (ζ k ) + V k +1 ( x, y)Ψ2k (ζ k ) + α yk ( x, y)Ψ3k (ζ k ) + λky Ψ4k (ζ k )

AC C

w k ( x, y , ζ k ) = w 0 ( x, y )

(A1.12)

The superscripts k and 0 denote the discrete layer and the mid-plane, respectively. Following expressions hold Ψ1 = (1 − ζ k ) / 2 , Ψ2 = (1 + ζ k ) / 2 , k

k

Ψ3k = hk / 2(ζ k − 1) , Ψ4k = hk / 2ζ k (ζ k − 1) ,

hk being the thickness of the discrete layer and ζ k =

2

z1k + z 2k 2 k k z− . Symbols z1 , z 2 denote the hk hk

thickness coordinates of bottom and top surfaces of the layer, while displacements of the layer bottom surface and

2

U k , V k denote the in-plane

U k +1 , V k +1 their counterparts of the top surface that

describe extension and rotation of the layer. The last two terms of in-plane displacements, that describe quadratic and cubic variations, vanish at the bottom and top surfaces of the layer to satisfy kinematic

42

ACCEPTED MANUSCRIPT continuity across the thickness of the laminate. Higher-order terms

α xk , λkx , α yk , λky are additional

degrees of freedom of the discrete layer and represent hyper-rotations.

Zhen and Wanji [45]

RI PT

A layer-by-layer representation of in-plane displacement components is assumed as the sum of a global component for each displacement uG ,

vG , which is expanded up to the third order, and two local )k )k k k contributions for each component uL , vL , uL , vL . The mid-plane of the laminate z = 0 is assumed

ζk

(−1 ≤ ζ k ≤ 1) . The displacement field is:

M AN U

uk ( x, y, z) = uG ( x, y, z) + uLk ( x, y, z) + uˆLk ( x, y, z)

is measured from the mid-plane of the layer k ;

SC

as the reference plane, the local thickness coordinate

vk ( x, y, z) = vG ( x, y, z) + vLk ( x, y, z) + vˆLk ( x, y, z) wk ( x, y, z) = wG ( x, y, z)

(A1.13)

wherein the single contributions are expressed as

TE D

uG ( x, y, z) = u0 ( x, y) + zu1 ( x, y) + z 2u2 ( x, y) + z3u3 ( x, y) vG ( x, y, z) = v0 ( x, y) + zv1 ( x, y) + z 2v2 ( x, y) + z3v3 ( x, y) wG ( x, y, z) = w0 ( x, y)

2

EP

 2  2 z +z  z +z  uLk ( x, y, z) =  z − k +1 k  u1k ( x, y) +  z − k +1 k  u2k ( x, y) zk +1 − zk  zk +1 − zk   zk +1 − zk  zk +1 − zk 2

 2  2 z +z  z +z  v ( x, y, z) =  z − k +1 k  v1k ( x, y) +  z − k +1 k  v2k ( x, y) zk +1 − zk  zk +1 − zk   zk +1 − zk  zk +1 − zk

AC C

k L

3

 2 z +z  uˆLk ( x, y, z ) =  z − k +1 k  u3k ( x, y) zk +1 − zk   zk +1 − zk 3

 2 z +z  vˆ ( x, y, z) =  z − k +1 k  v3k ( x, y) zk +1 − zk   zk +1 − zk k L

(A1.14) being assumed

ζ k = ak z − bk and ak = 2 /( zk +1 − zk ) , bk = ( zk +1 + zk ) /( zk +1 − zk ) . After enforcing

the stress-free boundary conditions on transverse shear stresses and the continuity of transverse shear stresses at layer interfaces, the final displacement field is obtained in the following form:

43

ACCEPTED MANUSCRIPT u k = u0 + (1 Φ k ( z )u1k + 2 Φ k ( z )u 2k + 3 Φ k ( z )u3k ) + [Φ k2 ( z )u1 + Φ 3k ( z )u2 + Φ 4k ( z )u3 ] + Φ5 ( z ) w0, x k

v k = v0 + (1Ψ k ( z )v1k + 2Ψ k ( z )v 2k + 3Ψ k ( z )v3k ) + [ Ψ2k ( z )v1 + Ψ3k ( z )v2 + Ψ4k ( z )v3 ] +

Ψ5k ( z)w0, y wk = w0

(A1.15) k

k

k

k

k

k

k

k

k

k

RI PT

The expressions of i Φ , Φ i , Ψ , Ψi just depend on material properties and thickness of layers. The

AC C

EP

TE D

M AN U

SC

unknown variables are u 0 , u1 , u 2 , u3 , u1 , u 2 , u3 , v0 , v1 , v2 , v3 , v1 , v2 , v3 , w0 .

44

ACCEPTED MANUSCRIPT Appendix A2: expressions of variable contributions ∆i The expressions of coefficients of high order contributions to the displacement field

Ax 2 ...Axn ,

Ay 2 ... Ayn , Az 2 ...Azn by Eq. (9) are derived as closed form expressions from the enforcement of local equilibrium equations (3) at selected points across the thickness. These terms are here referred as

RI PT

hierarchical terms, as they enable a variable representation in different regions across the thickness. The fulfilment of local equilibrium equations (3) should be enforced at np=Nlay · ord_u – 2 points across the thickness, where Nlay is the number of computational layers, while ord_u is the order of the expansion chosen for in-plane displacements.

To illustrate the procedure, consider a displacement field with a through-thickness representation for the

SC

in-plane displacement components u and v up to the third order and for the transverse component w up to the fourth order. The enforcement of the stress-free condition (4) for the transverse shear stress

σ xz

at

M AN U

the upper and lower faces, respectively, give rise to the following equations:

2 2 3 4 ni  0 i h i h k i h i h i h i h C γ x + 2C x + 3Dx   + ∑ Φ x H k + b, x + c, x   + d, x   + e, x   + 2 2  2  k =1 2 2 2  2 ni ni ni  h   k k h + ∑ Ψ , x ( x, y )  − zk  H k + ∑ Ω, x  − zk  H k + ∑ Cwk , x H k  + 2  2  k =1 k =1 k =1  2 2 3 4 ni  0 nl i h i h k i h i h i h i h +C45 γ y + 2C y + 3Dy   + ∑ Φ y H k + b,y + c,y   + d ,y   + e,y   + 2 2  2  k =1 2 2 2  2 ni ni ni  h   k k h + ∑ Ψ ,y ( x, y )  − zk  H k + ∑ Ω,y  − zk  H k + ∑ Cwk ,y H k  = 0 2  2  k =1 k =1 k =1 

EP

TE D

nl 44

2 2 3 4  0  h i h i h i h i h i  i h C γ x − 2C x + 3Dx   − b, x + c, x   + d , x  −  + e, x    + 2 2 2 2  2  2    2 2 3 4  0  1 i h i h i h i h i h i h +C45 γ y − 2C y + 3Dy   + b,y + c,y   + d ,y   + e,y    = 0 2 2 2 2 2  2   

AC C

1 44

(A2.1)

Similar equations are obtained from the enforcement of the stress-free boundary conditions for the other transverse shear stress

σ yz . By enforcing the transverse normal stress σ zz to equal the applied

distributed load acting on the upper and lower faces of the laminate, the following equations are obtained:

45

ACCEPTED MANUSCRIPT 2 3 ni ni  0 h 0   i h i h k h 0 C u,x + (γ x,x − w ,xx ) + Cx,x   + Dx,x   + ∑Φx,x  − zk  Hk + ∑Cuk,x Hk  + 2  2  2  k =1 2  k =1   nl 13

2 3 ni ni 0 h 0  h  i h i  h 0 +C v,y + (γ y,y − w ,yy ) + Cy,y   + Dy,y   + ∑Φky, y  − zk  Hk + ∑Cvk, y Hk  + 2  2  2  k=1 2  k =1   2 3 ni ni i h   i h i h i  h k +C b + 2c + 3d   + 4e   + ∑Ψ Hk + 2∑Ωk  − zk  Hk  + 2  2  2  k=1 2   k =1  nl 33

RI PT

nl 23

2 3 ni ni  0 h 0 h  i  h i  h 0 +C u,y + (γ x,y − w ,xy ) + Cx,y   + Dx,y   + ∑Φkx,y  − zk  Hk + ∑Cuk,yHk + 2  2  2  k=1 2  k =1  2 3 ni ni  h 0   h i  h i  h k h 0 0 +v,x + (γ y,x − w ,yx ) + Cy,x   + Dy,x   + ∑Φy,x  − zk  Hk + ∑Cvk,x Hk  = p0   2  2  2  k=1 2   2 k =1 

(A2.2)

M AN U

2 3  0 h 0 h  i h i  0 C u, x − (γ x , x − w , xx ) + C x , x   + Dx , x  −   + 2 2  2    2 3  0 h 0  i h i h 1 0 +C23  v,y − (γ y,y − w ,yy ) + Cy,y   + Dy,y    + 2 2  2   

SC

nl 36

1 13

2 3  i h  i h ih i +C b − 2c + 3d   + 4e  −   + 2 2  2   

TE D

1 33

2 3  0 h 0 h i h i  0 +C u,y − (γ x ,y − w , xy ) + Cx ,y   + Dx ,y  −  + 2 2  2  2 3 h 0 h  h i h i  0 + v,x0 − (γ y,x − w0,yx ) + Cy,x + D − y,x      = p −  2 2  2    2

EP

1 36

σ zz, z to vanish at the upper and lower faces, the equations obtained are:

AC C

while by enforcing the gradient

(A2.3)

2 2   h h  1  0 0 i h i h C131 (γ x0, x − w0, xx ) − 2Cxi , x + 3Dxi , x    + C23 ( γ − w ) − C + 3 D  y,y ,yy y,y y,y    + 2 2  2    2     2 2  i  1  0 1 i h i h 0 i h i h +C33 2c − 6d +12e    + C36 (γ x,y − w , xy ) − 2Cx,y + 3Dx,y   + 2 2  2   2   2

0 i h i h +(γ y,x − w0,yx ) − 2Cy,x + 3Dy,x   2  2

  =0 

46

ACCEPTED MANUSCRIPT 2 ni  0  0 i h i h C (γ x, x − w , xx ) + 2Cx, x + 3Dx, x   + ∑Φkx, x Hk  + 2  2  k =1   nl 13

2 ni  0  0 i h i h +C (γ y,y − w ,yy ) + 2Cy,y + 3Dy,y   + ∑Φky, y Hk  + 2  2  k =1   nl 23

2 ni  i  i h i h +C 2c + 6d +12e   + 2∑Ωk Hk  + 2 2 k =1   2 ni  0 0 i h i h +C (γ x,y − w , xy ) + 2Cx,y + 3Dx,y   + ∑Φkx,y Hk + 2  2  k =1  2 ni  0 0 i h i h +(γ y,x − w ,yx ) + 2Cy,x + 3Dy,x   + ∑Φky,x Hk  = 0 2  2  k =1 

RI PT

nl 33

SC

nl 36

(A2.4)

following equations:

M AN U

For example, the enforcement of equilibrium equations at a point zp in the layer j of a beam leads to j −1 j   C11j U ,ixx +C 44j U ,izz + C13j + C 44j W ,ixz = −C11j  U ,0xx + ∑ Φ,kxx ( z p − zk ) + ∑ Cu ,kxx  + k =1 k =1   j −1 j −1   − C13j + C 44j  ∑ Ψ,kx +2∑ Ω,kx ( z p − zk )  − C 44j U ,0zz k =1  k =1 

(A2.5)

q q   C 44j U ,izz +W ,ixz + C13j U ,ixz +C 33j W ,izz = −C 44j  U ,0zz + ∑ Ψ ,kx +2∑ Ω,kx ( z p − zk )  + k =1 k =1   q q   −C13j  U ,0xz + ∑ Φ x ,kx  − 2C 33j ∑ Ω k k =1 k =1  

(A2.6)

(

)

)

EP

(

TE D

(

)

The system of equations obtained from the enforcement of local equilibrium equations is solved using symbolic calculus (MATLAB® symbolic software package). Please note that, the solution is obtained

AC C

starting from the bottom layer of the structure and then moving up across the thickness. It is also remarked that all the algebraic operations required for determining the expressions of the coefficients of this appendix and those of next two appendices are carried out once for all. The location of the np points being arbitrary, it was chosen to enforce local equilibrium equations at all points across the thickness of the structure where the largest errors are shown. Errors are is preliminary checked once for all for each problem (i.e. lay-up, constituent materials, loading and boundary conditions) starting from points initially uniformly distributed across the thickness. If the residual of each equation is equal or lower than a fixed tolerance the point is chosen, otherwise it is moved to another position.

47

ACCEPTED MANUSCRIPT Appendix A3: expressions of continuity amplitudes Φxk, Φyk, Ψk, Ωk Following expressions of strain tensor components hold for the displacement field (7): ni

ni

ε xx ( x, y , z ) = U + U + ∑ Φ ( z − zk ) H k + ∑ Cuk, x H k i ,x

0 ,x

k x, x

k =1

k =1

ni

ni

ε yy ( x, y, z ) = V,y0 + V,yi + ∑ Φ ky , y ( z − zk ) H k + ∑ Cvk, y H k k =1

RI PT

k =1

ni

ni

ε zz ( x, y, z ) = W, zi + ∑ Ψ k H k + 2∑ Ω k ( z − zk ) H k k =1

k =1

ni

ni

k =1

k =1

ε xz ( x, y, z ) = U ,z0 + U ,zi + ∑ Φ kx H k + W, x0 + W,xi + ∑ Ψ ,kx ( z − zk ) H k + ni

ni

k =1

k =1

ni

ni

k =1

k =1

SC

+ ∑ Ω,kx ( z − zk ) 2 H k + ∑ Cwk , x H k

ni

ni

k =1

k =1

+ ∑ Ω,yk ( z − zk ) 2 H k + ∑ Cwk ,y H k ni

M AN U

ε yz ( x, y, z ) = V,z0 + V,zi + ∑ Φ ky H k + W,y0 + W,yi + ∑ Ψ ,yk ( z − zk ) H k +

ni

ε xy ( x, y , z ) = U ,y0 + U ,yi + ∑ Φ kx ,y ( z − zk ) H k + ∑ Cuk,y H k + k =1

k =1

ni

ni

k =1

k =1

TE D

+V,x0 + V,xi + ∑ Φ ky ,x ( z − zk ) H k + ∑ Cvk,x H k

(A3.1)

ùk being the number of interfaces crossed at z starting from the laminate lower face (k is the interface from layers k and k+1). From previous expressions of strains and from stress-strain relations (1), the expression of stresses σxz, σyz, σzz, and of the stress gradient σzz,z are in the generic layer k. The following

EP

analytic expressions of the continuity amplitudes are obtained for the layer k from the enforcement of stress compatibility equations (4), (5):

(

) + Φ (U + U ) + Φ (U + U ) + Φ (V + V ) + Φ (V + V ) + Φ (V + V ) + Φ (W + W ) + Φ (W + W ) + + Φ (W + W ) Φ = Φ (U + U ) + Φ (U + U ) + Φ (U + U ) + Φ (V + V ) + Φ (V + V ) + Φ (V + V ) + Φ (W + W ) + Φ (W + W ) + + Φ (W + W )

Φ kx = Φu1 U 0 + U i (k )

0

(k )

i

,y

u5

(k )

0

0

i

,x

0

(k )

i

,y

v5

(k )

,y

,z

0

(k )

i

,z

u3

(k )

i

u6

v1

v9

0

(k )

i

0

(k )

i

,x

u7

0

i

,x

u4

0

+

i

(A3.2)

,y

u8

,z

(k )

(k )

0

u2

i

u9

k y

(k )

,x

AC C

(k )

0

v6

(k )

0

0

(k )

i

,y

v2

(k )

i

,z

0

v7

0

(k )

i

,z

v3

(k )

i

,x

0

i

,x

v4

v8

0

i

,y

+ (A3.3)

i

,z

48

ACCEPTED MANUSCRIPT

(

) + Ψ (U + U ) + Ψ (U + U ) + Ψ (V + V ) + + Ψ (V + V ) + Ψ (V + V ) + Ψ (W + W ) + Ψ (W + W ) + (A3.4) + Ψ (W + W ) Ω = Ω (U + U ) + Ω (U + U ) + Ω (U + U ) + Ω (U + U ) + + Ω (U + U ) + Ω (U + U ) + Ω (V + V ) + Ω (V + V ) + + Ω (V + V ) + Ω (V + V ) + Ω (V + V ) + Ω (V + V ) + (A3.5) + Ω (W + W ) + Ω (W + W ) + Ω (W + W ) + Ω (W + W ) + + Ω (W + W ) + Ω (W + W ) Ψ k = Ψ1 U 0 + U i (k )

(k )

0

,y

(k )

0

,z

(k )

i

, xx

,z

0

,x

0

i

,x

4

(k )

i

7

, yz

0

, xz

0

i

,y

8

0

i

, xx

13

, yy

0

(k )

i

, yz

0

i

, zz

8

(k )

i

, yz

(k )

, xy

18

, xx

0

15

0

12

, xz

i

, yy

0

i

, xy

0

i

, zz

(k )

i

0

4

(k )

i

0

11

i

14

, xz

7

(k )

i

(k )

i

0

16

i

, yy

SC

0

, zz

0

0

3

(k )

i

0

10

(k )

(k )

i

, xy

6

(k )

i

0

2

(k )

i

0

9

17

(k )

i

(k )

i

RI PT

0

5

(k )

0

0

3

,z

(k )

(k )

,y

6

1

(k )

(k )

i

i

9

(k )

0

2

(k )

i

5

k

(k )

,x

which involve derivatives of various order of the functional d.o.f. and unknown coefficients here termed as the continuity constants. At every interface k, 45 unknowns Φu1 (k), …., Φu9 (k), Φv1 (k), …., Φv9 (k), Ψ1 , …., Ψ9 (k), Ω1 (k), …., Ω18 (k) should be computed. The arbitrary nature of displacements enables the

M AN U

(k)

decomposition of stress contact conditions (4), (5), which turns (A3.2) to (A3.5) into 45 independent conditions at each interface, since the homologous displacements and displacement derivatives have to be collected separately at both hands as an independent condition to be fulfilled. This system of 45 equations in 45 unknowns is solved with the symbolic calculus tool of MatLab® without introducing any

AC C

EP

TE D

simplifying assumption.

49

ACCEPTED MANUSCRIPT Appendix A4: expressions of continuity amplitudes Cuk, Cvk , Cwk These expressions are obtained in a straightforward way from enforcement of kinematic continuity conditions (2) at the interfaces between the mathematical layers i and i+1 where representation (7) is changed as:

RI PT

C uk ( x, y ) = U i ( x, y , z ℑ ) − U i +1 ( x, y , z ℑ )

(A4.1)

C vk ( x, y ) = V i ( x, y, z ℑ ) − V i +1 (x, y , z ℑ )

(A4.2)

C wk ( x, y ) = W i ( x, y , z ℑ ) − W i +1 ( x, y , z ℑ )

SC

(A4.3)

Differently to stress continuity conditions, the former equations do not involve derivatives of the functional d.o.f. of any order. Considering for example a beam, (A4.1) and (A4.3) write in explicit form

ℑ

M AN U

at the generic interface ℑ from computational layers i and i+1:

ℑ

U 0 ( x, y, z ℑ ) + U i ( x, y, z ℑ ) + ∑ Φ kx ( x, y )( z ℑ − z k ) H k + ∑ C uk ( x, y ) H k = k =1

= U ( x, y , z ℑ ) + U 0

i +1

k =1

ℑ+1

( x, y , z ℑ ) + ∑ Φ k =1

ℑ

k x

ℑ+1

(A4.1’)

( x, y )( z ℑ − z k ) H k + ∑ C ( x, y ) H k k u

k =1

ℑ

ℑ

W 0 (x, y, zℑ ) +W i (x, y, zℑ ) + ∑Ψk (x, y)(zℑ − zk )Hk + ∑Ωk (x, y)(zℑ − zk ) 2 Hk +∑Cwk (x, y)Hk = k =1

TE D

k =1

ℑ+1

k =1

ℑ+1

(A4.3’)

ℑ+1

= W (x, y, zℑ ) +W (x, y, zℑ ) + ∑Ψ (x, y)(zℑ − zk )Hk + ∑Ω (x, y)(zℑ − zk ) Hk +∑C (x, y)Hk 0

i

k

k =1

2

k =1

k w

AC C

EP

k =1

k

50

ACCEPTED MANUSCRIPT

mm

2500

2750

3000

5600

6250

8800

12500

B

“

2500

1750

3000

4400

6350

6200

12500

h

“

75

18.5

8.75

32.5

100

65

115

2λ

“

22

8

2.5

12

21.5

44

18.5

28

650

2

N/mm

750

8750

525

45000

37.5

E

“

1850

105000

207000

150000

18.65

υ

/

0.18

0.24

0.33

0.26

0.11

kN/mm

12

8.25

4

9.25

7.5

mm

2.3657

0.7518

8.5850

3.5240

3.2606

2.3794

0.7541

0.5824

3.5257

3.2739

2.3798

0.7543

0.5839

3.5268

2.3796

0.7543

0.5836

3.5263

0.6646

6.2808

26.6102

12.0765

ADZZ

0.6714

6.2931

26.6075

ADL

0.6741

6.3007

ASB

0.6727

ADZZ ADL ASB

σ Max

tum

2

N/mm

N/mm2

ADZZ ADL ASB

tvm ADZZ

2

N/mm

0.15

0.22

10

7.25

28.033

33.5425

28.420

33.8657

3.2745

28.435

33.8719

3.2738

28.428

33.8616

1.5450

5.8229

4.9602

12.0982

1.5502

5.8582

5.0075

26.6223

12.1004

1.5505

5,8654

5.0089

6.2899

26.618

12.0986

1.5501

5.8576

5.0074

0.0818

0.1565

0.2810

0.2575

0.0898

0.1610

0.1729

0.8146

0.1573

0.2825

0.2581

0.0906

0.1619

0.1745

0.8209

0.1574

0.2829

0.2584

0.9049

0.1621

0.1746

0.8191 0.0818 0.8146 0.8209

ASB

0.8191

0,1573

0.2829

0.2383

0.9046

0.1619

0.1746

0.1934

0.2810

0.2897

0.0898

0.1898

0.1729

0.1936

0.2826

0.2581

0.0906

0.1908

0.1745

0.1938

0.2829

0.2584

0.9049

0.1910

0.1746

0.1937

0.2828

0.2383

0.9046

0.1908

0.1745

AC C

EP

ADL

SC

wmax

16000

M AN U

P0

2

10000

TE D

G

RI PT

L

Table 1. Three-layered plate benchmark: comparison among the elasticity solution by Faraboschi [11] and approximate solutions by the present adaptive zig-zag (ADZZ), discrete-layer (ADL) and sublaminate (ASB) models.

-3

E1 [GPa]

172.4

0.1

1.72 10

E2 [GPa]

6.89

0.1

6.89 10-5

0.1

6.89 10

-5

3.45 10

-5

3.45 10

-5 -5

E3 [GPa] G12 [GPa] G13 [GPa]

6.89 3.45 3.45

0.04 0.04

G23 [GPa]

1.378

0.04

1.38 10

υ12

0.25

0.25

υ13

0.25

υ23

0.25

Gr-Ep

MAT bl

MAT g

137.92

5.9

132.38

5.512

10

10.76

5.512 2.76 2.76

10 5.9 0.2

Foam

1.076

Carbon-

PVC

MAT

MAT

MAT

Epoxy

foam

1

2

3

0.035

157.9

0.104

1

33

25

0.05

0.035

9.584

0.104

1

1

1

0.05

0.035

5.65

0.0123

5.65

0.0123

9.584 5.93 5.93

SC

MAT bc

0.104 0.04 0.04

1

0.2 0.2

M AN U

MAT c

1

0.8 0.8

1 0.5 0.5

p

pvc

h

25x106 32.57x106

25x104

250

1x106 1x106

25x104

250

MAT 4

0.05 0.0217 0.0217

m

1x10

6

5x10

5

5x10

5 5

1x10

6

6.5x10

25x10 5

8.21x10

6

3.28x10

6

4

2500

9.62x10

4

9.62x10

4

875

9.62x10

4

1750

1

1.1024

0.7

3.61

0.0123

3.227

0.04

0.2

0.8

0.5

0.0217

2x10

0.25

0.25

0.25

0.24

0.4

0.32

0.3

0.25

0.25

0.25

0.15

0.25

0.25

0.3

0.9

0.25

0.25

0.25

0.25

0.24

0.4

0.32

0.3

0.25

0.25

0.25

0.15

0.25

0.25

0.3

3x10-5

0.25

0.25

0.25

0.25

0.49

0.4

0.49

0.3

0.25

0.25

0.25

0.15

0.25

0.25

0.3

3x10-5

TE D

MAT p

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

Table 2. Mechanical properties of laminates and sandwiches. Properties of material p,m,pvc and h are normalized by the-in-plane shear modulus

G12( h) of material h

ACCEPTED MANUSCRIPT

T5

σ xz

σ zz

HR-RZT

-0.020

0.020

0.040

---

ADZZ

-0.011

0.01

0.036

-0.005

ASB

-0.018

-0.009

0.052

-0.086

ADL

-0,015

0.009

0.041

-0.003

HR-RZT

0.360

-0.94

0.060

ADZZ

-0.064

-0.071

0.073

ASB

-0.087

-0.093

0.077

ADL

-0.062

-0.058

0.071

-0.120

---

0.136

-0.159

0.141

-0.127

0.123

-0.124

HR-RZT

0.290

0.610

ADZZ

-0.049

-0.098

ASB

-0.066

-0.128

ADL

-0.572

0.074

M AN U

T6

σ xx

RI PT

T4

w

Model

SC

Laminate

---

-0.184 0.136 0.128

Table 3. Normalized maximum deflection, maximum absolute axial stress and maximum absolute transverse shear stress, for benchmarks T4, T5 and T6, as predicted by the present adaptive (ADZZ),

AC C

EP

with respect to exact solution

TE D

sublaminate (ASB) and discrete-layer (ADL) models, by the RZT model [8] (HR-RZT) as relative percentage errors

ACCEPTED MANUSCRIPT

% error

Sandwich plate S4

z= 0.41h

EEyz Exact

7.6720

EEyz theory [43]

7.6424

-0.38%

EEyz Global local [45]

8.6609

12.89%

Sandwich plate S5

z=-0.46 h

0.19284

RI PT

EExz Exact EExz theory [43]

0.21749

EExz theory [34]

0.46332 140.25%

EExz Exact EExz theory [43]

z= 0.03h

794.098

1.19%

EExz Sublaminate theory [38]

614.882

-21.65%

EExz Global local theory [45]

761.419

-2.98%

EExz FEM 3D

Laminated plate T2

z= -0.4 h

z=0

Sandwich plate T3

z=0.46

T5

AC C

Laminate J

z=0.3

EP

Laminate G T4

2.2303

12.54%

EExz RZT theory [48]

2.3849

20.34%

EExz Exact

3.4885

EExz theory [43]

3.4821

-0.18%

EExz RZT theory [51]

4.2254

21.12%

EExz Exact

1.3790

EExz theory [43]

1.3798

0.06%

EExz theory [26]

1.2428

-9.87%

TE D

z=0.04

EExz theory [43]

1.9818

M AN U

Cantilever sandwich beam S2

784.786

SC

Laminated plate T1

Laminate L T6

Propped beam S3 at x/Lx = 7/8

Propped beam S3 at x/Lx = 1

z=0

z=0.2

12.78%

EExz Exact

1.8362

EIxz theory [43]

1.8731

2.01%

EIxz RZT theory [8]

1.7761

-3.27%

EExz Exact

3.8019

EExz theory [43]

3.8194

0.46%

EExz RZT theory [8]

3.7671

-0.92%

EExz Exact

2.5633

EExz theory [43]

2.6148

2.01%

EExz RZT theory [8]

2.6354

2.81%

EExz Reference solution [15]

0.2525

EExz theory [43]

0.2588

EExz Reference solution [15]

2.3256

EExz theory [43]

2.3639

2.50%

1.65%

Table 4. Errors in the strain energy density computation by the theories considered in this paper with respect to exact or reference solutions. Index EEij computed where the largest error on stress

σ ij

is made across the thickness.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 1. a) Reference mid-plane

zk

Ω and coordinate system of the multi-layered plate (k denotes the layer numbering,

the position of the layer across the thickness ; b) through-the-thickness slope changes of the normal as postulated by

Murakami’s kinematic-based zig-zag function (ZZ-K) and Di Sciuva’s physically based zig-zag function (ZZ-F). On the left: the through-thickness variation of the normal as predicted by first- (FSDT) [29] and third-order (HSDT) [30] shear deformation C° theories .

Figure 2. Positive sense of bending

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

w,x0 (counter clockwise)

and shearing rotations

γ x0

(clockwise) of the normal. The

horizontal dashed line denotes the reference mid-plane, the bold lines the upper and lower faces and parallel lines the layer interfaces. The slope of the normal is changed at the interfaces by the zig-zag contributions to displacements as in

AC C

EP

TE D

Figure 1b.

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 3. Scheme of the points across the thickness at which displacement and stress continuity conditions, stress boundary conditions and equilibrium equations are enforced for the discrete-layer model of this paper. J and M are

AC C

EP

external layers whose faces constitutes the upper and lower bounding faces of the laminate, L is an intermediate layer.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 4. Through-thickness transverse shear stress

σ xz

and the in plane stress

σ xx

distributions at

(Lx/5, Ly/2) for a cantilever, uniformly loaded sandwich beam with a length-to-thickness ratio Lx/h=10, as predicted by the present adaptive, discrete-layer and sublaminate models, the 3D FEA solution [48] and the RZT theory [48]

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 5. Through-thickness distribution of normalized shear stress

σ xz

at: a) x/Lx=7/8 and x/Lx =0 (clamped edge); b)

at x/Lx=1 (supported edge), for a propped-cantilever beam, as predicted by the reference solution [15] and by the present zig-zag adaptive, discrete-layer and sublaminate models (Lx/h=5.714). Solution S=0 by [15] disregards the transverse normal deformability.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 6. Through-thickness variation of transverse shear stresses

σ xz

(a) and

σ yz (b) for a simply-supported sandwich

plate undergoing bi-sinusoidal loading with Lx/h=4, as predicted by the present adaptive, discrete-layer and sublaminate models, by the global-local zig-zag model [45] and exact solution.

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 7. Through-thickness variation of the normalized transverse shear stress σ xz and in-plane displacement

u , as

EP

predicted by the present models, the zig-zag theory [34] based on Murakami’s zig-zag function (here termed geometrically-based) and the exact solution for a rectangular three-layer simply-supported sandwich plate (Lx/h=4,

AC C

Ly/Lx=3 ) under sinusoidal loading

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

TE D

Figure 8. Through-thickness stress and displacement fields for a simply-supported sandwich beam with undamaged (UNDAM) or damaged laminated faces ( E3 of upper face and core reduced by a factor 10-2) under sinusoidal loading,

AC C

(Lx/h=4)

EP

as predicted by the PWM 3 model [41], the present adaptive, discrete-layer and sublaminate models and exact solution

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 9. Comparison between transverse shear stress, in-plane stress and in-plane displacement predicted by the present adaptive, discrete-layer and sublaminate models, by the sublaminate theory[38], the global local theory [45] and the exact solution [38] for a simply supported [0°- 90°- 0°- 90°] laminated beam (Lx/h=4). The damaged case

AC C

EP

TE D

considers E3 of the upper layer reduced by a factor 10-2

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 10. Through-thickness distribution of displacements and stresses for simply-supported [0°/-90°/0°] rectangular plate with bi-sinusoidal loading. Results by the present adaptive, LIN-ADZZ, CUB-ADZZ, sublaminate and discretelayer models derived from it, the RZT theory with variable transverse displacement [51] and exact solution.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 11. Comparison between transverse shear stress and in-plane stress displacement by the present adaptive, discrete-layer and sublaminate models ( indicated as present model, present model with discrete-layer formulation and present model with sublaminate formulation, respectively), the discrete-layer model [26] and the exact solution for a sandwich plate with L/h=10.

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 12. Through-thickness distribution of the transverse shear stress for the simply-supported laminated beam T4

TE D

under sinusoidal loading, as predicted by the present adaptive, sublaminate and discrete-layer models, by the RZT

AC C

EP

model [8] (HR-RZT) and exact solution

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

TE D

Figure 13. Through-thickness distribution of the transverse shear stress for the simply-supported laminated beam T5 under sinusoidal loading, as predicted by the present adaptive, sublaminate and discrete-layer models, by the RZT

AC C

EP

model [8] (HR-RZT) and exact solution

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 14. Through-thickness distribution of the transverse shear stress for the simply-supported laminated beam T6

TE D

under sinusoidal loading, as predicted by the present adaptive, sublaminate and discrete-layer models, by the RZT

AC C

EP

model [8] (HR-RZT) and exact solution