Structural theory and finite element modelling of linear elastic sandwich beams subject to severe boundary conditions

Accepted Manuscript Structural theory and finite element modelling of linear elastic sandwich beams subject to severe boundary conditions Andrea Panteghini, Lorenzo Bardella PII:

S0997-7538(16)30154-1

DOI:

10.1016/j.euromechsol.2016.10.012

Reference:

EJMSOL 3372

To appear in:

European Journal of Mechanics / A Solids

Received Date: 9 August 2016 Revised Date:

17 October 2016

Accepted Date: 18 October 2016

Please cite this article as: Panteghini, A., Bardella, L., Structural theory and finite element modelling of linear elastic sandwich beams subject to severe boundary conditions, European Journal of Mechanics / A Solids (2016), doi: 10.1016/j.euromechsol.2016.10.012. 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

RI PT

Structural theory and finite element modelling of linear elastic sandwich beams subject to severe boundary conditions Andrea Panteghini, Lorenzo Bardella∗

M AN U

SC

Department of Civil, Environmental, Architectural Engineering and Mathematics (DICATAM) University of Brescia, Via Branze 43 – 25123 Brescia, Italy

Abstract

We further develop and improve a structural theory recently proposed by our group, with the aim of determining the simplest kinematics which allows the accurate modelling of any plane sandwich beam in the linear elastic regime. The model builds on Yu–Krajcinovic zig-zag warping, in which each

TE D

layer, of arbitrary thickness and modulus, is allowed to shear through an independent cross-section rotation. Moreover, the core kinematics is enriched by allowing for a quadratic variation along the core thickness of both the

EP

longitudinal and the transverse displacement components. By implementing the proposed theory in a structural finite element, we discuss the contribution to the modelling capability of each independent term entering the chosen

AC C

core kinematics. Such kinematics, along with a Jourawski-like approach to evaluate the shear stress, leads to a model which can accurately describe the stress state for any relative stiffness between the sandwich layers, even in ∗

Corresponding author, Tel.: +390303711238; fax: +390303711312 Email addresses: [email protected] (Andrea Panteghini), [email protected] (Lorenzo Bardella) URL: http://lorenzo-bardella.unibs.it/index_eng.htm (Lorenzo Bardella) Preprint submitted to European Journal of Mechanics -A/Solids

October 17, 2016

ACCEPTED MANUSCRIPT

the case of “severe boundary conditions”, including loading on a specific skin

RI PT

coupled with constraints realised, on certain cross-sections, on the opposite skin only. We demonstrate this claim by considering many benchmarks and

by a thorough comparison with the results obtained from continuum plane stress Finite Element (FE) simulations. From such a comparison we also

SC

clearly establish the superior computational efficiency of the new structural finite element with respect to the continuum FE analyses.

M AN U

Keywords:

Sandwich beams, Warping, Finite Element Method 1. Introduction

We propose a structural model for plane sandwich beams to accurately describe their linear elastic mechanical response due to bending accompanied

TE D

with non-uniform shear, also including eccentric axial loading. We focus on sandwiches with perfect interfaces, whereas we allow for any geometry (from thin to thick skins, Allen, 1969) and for any relative stiffness between core and skins, spanning the whole range from antiplane sandwiches (Allen, 1969) to

EP

cross-sections whose bending behaviour is well represented by the Bernoulli– Navier kinematics.

AC C

Our main concern is the capability of predicting the sandwich behaviour

under Severe Boundary Conditions (SBCs), as introduced by Mattei and Bardella (2016). Modelling the stress state ensuing from SBCs, for instance encompassing any case in which a support is realised on the skin opposite to that a transverse load is applied on, is a difficult task in the context of a structural theory, for the large stress gradients involved, combined with the

2

ACCEPTED MANUSCRIPT

lack of Saint-Venant’s principle if the core is not stiff enough. To address

RI PT

this issue, a suitably rich cross-section kinematics must be assumed, in order to accurately describe the warping.

In most of the literature, usually, each proposed structural model is reliable in a limited and problem-dependent range of relative stiffness (as for

SC

instance discussed in Serpilli and Lenci, 2008; Berdichevsky, 2010; Tonelli

et al., 2012). As an example relevant for the present study, the model of

M AN U

Frostig et al. (1992) may be excellent for sandwiches having thin skins and very soft core. In fact, in the theory of Frostig et al. (1992) the core is modelled as a plane stress continuum under the constraint of vanishing longitudinal stress (that is the antiplane sandwich assumption), implying shear stress independent of the through-the-thickness coordinate y. At the opposite end, the computationally far less expensive First Order Shear Deformation

TE D

theory (FOSD, Timoshenko, 1921; Bert, 1973; Gordaninejad and Bert, 1989; Bardella, 2008), relying on Bernoulli–Navier’s kinematics (and, therefore, on Saint-Venant’s principle), cannot account for SBCs and may be appropriate for not-too-soft cores, in cases in which warping may be neglected.

EP

In between the aforementioned approaches, there is a huge body of literature on structural theories to appropriately model sandwich and laminated

AC C

beams, with emphasis on the description of warping (see, e.g., Allen, 1969; Ghugal and Shimpi, 2001; Yu and Hodges, 2005; Tessler et al., 2009; Icardi and Sola, 2016, and references therein). Among such theories, we consider those in which the displacement field is just C 0 -continuous with respect to y, whereas they rely on a zig-zag warping due to the deformation jump at the interfaces. Structural models based on zig-zag warping are known to be

3

ACCEPTED MANUSCRIPT

able to accurately describe both the deflection and the longitudinal stress

RI PT

field, at least when SBCs are not concerned with. This may hold both for discrete layer theories, in which the number of structural unknowns is proportional to the number of layers (see, e.g., Yu, 1959; Heller, 1969; Krajcinovic,

1972, 1975; Sharma and Rao, 1982; Bardella and Tonelli, 2012; Galuppi and

SC

Royer-Carfagni, 2012; Bardella and Mattei, 2014; Bardella et al., 2014), and for equivalent single layer theories, aiming at maintaining the number of un-

M AN U

knowns independent of the number of layers (see Zuo and Hjelmstad, 1998; Ghugal and Shimpi, 2001; Tessler et al., 2009; Carrera et al., 2013; Icardi and Sola, 2016, and references therein). Though the distinction between discrete layer and equivalent single layer theories may be crucial for laminated structures, it is not much relevant for sandwiches.

It is well-known that the shear stress estimate is in general poor if simply

TE D

obtained on the basis of the zig-zag kinematics and the constitutive law. A shear stress recovery technique is therefore required (Matsunaga, 2002) and here we propose a Jourawski-like approach (Jourawski, 1856) which has been proved to provide accurate estimates in our previous investigations on other

EP

structural beam models (Bardella and Tonelli, 2012; Bardella and Mattei, 2014; Bardella et al., 2014; Mattei and Bardella, 2016).

AC C

In order to reach our goal, we build our structural model on the the-

ory of Mattei and Bardella (2016), that, in turn, is an extension of Yu– Krajcinovic zig-zag theory (Yu, 1959; Krajcinovic, 1972, 1975) to account for the transverse core strain, that is the normal strain εy in the core along the thickness direction y. Such extension has been proposed as a first step towards properly accounting for SBCs. Note that if the boundary conditions

4

ACCEPTED MANUSCRIPT

(both constraints and loads) are instead smoothly applied all along the cross-

RI PT

section, the Jourawski-like approach applied to the Yu–Krajcinovic zig-zag theory (neglecting εy ) allows one to obtain explicit analytic solutions that have been shown to be very accurate (Bardella and Tonelli, 2012; Bardella

and Mattei, 2014). In order to obtain the same accuracy when SBCs are con-

SC

cerned with, it is necessary to account for the transverse core strain. After

the seminal work of Frostig et al. (1992), other studies have recognised the

M AN U

importance of modelling the transverse core strain (see, e.g., Vidal and Polit, 2009; Phan et al., 2012; Jedari Salami et al., 2016; Mattei and Bardella, 2016, and references therein). Mattei and Bardella (2016) assume the transverse core strain to be independent of y. This leads to a model characterised by six independent structural functions, required to have all layers’ rotations independent of the beam axis transverse displacement. By resorting to the

TE D

Rayleigh–Ritz method to obtain numerical solutions, Mattei and Bardella (2016) have shown that their model can predict some relevant features of the flexure due to non-uniform shear under SBCs. Here, we improve the model of Mattei and Bardella (2016) by enriching the core kinematics. In fact, for

EP

both the displacement components of the core we consider second-order polynomials of the coordinate y, in which each monomial turns out to be given

AC C

by the product of a shape function and an independent structural function of the beam axis coordinate, x. Hence, the present model, with respect to that of Mattei and Bardella (2016), has two further unknown structural functions, having eight independent of them, overall. In order to obtain predictions of this novel structural model, we imple-

ment a new one-dimensional finite element, having 8 Degrees Of Freedom

5

ACCEPTED MANUSCRIPT

(DOFs) per node, henceforth referenced to as the “8DOFs model”. By con-

RI PT

straining the DOFs related to the quadratic shape functions in the core, our

implementation particularises to the model of Mattei and Bardella (2016) (having 6 DOFs per node, referenced to as the “6DOFs model”) and also

to a model in which only the transverse core displacement component is a

SC

quadratic function of y, while the longitudinal core displacement component is linear in y (that is, the “7DOFs model”). This allows a thorough discus-

M AN U

sion of the influence, on the model predictions, of each structural function added to the basic Yu–Krajcinovic zig-zag kinematics.

By applying the new model to benchmarks already analysed in literature and also to a new one, we demonstrate the modelling capability and computational efficiency of our proposal. In particular, for what concerns the boundary value problems already considered in literature, we focus on simply-

TE D

supported and propped-cantilever sandwiches subject to uniform transverse load on the skin opposite to that the supports are applied on. Instead, the benchmark introduced for the first time in this work is characterised by a highly inhomogeneous core deformability. The discussion is based on

EP

the comparison with reference results obtained by two-dimensional, plane stress, continuum Finite Element (FE) analyses and also with other litera-

AC C

ture models (FOSD theory; Krajcinovic, 1972; Frostig et al., 1992; Mattei and Bardella, 2016). We show that the present 8DOFs model is the sole one which can accurately predict the stress state in most of the beam for any sandwich under SBCs. In particular, with respect to the Mattei and Bardella (2016) 6DOFs model, the 8DOFs model is significantly more accurate both for soft cores and in the limit of homogeneous cross-section.

6

ACCEPTED MANUSCRIPT

Outline of the paper. In Sect. 2 we define the cross-section kinematics char-

RI PT

acterising the new structural theory, whose Total Potential Energy functional

is established in Sect. 3. Sect. 4 is devoted to the Jourawski-like analytic

estimates for the shear stress. In Sect. 5 we provide the details to implement the new one-dimensional finite element and in Sect. 6 we demonstrate the

SC

modelling capability of the proposed structural theory. Finally, in Sect. 7 we offer some concluding remarks.

M AN U

2. Sandwich beam kinematics and notation

The beam axis x is chosen to coincide with the core centre-line and y is the coordinate along the cross-section height. We consider the static plane stress mechanical response, in the (x, y) plane, of sandwich beams subject to concentrated moments, transverse and axial forces, and distributed loads.

TE D

We denote the core thickness as c and the skins’ thicknesses as tl and tu ; here and henceforth, the subscripts l and u are used to refer to the lower and upper layers, respectively. The symbols El and Eu and Gl and Gu are employed for

EP

the Young and shear moduli of the skins, respectively. Independently of the possible anisotropy of the skins, the Young modulus refers to the longitudinal modulus (along x), whereas the shear modulus refers to the xy plane. Though

AC C

we formulate the theory for anisotropic core behaviour, in the benchmarks we consider isotropic core of Young modulus Ec and Poisson’s ratio νc . Along with the shear and longitudinal normal strains and stresses in each

layer, the proposed model accounts for the transverse normal strain, εy , and stress, σy , in the core. This is necessary to accurately predict the stress state of sandwich beams under SBCs, mostly if the skins are significantly stiffer 7

ACCEPTED MANUSCRIPT

than the core, as in the case of cores constituted by foam. In Mattei and

RI PT

Bardella (2016) εy is assumed to be independent of y: εy (x, y) = εy (x) ≡ sC (x), with sC (x) denoting a structural function governing the mechanical problem. In this work, we show that, in order to obtain a model able to well

describe the sandwich flexure for almost any relative stiffness of the sandwich

SC

layers, one has to enrich the transverse core strain description, introducing a linear dependence of εy on y:

2y c

M AN U

εy (x, y) = sC (x) + sL (x)

(1)

in which sL (x) is a further structural function. Furthermore, with respect to the kinematics adopted in Mattei and Bardella (2016), we add a quadratic term for the longitudinal displacement in the core, governed by the structural function p(x).

1

The fundamental assumptions for the adopted kinematics,

TE D

relying on a zig-zag warping, are illustrated in Fig. 1. In particular, Fig. 1 displays the eight structural functions employed in this investigation as independent functions, since they are suitable DOFs for the FE implementation

ul (x) ,

EP

of the proposed structural theory: uu (x) ,

v(x) ,

φl (x) ,

φu (x) ,

vl (x) ,

vu (x) ,

p(x) (2)

AC C

in which ul (x) and uu (x) are the axial displacements of the centre-lines of 1

In this investigation we do not seek to provide the most efficient structural model, for

a given computational cost, to describe the sandwich behaviour under SBCs. If this were the goal, there would be several alternative models to analyse, by considering different shape functions (e.g., we might substitute any parabolic shape function in the core with a sinusoidal one). See, e.g., Tessler et al. (2009); Vidal and Polit (2009); Carrera et al. (2013); Icardi and Sola (2016), and references therein.

8

ACCEPTED MANUSCRIPT

the bottom and top layers, v(x) is the deflection of the sandwich longitudinal

RI PT

axis x, φl (x) and φu (x) are the rotations of the bottom and top faces, and vl (x) and vu (x) are the deflections of the lower and upper skins.

The relations between the alternative structural variables sC (x) and sL (x), introduced above in discussing the transverse core strain description, and the

 1 vl (x) − vu (x) , c

sL (x) =

 2 vl (x) + vu (x) − 2v(x) c

(3)

M AN U

sC (x) =

SC

transverse displacements chosen as independent variables are:

The inverse of expressions (3) are reported in Fig. 2, together with other structural functions that are conveniently defined in order to write, in the following, clearer and more compact relations. Such functions are the axial displacements of the lower and upper interfaces, ucl (x) and ucu (x), respec-

TE D

tively, and the average core rotation, φc (x), defined as: ucl (x) = ul (x) + φl (x)

tl 2

tu 2  1 φc (x) = ucu (x) − ucl (x) c

EP

ucu (x) = uu (x) − φu (x)

(4) (5) (6)

AC C

3. The Total Potential Energy functional On the basis of the kinematics (2)–(6) depicted in Figs. 1 and 2, the

pointwise displacement field components dx (x, y) and dy (x, y) along the x

9

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

and y directions, respectively, read:   c tl   u (x) − φ (x) y − − , if y ∈ Al  l l   2 2       
 1   ucu (x) + ucl (x) 2 dx (x, y) = (7)  4y 2    −φ (x) y + p(x) 1 − , if y ∈ A  c c  c2             u (x) − φ (x) y + c + tu , if y ∈ A u u u 2 2 where Al , Ac , and Au are the cross-section regions occupied by the lower skin,

TE D

the core, and the upper skin, 2 and  c c   vl (x) ≡ v(x) + sC (x) + sL (x) , if y ∈ Al   2 4        y2 dy (x, y) = (8) v(x) + sC (x) y + sL (x) , if y ∈ Ac  c           vu (x) ≡ v(x) − sC (x) c + sL (x) c , if y ∈ Au 2 4 Then, the strain field assumes the following form (henceforth, the symbol f 0

AC C

EP

denotes the derivative of the function f with respect to x):   c tl  0 0  u (x) − φ (x) y − − , if y ∈ Al  l l   2 2       
 1 0   ucu (x) + ucl (x)0 2 εx (x, y) =  4y 2   0 0  −φ (x) y + p (x) 1 − , if y ∈ Ac  c  c2             u0 (x) − φ0 (x) y + c + tu , if y ∈ A u u u 2 2 2

Henceforth, Al , Ac , and Au will be employed for the areas of those regions as well.

10

ACCEPTED MANUSCRIPT

RI PT

εy (x, y) =

 2y   sC (x) + sL (x) , if y ∈ Ac   c

M AN U

SC

    0 , if y ∈ A ∪ A l u    vl0 (x) − φl (x) , if y ∈ Al          2    v 0 (x) − φc (x) + s0 (x)y + s0 (x) y C L c γxy (x, y) = 8y   −p(x) 2 , if y ∈ Ac   c           v 0 (x) − φ (x) , if y ∈ A u u u

(9)

By considering possible core orthotropy, the relevant constitutive relations read

3

    E ε (x, y) , if y ∈ A , with i = l, u i x i    E (xy) εx (x, y) + Ec(y) εy (x, y)] , if y ∈ Ac   c

TE D

σx (x, y) =

   E (x) ε (x, y) + Ec(xy) εy (x, y) , if y ∈ Ac   c x

σy (x, y) =

(10)

(11)

EP

    0 , if y ∈ (A ∪ A ) l u

τxy (x, y) = Gi γxy (x, y) , if y ∈ Ai , with i = l, c, u

(12)

AC C

where the choice

Ec(x) = Ec(y) =

3

Ec , 1 − νc2

Ec(xy) =

E c νc 1 − νc2

(13)

The modification of the constitutive relations (10) and (11) to model the cylindrical

bending of sandwich plates is straightforward, though not considered here.

11

ACCEPTED MANUSCRIPT

refers to plane stress isotropic core behaviour of Young modulus Ec , Poisson’s (x)

the choice Ec

RI PT

ratio νc , and shear modulus Gc = Ec /[2(1 + νc )]. In the limit εy → 0,

= Ec allows the particularisation of the present model to

that of Krajcinovic (1972). Under the further constraint φu = φc = φl , the Timoshenko beam model is recovered.

SC

The foregoing expressions allow the determination of the elastic strain

energy, so that the Total Potential Energy (TPE) functional turns out to be:

2

M AN U


Ψ ul (x), uu (x), φl (x), φu (x), v(x), vl (x), vu (x), p(x) Z ( 1 L 2 2 = Eu Iu (φ0u ) + Eu Au (u0u ) 2 0 2

2

2

AC C

EP

TE D

+ Gu Au (vu0 − φu ) + El Il (φ0l ) + El Al (u0l ) + Gl Al (vl0 − φl ) n 1 h 2  2   i (x) 0 0 tl 0 0 tu 0 0 tu 0 0 tl + Ec Ac ul + φl uu − φ u + uu − φu + ul + φl 3 2 2 2 2  o 8 0 2 2 0 0 0 tu 0 0 tl (p ) + p uu − φu + ul + φl 15 3 2 2 (y) h  
i Ec Ac 2 2 + 7 vl + vu − 2vl vu + 16v v − vu − vl 3c2  h 0 0 tu vl − vu (xy) 0 0 tl + Ec Ac ul + φl + uu − φu 2 2 c   i
2 0 4 0 0 tu 0 0 tl − u − φu − ul − φl vl + vu − 2v + p (vl − vu ) 3c u 2 2 3c h  i2  1 t t 1 8 2 l u + Gc Ac v 0 − vl0 − vu0 − p uu − φu − ul − φl + c 2 2 12 c   i h
tu tl 1 uu − φu − ul − φl vl0 + vu0 − 2v 0 + 4 v0 − c 2 2 )
1 0 2 + v + vu0 − 2v 0 dx − V (14) 20 l

in which L is the beam length, Ii = b t3i /12 (with i = l, u and b denoting the cross-section width), and V accounts for the work of the applied loads. For 12

ACCEPTED MANUSCRIPT

the sake of brevity, here and henceforth, we drop the dependence on x of the

the minimum of the TPE functional is attained: min (ul ,uu ,φl ,φu ,v,vl ,vu ,p)∈K

RI PT

structural unknown functions. The equilibrium conditions are satisfied when

Ψ(ul , uu , φl , φu , v, vl , vu , p)

(15)

SC

with K the set of all the compatible fields (ul , uu , φl , φu , v, vl , vu , p), which

M AN U

must be continuous and must respect the essential boundary conditions. 4. The shear stress evaluation

The shear stress field obtained by substituting Eqs. (9) in Eq. (12) can be largely improved by following Bardella and Tonelli (2012), who have proposed to take advantage of Jourawski’s approach (Jourawski, 1856), that is classically applied to Bernoulli-Navier beams. As is well-known, this “stress

TE D

recovery” approach consists of imposing the force balance along the x direction of any region of length dx, obtained by cutting such an infinitesimal slice along any plane orthogonal to the y axis. From such a procedure, we obtain

EP

the average over b of τxy (x, y, z), that is the Jourawski-like estimate τ zz xy (x, y).

AC C

The analytic result of this procedure for the structural model developed in

13

ACCEPTED MANUSCRIPT

RI PT

this work leads to the following expressions:  1   El (c + 2 tl − 2 y)[c φ00l + 4 u00l − 2 φ00l y] , if y ∈ Al   8           c − 2 y n (xy) h 0 y i 0  00  + E t u E 4 s + 2 s 1 + 2 l l l  c C L  8 c  h (x) zz 00 00 00 00 +Ec φc c + 2 φl tl + 4 ul − 2 φc y τ xy (x, y) =    y y 2 io 8 00    1 − p − 2 , if y ∈ Ac +    3 c c2         1   Eu (c + 2 tu + 2 y)[c φ00u − 4 u00u + 2 φ00u y] , if y ∈ Au 8

M AN U

SC

(16)

In the following, we employ the terminology introduced in our previous works (see, e.g., Mattei and Bardella, 2016 and references therein) and indicate relations (16) as the Jourawski-zigzag estimate.

TE D

5. The new sandwich beam finite element In this work, in order to demonstrate the modelling capability of the proposed structural theory, to discuss it, and to provide a computationally

EP

efficient tool, we implement the new theory in a one-dimensional finite element. Hence, for the sake of brevity, we skip the derivation and discussion

AC C

of the (Euler–Lagrange) balance equations and natural boundary conditions ensuing from the minimum (15). They would be similar to (though, slightly richer than) those analysed by Mattei and Bardella (2016) for the 6DOFs model obtainable from the present model by constraining sL (x) and p(x) to

vanish. In the following, we limit our attention to writing the weak form of the balance equations (Sect. 5.1) and, then, the FE algorithm (Sect. 5.2 and Appendix A). 14

ACCEPTED MANUSCRIPT

5.1. The weak form of the balance equations

RI PT

By neglecting concentrated loads for the sake of brevity, the stationarity of the TPE functional (14) requires the following (balance) equations to be satisfied: Z L ("

(x)

  1 0 0 0 0 0 2uu − φu tu + ul + φl tl + 2p 2 0 #   (xy) Ec Ac vl − vu 2 vl + vu − 2 v − + δu0u 2 c 3 c   uu − 12 φu tu − ul − 12 φl tl Gc Ac 0 − v − c c )  Z L 1 0 0 0 quu δuu dx = 0 + (vl + vu − 2 v ) δuu dx − 6 0

SC

Ec Ac + 6

(17)

M AN U

Eu Au u0u

TE D

with quu (x) the distributed axial load applied on the upper skin.

Z L (

EP

 Gc Ac h 1 0 5 − vl + 2 vu0 + v 0 + (φu tu 15 2 4c 0  i 0 +2 ul + φl tl − 2 uu + 8p) + Gu Au (vu − φu ) δvu0    5 1 Ac − 2 Ec(xy) − φ0u tu + φ0l tl + u0l + 5u0u + 4p0 c 6c 2 2  ) Z L

AC C

+ 2Ec(y) (8 v − vl − 7 vu ) δvu

dx −

qvu δvu dx = 0 0

with qvu (x) the distributed transverse load applied on the upper skin.

15

(18)

ACCEPTED MANUSCRIPT

RI PT

Z L (  tu Ac h (x)  0 1 − Ec 2 uu + u0l + φ0l tl − φ0u tu + 2p0 12 2 0  (xy) i Ec 0 ( vl + 4v − 5 vu ) + Eu Iu φu δφ0u + c  tu Gc Ac  + 3 φu tu + (vu0 + vl0 + 4v 0 )c + 6 ul 12 c2  ) Z L  0 − 6 uu + 3φl tl − Gu Au (vu − φu ) δφu dx − qφu δφu dx = 0

SC

(19)

M AN U

0

with qφu (x) the distributed couple applied on the upper skin. Z L (

TE D

Ac h (x) 0 Ec (4ul + 2 φ0l tl + 2 u0u − φ0u tu + 4p0 ) 12 0  (xy) i 2Ec 0 + (5 vl − vu − 4 v) + El Al ul δu0l c Gc Ac h (4 v 0 + vl0 + vu0 )c + 3 φl tl − 6 uu + 6c2 ) Z L i + 3 φu tu + 6 ul δul dx − qul δul dx = 0

(20)

0

EP

with qul (x) the distributed axial load applied on the lower skin. Z L (

AC C

 5 Gc Ac h 0 1 0 2vl − vu + v 0 + (φu tu + 2ul − 8p 15 2 4c 0  i + φl tl − 2 uu ) + Gl Al (vl0 − φl ) δvl0  1  Ac h 5 + 2 Ec(xy) − φ0u tu + φ0l tl + 5u0l + u0u + 4p0 c 6c 2 2 ) Z L i (y) + 2Ec (7 vl + vu − 8 v) δvl dx − qvl δvl dx = 0 0

16

(21)

ACCEPTED MANUSCRIPT

RI PT

with qvl (x) the distributed transverse load applied on the lower skin. Z L (

(22)

SC

 1 tl Ac h (x)  1 0 1 Ec φl tl − φ0u tu + u0l + u0u + p0 6 2 4 2 0  (xy) i Ec  5 1 0 + vl − vu − 2 v + El Il φl δφ0l c 2 2  Gc tl Ac  + 3 φu tu + (vu0 + vl0 + 4 v 0 ) + 6 ul 12 c2  ) Z L  0 − 6 uu + 3 φl tl − Gl Al (vl − φl ) δφl dx − qφl δφl dx = 0

M AN U

0

with qφl (x) the distributed couple applied on the lower skin. Z L(

(23)

TE D

Ac Gc  (8 v 0 + vl0 + vu0 )c + 5 φl tl − 10 uu 15c 0 1  2Ac h + 5 φu tu + 10 ul δv 0 − 2 Ec(xy) φ0u tu − u0u 3c 2 ) Z L  i 1 0 0 (y) + ul + φl tl c + 4Ec ( vl + vu − 2v) δv dx − qv δv dx = 0 2 0

AC C

EP

with qv (x) the distributed transverse load applied on the core centre (x axis). #  Z L (" (x)  Ec 8 0 1 v − v 1 2 l u p + u0l + φ0l tl + u0u − φ0u tu + Ec(xy) Ac δp0 3 5 2 2 3 c 0 (24)  0  ) Z L 2 vl − vu0 8p − Gc Ac − 2 δp dx − qp δp dx = 0 3 c c 0 with qp (x) a distributed axial load applied along the beam axis x and working for p.

5.2. The Finite Element discretisation By denoting with ξ(r) a generic spatial function written in terms of the intrinsic coordinate r ∈ [1, −1] on the parent element, we assume the following 17

ACCEPTED MANUSCRIPT

isoparametric FE interpolation: X

ξ (i) N (i) (r)

i

(25)

RI PT

ξ(r) =

in which ξ (i) (with i = 1, 2, 3, · · · ) are the nodal values of ξ and N (i) (r) are

nodes.

SC

the interpolation functions, assuming value 1 at the node i and 0 at the other

Let nG be the total number of Gauss points, each one characterised by a (i)

M AN U

position rn and a weight ωn . Let Nn be the shape function referred to the i-th node, evaluated at n-th Gauss point, i.e., at r = rn . It results: 1 dN (i) dN (i) (r(x)) 0(i) = Nn = dx Jn dr r=rn r=rn where

(i) X dx(r) (i) dN Jn = = x dr r=rn dr r=rn i

(26)

(27)

TE D

is the jacobian operator, x(i) being the nodal coordinate of the i-th node. The above definitions allow us to write the FE stiffness matrix as reported in Appendix A. This is basically all is needed for a linear finite element, in

EP

this investigation implemented in a user subroutine (uel) for the commercial code ABAQUS (Dassault Syst`emes, 2013).

AC C

6. Results and discussion 6.1. The FE models: discretisations and computational times All the FE analyses are performed with the commercial code ABAQUS

(Dassault Syst`emes, 2013). The continuum plane stress FE models, taken as reference for all the benchmarks discussed in this work, are totally similar to those described in 18

ACCEPTED MANUSCRIPT

detail in Tonelli et al. (2012) and Mattei and Bardella (2016). They consist

RI PT

of 140800 fully integrated, linear plane stress elements, having 142489 nodes and 284978 DOFs overall. Each analysis requires ≈ 25 s of computational

time on a workstation endowed with 24 GB RAM and CPU characterised by a 3.2 GHz clock and 6 cores.

SC

In comparison, the FE structural models consist of 200 cubic (4-noded)

finite elements of equal size, for 601 nodes overall. Such FE analyses requires, 4

M AN U

on the same workstation, only less than 0.3 s of computational time.

6.2. Simply-supported beam subject to uniform load under SBCs: A critical assessment of some literature models

Here we refer to the benchmark of Tonelli et al. (2012), illustrated in Fig. 3. Through this benchmark we seek to demonstrate that the proposed

TE D

8DOFs model can accurately predict the stress state under SBCs in most of the sandwich domain for any relevant relative stiffness between core and skins. As shown in the following, this is not the case of some other literature models, whose predictions are accurate in a limited range of relative stiffness

EP

only.

As in Tonelli et al. (2012), we consider a sandwich beam having identical

AC C

skins, that is tu ≡ tl ≡ t (so that h = c + 2t) and Eu ≡ El ≡ Es (with the subscript s denoting “skins”). Each layer has Poisson’s ratio equal to 4

Note that a precise comparison in terms of computational cost between the contin-

uum and structural FE analyses could be performed only after optimising the FE meshes employed in the continuum analyses, which would turn out to be non-uniform and problemdependent. Such a study is beyond the scope of the present investigation.

19

ACCEPTED MANUSCRIPT

0.3. The beam slenderness is h/L = 3/20, with h denoting the cross-section

RI PT

height (here, h = c + 2t).

The two-dimensional boundary value problem consists of a simply-supported sandwich beam subject to a uniformly distributed load applied at the

top skin (y = −c/2 − t) and supports realised at the bottom corner nodes

SC

(y = c/2 + t). Because of the symmetry, we model half of the problem and

L refers to the length between the support and the mid cross-section, where

M AN U

the shear force vanishes. We provide the results in terms of the distance d from the support (see Fig. 3) limiting the region x ∈ [d, L] in which the considered theoretical model predicts a stress field almost coincident with that predicted by the two-dimensional plane stress FE analysis. To check this we follow Tonelli et al. (2012) and consider, in the plane stress analysis, the longitudinal stress σx in the Gauss points closest the outer skins’ boundaries:

TE D

we assume a tolerance of 1% on the relative error of the theoretical model with respect to this quantity.

The main goal of Tonelli et al. (2012) is to provide information on the accuracy of various theoretical predictions as a function of the relative stiffness

EP

between core and skins, the latter being somehow proportional to the sandwich heterogeneity H, which should be a function, among other parameters,

AC C

of Ec /Es and t/c. The sandwich response depends in a complicated way on

heterogeneity H, the slenderness h/L, and the boundary conditions. We consider the following models:

• the classical FOSD theory (Timoshenko, 1921; Bert, 1973; Gordaninejad and Bert, 1989; Bardella, 2008), which just requires, for what concerns the stress evaluation in a statically determinate case as that here 20

ACCEPTED MANUSCRIPT

dealt with, the straightforward cross-sectional analysis inherent to the

RI PT

Bernoulli-Navier kinematics;

• the Krajcinovic (1972) zig-zag theory, with Jourawsky-zigzag estimate for the shear stress proposed by Bardella and Tonelli (2012);

SC

• the Frostig et al. (1992) theory for antiplane sandwiches having thin skins;

and 6DOFs models.

M AN U

• our new 8DOFs model, considering also its simplifications to the 7DOFs

The results for the first three models are taken from Tonelli et al. (2012). The obtained predictions are reported in Figs. 4 and 5. Let us remind that the 7DOFs model is obtained from the 8DOFs model by constraining

TE D

p(x) to be zero, thus neglecting the quadratic term in the core longitudinal displacement, and the 6DOFs model, that is the Mattei and Bardella (2016) model, is obtained by setting sL (x) = p(x) = 0, thus also neglecting the linear variation of the transverse core strain. Among the tested models, only

strain.

EP

the FOSD and Krajcinovic (1972) theories totally neglect the transverse core

AC C

It is clear that the 8DOFs model is the only one providing good predictions for every sandwich heterogeneity H. In order to reach this goal, it is important to consider the linear contribution sL (x) to the transverse core strain, neglected in Mattei and Bardella (2016). This contribution plays a role not only, as expected, for soft cores (see Fig. 4, concerned with thick skins), but also when the sandwich is close to a homogeneous beam, as clearly seen in Fig. 5, concerned with thin skins. In fact, when Ec /Es is close to 21

ACCEPTED MANUSCRIPT

1, the core modulus employed in our model, that is Ec /(1 − νc )2 (see Eqs.

RI PT

(10)–(13)), requires an appropriate compliance along the y direction to compensate for the fact that it is larger than Ec . Such a compliance is provided by

allowing a linear variation of εy , thus conciliating between two-dimensional plane stress and classical beam theory.

SC

Also accounting for p(x) seems to be worth, as the predictions of the 8DOFs model are significantly better than those of the 7DOFs model. This

M AN U

will be confirmed by the benchmarks of sections 6.3 and 6.4. The 8DOFs predictions are the best ones for almost every H, with two exceptions. The first exception is related to very soft cores, for which the Frostig et al. (1992) model provides slightly better predictions for certain mid values of t/c. On the one hand, this is expected as the Frostig et al. (1992) model describes the core as a plane stress continuum (though neglecting σx and,

TE D

then, predicting τxy independent of y, because of the antiplane assumption). On the other hand, this is quite interesting because the Frostig et al. (1992) model is developed under the hypothesis of thin skins, while here it turns out to provide the best predictions for not-too-thin faces. Most likely, this

EP

is because those mid values of t/c are such that the difference in stiffness between core and skins is large enough to make the antiplane core assumption

AC C

accurate, while t/c is still small enough to make the hypothesis of thin skins acceptable.

The second exception is observed in Fig. 5 for t/c = 1/28 and 1/118. Such

exception is explained by analysing the stress state, as for instance illustrated in Fig. 6 in terms of the modulus of the relative error on the longitudinal stress, for the case t/c = 1/118 and Ec /Es = 0.9. First, let us emphasise

22

ACCEPTED MANUSCRIPT

that, overall, the 8DOFs model must always be more accurate than the

RI PT

7DOFs model, in the sense that the TPE minimum obtained when adding the eighth DOF (p(x)) to the 7DOFs model may not increase. However,

pointwise, the improvement of the 8DOFs model turns out to be mostly concerned with the stress state at the lower skin (void circles in Fig. 6),

SC

where the stress is larger. Instead, in the upper skin, the errors of the 8DOFs and 7DOFs models have opposite sign, with the modulus of the error of the

M AN U

latter going below the 1% tolerance at a slightly shorter distance, from the support, with respect to that of the 8DOFs model. Finally, there is also a less relevant issue concerned with setting the tolerance: from Fig. 6 it is clear that slightly diminishing the tolerance, in this case, would make the distance d for the 8DOFs model shorter than that for the 7DOFs model. However, the important outcome of this analysis is that in a region quite close to the

TE D

support, that is the one mostly affected by SBCs, the 8DOFs model seems to provide significantly better stress state predictions; since in this region all the models give an error larger than the tolerance, this information cannot

EP

be included in graphics such as those reported in Figs. 4 and 5. 6.3. Propped-cantilever beam subject to uniform load under SBCs

AC C

As discussed by Bardella and Mattei (2014) and Mattei and Bardella (2016), the propped-cantilever benchmark is particularly interesting for it is a statically indeterminate case and because of the fully clamped beam end condition. Such condition, imposing all the displacement components to be zero pointwise at the beam end, in fact, leads to difficulties in the prediction of the stress state close to the clamped section in the two-dimensional FE analysis, due to the unboundedness of the longitudinal stress σx at the 23

ACCEPTED MANUSCRIPT

well, 1982; Tullini and Savoia, 1995).

RI PT

clamped section corners in the linear elastic solution (see Gregory and Glad-

In this benchmark, we follow Mattei and Bardella (2016) and apply the support condition in such a way as to set to zero the transverse displacement

component of each point of the lower skin, while the uniform load is applied

SC

on the upper skin. We choose mechanical properties of the layers such that

Eu /El = 1.6 and, unless otherwise specified, Eu /Ec = 166.6. The shear

M AN U

moduli are obtained by assuming, for all the layers, isotropy and Poisson’s ratio equal to 0.3. The geometry is characterised by tl = c/2, tu = c/4, and h/L =0.175, with h = tl + c + tu .

For the sake of brevity, the results are provided only in terms of the shear stress prediction, obtained, in the structural models, from the longitudinal stress through the Jourawski-zigzag estimate (16). The results are displayed

TE D

in Figs. 7–13. We first point out that the results related to the Mattei and Bardella (2016) model are here re-computed by using our FE implementation (constrained to obtain the 6DOFs model). Such results are in fact slightly more accurate than those obtained by Mattei and Bardella (2016), who em-

EP

ployed the Rayleigh-Ritz method to obtain a numerical solution (the most evident case is that for Eu /Ec = 16 in Fig. 13, to be compared with figure

AC C

20 of Mattei and Bardella, 2016). Note that the results obtained for the 7DOFs model are plotted only in the cross-sections where their difference with respect to the results of the 8DOFs model is relevant. Close to the clamped section (Figs. 7 and 8) the comparison is difficult

and shows the well-known issues recently highlighted by Bardella and Mattei (2014) and Mattei and Bardella (2016): the continuum FE solution is

24

ACCEPTED MANUSCRIPT

unreliable (as it predicts non-vanishing shear stress at the top and bottom

RI PT

boundaries), while the Jourawski-zigzag estimate may be useful. As shown

in Fig. 9, in the fourth of beam on the clamped section side, x ∈ [0, L/4],

overall, there is not appreciable improvement in the predictions of the 8DOFs model with respect to the 6DOFs model. However, in the beam region be-

SC

tween the mid-section to the simply-supported section (Figs. 10, 11, and 12)

the superior modelling capability of the 8DOFs model clearly emerges. In

M AN U

this region, the quality of the shear stress prediction strongly depends on the accuracy in evaluating the support reaction force. In fact, in this regard, the percentage relative error with respect to the two-dimensional continuum model is ≈ 0.037% for the 8DOFs model, ≈ 0.13% for the 7DOFs model, and ≈ 0.72% for the 6DOFs model. As is clear from Fig. 12, the 8DOFs model can partly mitigate the issue, physiological for a structural theory, of

TE D

being unable to predict zero shear stress pointwise in the core and in the upper skin of the propped section. In this section, by the way, the typical spurious stress values at the interface provided by the continuum FE analysis are clearly visible.

EP

Summing up, far enough from the beam ends, for about x ∈ [L/8, 7L/8], the match between the two-dimensional FE results and the Jourawski-zigzag

AC C

predictions ensuing from the 8DOFs model is excellent. The above results are further confirmed by Fig. 13, focussing on the section x = 7L/8 and showing the results for different values of Eu /Ec .

6.4. A novel benchmark with particularly severe boundary conditions In order to further demonstrate the accuracy of the proposed formulation and to highlight the different modelling capability of the 6DOFs, 7DOFs, 25

ACCEPTED MANUSCRIPT

and 8DOFs models, here we refer to a sandwich beam subject to peculiar

RI PT

boundary conditions, as depicted in Fig. 14, whose purpose is to promote highly inhomogeneous core deformability. Constraints are applied only to

the upper skin of one beam end (x = 0), that is clamped. This implies that uu (0) = vu (0) = φu (0) = 0 in the structural model, while the two

SC

displacement components of all the nodes belonging to that skin side are set to zero in the two-dimensional continuum FE analysis. The beam is

M AN U

loaded with both transverse and axial forces, as depicted in Fig. 14. In order to emphasise the SBCs, in the two-dimensional continuum FE model, such forces are applied on the central node of the skins’ sides. The material and geometrical data are the same as those for the proppedcantilever benchmark of section 6.3 (in particular, Eu /Ec = 166.6). Figs. 15–18 display the comparison, again in terms of shear stress, be-

TE D

tween the two-dimensional continuum model and the structural model, in its three versions implemented in this work: 6DOFs, 7DOFs, and 8DOFs. The 8DOFs model is the sole one able to almost perfectly match the twodimensional continuum FE results in all the considered cross-sections. This

EP

clearly emerges in the cross-section at x = L/4 (Fig. 16), in which the 8DOFs model can even reproduce a shear stress which changes second derivative with

AC C

respect to y along the core.

7. Concluding remarks We have extended the structural model for sandwich beams recently pro-

posed by Mattei and Bardella (2016) in such a way as to be able to describe the mechanical response of sandwich beams of any relative stiffness between 26

ACCEPTED MANUSCRIPT

core and skins, even when the sandwich is subject to Severe Boundary Condi-

RI PT

tions (SBCs). More specifically, the model aims at accurately predicting the

stress state of linear elastic plane sandwich beams of any cross-section, including thin or thick skins and soft or stiff cores (Allen, 1969). The model of Mattei and Bardella (2016) is based on the Yu–Krajcinovic zig-zag structural

SC

theory, in which all layers rotate independently of the beam axis (Yu, 1959; Krajcinovic, 1972), and accounts also for the transverse core strain, assumed

M AN U

to be independent of the through-the-thickness coordinate y by Mattei and Bardella (2016). Our extension, with respect to Mattei and Bardella (2016), consists of allowing both the transverse and longitudinal displacement components in the core to be quadratic in y, thus formulating the structural model in terms of eight independent functions, two more than those employed by Mattei and Bardella (2016).

TE D

Moreover, the shear stress predictions are improved, with respect to those directly obtainable by the structural model, by post-processing the longitudinal stress with a Jourawski-like technique. The new structural theory has been implemented in a one-dimensional

EP

finite element and its predictions have been compared both with those of other literature models (First Order Shear Deformation theory; Krajcinovic,

AC C

1972; Frostig et al., 1992; Mattei and Bardella, 2016) and with reference results obtained from accurate continuum plane stress FE analyses. In order to demonstrate the high accuracy and computational efficiency of

the new structural model, we have resorted to the following boundary value problems, all characterised by large stress gradients for the applied SBCs: (i) the simply-supported beam subject to uniform transverse load employed

27

ACCEPTED MANUSCRIPT

by Tonelli et al. (2012) to critically evaluate some literature models, (ii)

RI PT

the propped-cantilever beam subject to uniform transverse load analysed by

Mattei and Bardella (2016), and (iii) a novel benchmark in which the core deformation results to be highly inhomogeneous.

We have separately discussed the benefits of introducing the two new in-

SC

dependent structural functions, with respect those characterising the model

of Mattei and Bardella (2016), thus showing that allowing for a linear varia-

M AN U

tion of the transverse core strain is crucial to obtain a model able to predict the stress state for both soft and stiff cores. However, we have also shown that, in order to obtain very accurate stress state estimates in the case of soft cores, one has to further enrich the core kinematics. We have accomplished this task by describing the longitudinal displacement component in the core with a quadratic function of y. Future work might be concerned with the

TE D

optimisation of the core kinematics by considering different shape functions (see, e.g., Tessler et al., 2009; Vidal and Polit, 2009; Carrera et al., 2013; Icardi and Sola, 2016, and references therein.)

EP

8. Acknowledgements

Work done within a research project financed by the Italian Ministry

AC C

of Education, University, and Research (MIUR). The finite element code ABAQUS has been run at the Department of Civil, Environmental, Architectural Engineering and Mathematics (DICATAM), University of Brescia, Italy, under an academic license.

28

RI PT

ACCEPTED MANUSCRIPT

¥u(x)

SC

uu(x)

tu c 2

p(x)

2

tl

ul (x)

vu(x)

¥l (x) y

v(x)

¥u(x)

TE D

v l (x)

x

M AN U

c

x

¥l (x)

AC C

EP

p(x)

Figure 1: Enriched zig-zag kinematics for a sandwich beams with core deformable along the y direction and admitting quadratic displacement field: independent structural functions listed in Eq. (2).

29

RI PT

ACCEPTED MANUSCRIPT

tu

SC

u cu(x)

c 2

x

2

u cl (x) tl

M AN U

¥c(x)

c

v(x)- c s C (x)+ c sL(x) 2

y

4

TE D

¥c (x)

v(x)+ c s C (x)+ c sL(x) 4

x

AC C

EP

2

v'(x)

Figure 2: Enriched zig-zag kinematics for a sandwich beams with core deformable along the y direction and admitting quadratic displacement field: definition of structural functions dependent on the independent ones listed in Eq. (2).

30

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 3: Sketch of the benchmark with the notation employed to indicate the region x ∈ [d, L] where the models analysed accurately predict the stress state — from Tonelli

AC C

EP

et al. (2012).

31

6

5

5

4

4

3

3 2

1

d/h

0 0.0001

0.01

7

TE D

0 0.0001

0.01

1 Ec/Es

1 Ec/Es Frostig et al., 1992 FOSD Krajcinovic, 1972 6 DOFs 7 DOFs 8DOFs

7

5

EP

1

0.01

6

5

2

Ec/Es

d)

6

3

0 0.0001

1

d/h

1

4

RI PT

6

2

c)

7

SC

b)

d/h

7

M AN U

a)

d/h

ACCEPTED MANUSCRIPT

4 3 2 1 0 0.0001

0.01

1 Ec/Es

AC C

Figure 4: Extent of the region x ∈ [0, d] where considered models inaccurately predict the sandwich stress state. a) t/c = 1/3, b) t/c = 1/4, c) t/c = 1/6, d) t/c = 1/8.

32

6

5

5

4

4

3

3 2

1

d/h

0 0.0001

0.01

7

TE D

0 0.0001

0.01

1 Ec/Es

1 Ec/Es Frostig et al., 1992 FOSD Krajcinovic, 1972 6 DOFs 7 DOFs 8 DOFs

7

5

EP

1

0.01

6

5

2

Ec/Es

d)

6

3

0 0.0001

1

d/h

1

4

RI PT

6

2

c)

7

SC

b)

d/h

7

M AN U

a)

d/h

ACCEPTED MANUSCRIPT

4 3 2 1 0 0.0001

0.01

1 Ec/Es

AC C

Figure 5: Extent of the region x ∈ [0, d] where considered models inaccurately predict the sandwich stress state. a) t/c = 1/10, b) t/c = 1/13, c) t/c = 1/28, d) t/c = 1/118.

33

RI PT SC

0.015

0.01

7DOFs, bottom 7DOFs, top 8DOFs, bottom 8DOFs, top Tolerance

0

0.02

0.04

0.06

EP

0

TE D

0.005

M AN U

Modulus of the relative error

ACCEPTED MANUSCRIPT

0.08

0.1

0.12 x/L

Figure 6: Modulus of the relative error on the longitudinal stress, for the case t/c = 1/118

AC C

and Ec /Es = 0.9.

34

3

2.5

RI PT

Aτxy /(qL)

ACCEPTED MANUSCRIPT

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 8DOFs

2

1.5

0.5

0 0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

M AN U

0

SC

1

Figure 7: Nondimensional shear stress at the fully clamped section (x = 0) of the proppedcantilever beam. Comparison among continuum plane stress FE, structural model of

2

1.5

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 8DOFs

EP

1

TE D

Aτxy /(qL)

Mattei and Bardella (2016) (6DOFs), and the present structural model (8DOFs).

AC C

0.5

0

0

0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

Figure 8: Nondimensional shear stress at the section x = L/80 of the propped-cantilever beam. Comparison among continuum plane stress FE, structural model of Mattei and Bardella (2016) (6DOFs), and the present structural model (8DOFs).

35

x=L/8 1

RI PT

Aτxy /(qL)

ACCEPTED MANUSCRIPT

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 8DOFs

0.8

0.6

x=L/4

0.2

0 0.2

0.4

0.6

0.8

1 (tld + c/2 - y)/h

M AN U

0

SC

0.4

Figure 9: Nondimensional shear stress at the sections x = L/8 and x = L/4 of the propped-cantilever beam. Comparison among continuum plane stress FE, structural model

0.12

0.1

0.08

EP

0.06

TE D

Aτxy /(qL)

of Mattei and Bardella (2016) (6DOFs), and the present structural model (8DOFs).

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 8DOFs

0.04

AC C

0.02

0

0

0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

Figure 10: Nondimensional shear stress at the section x = L/2 of the propped-cantilever

beam. Comparison among continuum plane stress FE, structural model of Mattei and Bardella (2016) (6DOFs), and the present structural model (8DOFs).

36

-0.5 2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

-0.4

RI PT

Aτxy /(qL)

ACCEPTED MANUSCRIPT

x=7L/8

-0.3

-0.2

0 0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

M AN U

0

SC

x=3L/4

-0.1

Figure 11: Nondimensional shear stress at the sections x = 3L/4 and x = 7L/8 of the propped-cantilever beam. Comparison among continuum plane stress FE, structural model of Mattei and Bardella (2016) (6DOFs), and the present structural model (7DOFs and

1

0.5

0

TE D

Aτxy /(qL)

8DOFs).

0

0.2

0.4

EP

-0.5

0.6

0.8

1 (tl + c/2 - y)/h

-1

-1.5

AC C

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

-2

-2.5

Figure 12: Nondimensional shear stress at the propped section (x = L) of the propped-

cantilever beam. Comparison among continuum plane stress FE, structural model of Mattei and Bardella (2016) (6DOFs), and the present structural model (7DOFs and 8DOFs).

37

-0.9

2D Plane Stress FE, Eu /Ec =1600 2D Plane Stress FE, Eu /Ec =400

-0.8

2D Plane Stress FE, Eu /Ec =16 Jourawski-zigzag, 6DOFs, Eu /Ec =1600

-0.7

Jourawski-zigzag, 6DOFs, Eu /Ec =400 Jourawski-zigzag, 6DOFs, Eu /Ec =16

-0.6

Jourawski-zigzag, 8DOFs, Eu /Ec =1600 Jourawski-zigzag, 8DOFs, Eu /Ec =400

-0.5

Jourawski-zigzag, 8DOFs, Eu /Ec =16

SC

-0.4

RI PT

Aτxy /(qL)

ACCEPTED MANUSCRIPT

-0.3 -0.2

0 0

0.2

M AN U

-0.1

0.4

0.6

0.8

1 (tl + c/2 - y)/h

Figure 13: Nondimensional shear stress at the section x = 7L/8 of the propped-cantilever beam as a function of the core compliance. Comparison among continuum plane stress

AC C

EP

model (8DOFs).

TE D

FE, structural model of Mattei and Bardella (2016) (6DOFs), and the present structural

Figure 14: Schematic of the benchmark employed to test the new structural theory against highly inhomogeneous core deformation.

38

SC

40

M AN U

Aτxy /F

RI PT

ACCEPTED MANUSCRIPT

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

30

20

10

0

0.2

0.4

0.6

TE D

0

0.8

1

(tl + c/2 - y)/h

Figure 15: Nondimensional shear stress at the section x = L/8 in the benchmark of Fig. 14. Comparison among continuum plane stress FE, structural model of Mattei and

AC C

EP

Bardella (2016) (6DOFs), the present structural model (7DOFs and 8DOFs).

39

Aτxy /F

ACCEPTED MANUSCRIPT

4

RI PT

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

3

1

0 0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

M AN U

0

SC

2

Figure 16: Nondimensional shear stress at the section x = L/4 in the benchmark of Fig. 14. Comparison among continuum plane stress FE, structural model of Mattei and

TE D

Bardella (2016) (6DOFs), the present structural model (7DOFs and 8DOFs).

Appendix A. The FE stiffness matrix In writing the FE stiffness matrix, we directly use the nodal DOFs as the

AC C

EP

indices to denote the matrix components. " ! # nG (x) X Ec Ac G A c c Ku(i) Eu Au + Nn0(i) Nn0(j) + Nn(i) Nn(j) Jn ωn (j) = 2 u uu 3 c n=1

# (xy) 1 Gc Ac (i) 0(j) 5 Ec Ac 0(i) (j) Kuu(i) vu(j) = − Nn Nn − Nn Nn Jn ωn 6 c 6 c n=1 " # nG (x) X Ec G c = − Nn0(i) Nn0(j) + 2 Nn(i) Nn(j) Ac tu Jn ωn Kuu(i) φ(j) u 6 2c n=1 nG X

"

40

2

x=L/2

x=3L/4

1.5

1

0 0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

M AN U

0

SC

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

0.5

RI PT

Aτxy /F

ACCEPTED MANUSCRIPT

Figure 17: Nondimensional shear stress at the sections x = L/2 and x = 3L/4 in the benchmark of Fig. 14. Comparison among continuum plane stress FE, structural model

TE D

Aτxy /F

of Mattei and Bardella (2016) (6DOFs), the present structural model (7DOFs and 8DOFs).

3

2D plane stress FE Jourawski-zigzag, 6DOFs Jourawski-zigzag, 7DOFs Jourawski-zigzag, 8DOFs

EP

2

AC C

1

0

0

0.2

0.4

0.6

0.8

1 (tl + c/2 - y)/h

Figure 18: Nondimensional shear stress at the section x = 7L/8 in the benchmark of

Fig. 14. Comparison among continuum plane stress FE, structural model of Mattei and Bardella (2016) (6DOFs), the present structural model (7DOFs and 8DOFs).

41

ACCEPTED MANUSCRIPT

Kuu(i) v(j) =

nG X 

SC

RI PT

" # nG (x) X Ec G c Ku(i) Nn0(i) Nn0(j) − 2 Nn(i) Nn(j) Ac Jn ωn (j) = u ul 6 c n=1 " # n (xy) G X 1 Gc Ac (i) 0(j) 1 Ec Ac 0(i) (j) Nn Nn + Nn Nn Jn ωn Kuu(i) v(j) = − l 6 c 6 c n=1 " # nG (x) X G Ec c Kuu(i) φ(j) = Nn0(i) Nn0(j) − 2 Nn(i) Nn(j) Ac tl Jn ωn l 12 2c n=1 −Gc Nn(i) Nn0(j) + Ec(xy) Nn0(i) Nn(j)

n=1

M AN U

" # nG (x) X Ec = Nn0(i) Nn0(j) Ac Jn ωn 3 n=1

Ku(i) (j) u p

Kvu(i) φ(j) = u

#  nG  (y) X 7 E A 2 c c Gc Ac + Gu Au Nn0(i) Nn0(j) + Nn(i) Nn(j) Jn ωn = 2 15 3 c n=1

nG X n=1

"

(xy)

5 Ec Ac tu (i) 0(j) Nn Nn + 12 c

TE D

Kvu(i) vu(j)

 2Ac Jn ωn 3c

nG X 

Kvu(i) u(j) = l

l

  1 Gc tu Ac 0(i) (j) − Gu Au Nn Nn Jn ωn 12 c

−Ec(xy) Nn(i) Nn0(j) + Gc Nn0(i) Nn(j)

n=1

EP

Kvu(i) v(j) =



nG X

"

Gc Ac 0(i) 0(j) Ec Ac (i) (j) − N N + Nn Nn Jn ωn 30 n n 3c2

AC C

n=1

(y)

 Ac Jn ωn 6c #

nG X   Ac tl Jn ωn −Ec(xy) Nn(i) Nn0(j) + Gc Nn0(i) Nn(j) 12c n=1 # " n (y) G X Ac Gc 0(i) 0(j) 8Ec Ac (i) (j) = Nn Nn − Nn Nn Jn ωn 2 15 3c n=1

Kvu(i) φ(j) = l

Kvu(i) v(j)

Kvu(i) p(j) =

nG X 

−Ec(xy) Nn(i) Nn0(j) + Gc Nn0(i) Nn(j)

n=1

42

 2Ac Jn ωn 3c

ACCEPTED MANUSCRIPT

RI PT

 nG  X 1 (x) 2 Kφu(i) φu(j) = Ec Ac tu + Eu Iu Nn0(i) Nn0(j) 12 n=1    1 Gc Ac t2u (i) (j) + + Gu Au Nn Nn Jn ωn 4 c2 # " nG (x) X G Ec c Kφ(i) Ac tu − Nn0(i) Nn0(j) + 2 Nn(i) Nn(j) Jn ωn (j) = u ul 12 2c n=1

M AN U

SC

nG X  A c tu  Kφu(i) v(j) = Gc Nn(i) Nn0(j) − Ec(xy) Nn0(i) Nn(j) Jn ωn l 12c n=1 " # nG (x) X Ec G c Ac tu tl − Kφu(i) φ(j) = Nn0(i) Nn0(j) + 2 Nn(i) Nn(j) Jn ωn l 24 4c n=1 nG X  A c tu  = Gc Nn(i) Nn0(j) − Ec(xy) Nn0(i) Nn(j) Jn ωn 3c n=1

Kφu(i) v(j)

l

l

TE D

Ku(i) u(j)

nG X  Ac tu  −Ec(x) Nn0(i) Nn0(j) Jn ωn Kφu(i) p(j) = 6 n=1   nG  X 1 Gc Ac (i) (j) (x) 0(i) 0(j) Ac Ec + El Al Nn Nn + 2 Nn Nn Jn ωn = 3 c n=1

EP

nG X  Ac  Gc Nn(i) Nn0(j) + 5Ec(xy) Nn0(i) Nn(j) Jn ωn Ku(i) v(j) = l l 6c n=1 " # nG (x) X Ec G c Ku(i) φ(j) = Ac tl Nn0(i) Nn0(j) + 2 Nn(i) Nn(j) Jn ωn l l 6 2c n=1

AC C

nG X  2Ac  Gc Nn(i) Nn0(j) − Ec(xy) Nn0(i) Nn(j) Jn ωn Ku(i) v(j) = l 3c n=1 # " nG (x) X Ec Ku(i) p(j) = Ac Nn0(i) Nn0(j) Jn ωn l 3 n=1 #   nG (y) X 2 7Ec Ac (i) (j) 0(i) 0(j) Kv(i) v(j) = Gc Ac + Gl Al Nn Nn + Nn Nn Jn ωn l l 15 3c2 n=1

43

ACCEPTED MANUSCRIPT

Kv(i) φ(j) l

l

Kv(i) v(j) = l

nG X

"

n=1

RI PT

"    nG (xy) X 5Ec Ac tl (i) 0(j) 1 Gc Ac tl 0(i) (j) = Nn Nn + − Gl Al Nn Nn Jn ωn 12c 12 c n=1 # (y) 1 8E A c c Gc Ac Nn0(i) Nn0(j) − Nn(i) Nn(j) Jn ωn 15 3c2

nG X  2Ac  (xy) (i) 0(j) Ec Nn Nn − Gc Nn0(i) Nn(j) Jn ωn 3c n=1  nG  X 1 (x) 2 Ec Ac tl + El Il Nn0(i) Nn0(j) = 12 n=1    1 Gc Ac t2l (i) (j) + Gl Al Nn Nn Jn ωn + 4 c2

Kφ(i) φ(j) l

l

Kφ(i) v(j) l

nG X  Ac tl  = Gc Nn(i) Nn0(j) − Ec(xy) Nn0(i) Nn(j) Jn ωn 3c n=1

Kφ(i) p(j)

nG X Ac tl  (xy) 0(i) 0(j)  Ec Nn Nn Jn ωn = 6 n=1

TE D

l

Kv(i) v(j) =

M AN U

l

SC

Kv(i) p(j) =

nG X n=1

"

# (y) 8 16E A c c Ac Gc Nn0(i) Nn0(j) + Nn(i) Nn(j) Jn ωn 15 3c2

EP

Kv(i) p(j) = 0

AC C

Kp(i) p(j)

 nG  X 8 16Gc Ac (i) (j) (xx) 0(i) 0(j) = Ac Ec Nn Nn + Nn Nn Jn ωn 15 3c2 n=1

References

Allen, H. G., 1969. Analysis and design of structural sandwich panels. Pergamon Press Ltd., Oxford.

Bardella, L., 2008. Reliability of first-order shear deformation models for sandwich beams. J. Mech. Mater. Struct. 3 (7), 1187–1206. 44

ACCEPTED MANUSCRIPT

Bardella, L., Mattei, O., 2014. On explicit analytic solutions for the accurate

RI PT

evaluation of the shear stress in sandwich beams with a clamped end.

Compos. Struct. 112, 157–168, corrigendum: Compos. Struct., 116:849, 2014.

Bardella, L., Paterlini, L., Leronni, A., 2014. Accurate modelling of the linear

SC

elastic flexure of composite beams warped by midlayer slip, with emphasis

M AN U

on concrete-timber systems. Int. J. Mech. Sci. 87, 268–280.

Bardella, L., Tonelli, D., 2012. Explicit analytic solutions for the accurate evaluation of the shear stresses in sandwich beams. J. Eng. Mech. 138 (5), 502–507, erratum: J. Eng. Mech., 138(10):1302, 2012.

Berdichevsky, V. L., 2010. An asymptotic theory of sandwich plates. Int. J.

TE D

Eng. Sci. 48, 383–404.

Bert, C. W., 1973. Simplified analysis of static shear factors for beams of nonhomogeneous cross section. J. Compos. Mater. 7, 525–529.

EP

Carrera, E., Filippi, M., Zappino, E., 2013. Laminated beam analysis by polynomial, trigonometric, exponential and zig-zag theories. Eur. J. Mech.

AC C

A-Solid. 41, 58–69.

Dassault Syst`emes, 2013. ABAQUS User’s & Theory Manuals — Release 6.13-1. Providence, RI, USA.

Frostig, Y., Baruch, M., Vilnay, O., Shelnman, I., 1992. High-order theory for sandwich-beam behavior with transversely flexible core. J. Eng. Mech. 118 (5), 1026–1043. 45

ACCEPTED MANUSCRIPT

Galuppi, L., Royer-Carfagni, G. F., 2012. Effective thickness of laminated

RI PT

glass beams: New expression via a variational approach. Eng. Struct. 38, 53–67.

Ghugal, Y. M., Shimpi, R. P., 2001. A review of refined shear deformation

theories for isotropic and anisotropic laminated beams. J. Reinf. Plast.

SC

Comp. 20 (3), 255–272.

M AN U

Gordaninejad, F., Bert, C. W., 1989. A new theory for bending of thick sandwich beams. Int. J. Mech. Sci. 31 (11/12), 925–934. Gregory, R. D., Gladwell, I., 1982. The cantilever beam under tension, bending or flexure at infinity. J. Elast. 12 (4), 317–343.

Heller, R. A., 1969. Interlaminar shear stress in sandwich beams. Exp. Mech.

TE D

9 (9), 413–418.

Icardi, U., Sola, F., 2016. Assessment of recent zig-zag theories for laminated and sandwich structures. Compos. Part B-Eng. 97, 26–52.

EP

Jedari Salami, S., Dariushi, S., Sadighi, M., Shakeri, M., 2016. An advanced high-order theory for bending analysis of moderately thick faced sandwich

AC C

beams. Eur. J. Mech. A-Solid. 56, 1–11. Jourawski, D., 1856. Remarques sur la r´esistance d’un corps prismatique et d’une pi`ece compos´ee en bois ou en tˆole de fer `a une force perpendiculaire

a` leur longeur. Annales des Ponts and Chauss´ees 12, 328–351.

Krajcinovic, D., 1972. Sandwich beam analysis. J. Appl. Mech.-T ASME 39 (3), 773–778. 46

ACCEPTED MANUSCRIPT

J. Eng. Ind.-T ASME 97 (3), 873–880.

RI PT

Krajcinovic, D., 1975. Sandwich beams with arbitrary boundary conditions.

Matsunaga, H., 2002. Interlaminar stress analysis of laminated composite beams according to global higher-order deformation theories. Compos.

SC

Struct. 55 (1), 105–114.

Mattei, O., Bardella, L., 2016. A structural model for plane sandwich beams

M AN U

including transverse core deformability and arbitrary boundary conditions. Eur. J. Mech. A-Solid. 58, 172–186.

Phan, C. N., Frostig, Y., Kardomateas, G. A., 2012. Analysis of sandwich beams with a compliant core and with in-plane rigidity — Extended highorder sandwich panel theory versus elasticity. J. Appl. Mech.-T ASME

TE D

79 (4), 041001.

Serpilli, M., Lenci, S., 2008. Limit models in the analysis of three different layered elastic strips. Eur. J. Mech. A-Solid. 27, 247–268.

EP

Sharma, S. R., Rao, D. K., 1982. Static deflections and stresses in sandwich beams under various boundary conditions. J. Mech. Eng. Sci. 24 (1), 11–20.

AC C

Tessler, A., Di Sciuva, M., Gherlone, M., 2009. A refined zig-zag beam theory for composite and sandwich beams. J. Compos. Mater. 43 (9), 1051–1081.

Timoshenko, S. P., 1921. On the correction for shear of the differential equation for transverse vibrations of prismatic bars. Philos. Mag. 38, 744–746.

Tonelli, D., Bardella, L., Minelli, M., 2012. A critical evaluation of mechanical models for sandwich beams. J. Sandw. Struct. Mater. 14 (6), 629–654. 47

ACCEPTED MANUSCRIPT

Tullini, N., Savoia, M., 1995. Logarithmic stress singularities at clamped-free

RI PT

corners of a cantilever orthotropic beam under flexure. Compos. Struct. 32 (1–4), 659–666.

Vidal, P., Polit, O., 2009. Assessment of the refined sinus model for the

SC

non-linear analysis of composite beams. Compos. Struct. 87 (4), 370–381.

Yu, W., Hodges, D. H., 2005. Generalized Timoshenko theory of the varia-

M AN U

tional asymptotic beam sectional analysis. J. Am. Helicopter Soc. 50 (1), 46–55.

Yu, Y. Y., 1959. A new theory of elastic sandwich plates — One dimensional case. J. Appl. Mech.-T ASME 26, 415–421.

Zuo, Q. H., Hjelmstad, K. D., 1998. Piecewise linear warping theory for

AC C

EP

TE D

multilayered elastic beams. J. Eng. Mech. 124 (4), 377–384.

48

ACCEPTED MANUSCRIPT Structural model including transverse core deformability for any plane sandwich beam

New Jourawski-zigzag shear stress estimates

RI PT

Capability of modelling "severe boundary conditions" in any sandwich beam

SC

Computationally efficient structural finite element for linear elastic sandwich beams

AC C

EP

TE

D

M AN U

Critical assessment against continuum plane stress finite element analyses