Generalized viscoelastic model for laminated beams using hierarchical finite elements

Composite Structures 235 (2020) 111794

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Generalized viscoelastic model for laminated beams using hierarchical finite elements

T

⁎

Ezequiel D. Sáncheza, , Liz G. Nallima, Facundo J. Bellomoa, Sergio H. Ollera,b a

ICMASa, INIQUI (CONICET), Engineering Faculty, National University of Salta, Avda. Bolivia 5150, 4400 Salta, Argentina CIMNE International Center for Numerical Method Engineering, Spain, UPC, Technical University of Catalonia (Barcelona Tech), Edif. C1, Campus Nord, Jordi Girona 1–3, Barcelona 08034, Spain

b

A R T I C LE I N FO

A B S T R A C T

Keywords: Viscoelastic model Laminated beams Higher order theory Hierarchical finite element

The formulation of an enriched hierarchical one-dimensional finite element suitable to analyze the rheological behavior of thick arbitrarily laminated beams is presented. The formulation is based on the Equivalent Single Layer (ESL) theory and was developed to allow the use of any high-order beam shear deformation theory (HBST) in a unified approach. A generalized Maxwell model was implemented to analyze the time-dependent behavior of the composite. The finite element employs local Lagrange and Hermitian support functions enriched with orthogonal Gram-Schmidt polynomials and is free of shear locking. The enriched macro elements can be used with very coarse meshes and the precision can be controlled without generating a new mesh. The formulation has been validated with numerical examples of symmetric and non-symmetric laminated beams using the threedimensional finite element program PLCD.

1. Introduction Sandwich structures are widely used as structural supporting elements due to their enhanced mechanical properties. These consist of two thin rigid layers that support much of the flexural load and a lightweight inner core, usually with energy dissipation capabilities. Consequently, these layered composite laminates are well suited for applications where vibration-induced stresses can produce problems such as fracture, fatigue, etc. These problems can be mitigated by the introduction of viscoelastic material core layers of different thickness [1]. The study of the vibration and dynamic response of viscoelastic structures and sandwich panels with viscoelastic cores has been studied over the years. There are several models based on layerwise theories for modelling of laminated structures with viscoelastic damping layers. Moreira et al. [2] developed a layerwise model which can describe the high shear pattern developed inside a thin viscoelastic soft core. Later, they implemented a 4-node facet type quadrangular shell finite element, where the bending stiffness of the facet shell element is based on the Reissner–Mindlin assumptions and the plate theory is enriched with a shear locking protection [3]. Mahmoodi et al. [4] investigated the non-linear free vibrations of viscoelastic beams using the KelvinVoigt model and a multi-scale method was used to analytically raise the non-linear modal shapes and natural frequencies in beams. Araújo et al. [5] presented an optimal design and parameter estimation of frequency ⁎

dependent passive damping of sandwich structures with viscoelastic core and a viscoelastic sandwich finite element model for the analysis of passive, active and hybrid structures [6,7], using a mixed approach, by considering a HSDT to represent the displacement field of the viscoelastic core, and a FSDT for the displacement fields of adjacent laminated face layers. The differential transformation method (DTM) in the frequency domain was employed by Arikoglu and Ozkol [1] to solve the motion equations of sandwich-composite beams including viscoelastic cores. Arvin et al. [8] developed a Finite Element Method (FEM) code to study the vibration frequencies of sandwich beams constituted by external layers of composite material and a viscoelastic core. Moita et al. [9] developed a finite element model for vibration analysis of active–passive damped multilayer sandwich plates, where the elastic layers are modelled using the classic plate theory and the core is modelled using the Reissener–Mindlin theory. They improved this formulation using Reddy’s third-order shear deformation theory [10]. The finite element is obtained by assembly of N ‘‘elements’’ through the thickness, using specific assumptions on the displacement continuity at the interfaces between layers. Galupppi et al. [11] analytically solved the time-dependent problem of a laminated beam with an intermediate viscoelastic layer, whose response is modeled by a Prony-Series. Later, they applied the analysis to a structure composed of elastic glass layers joined by intermediate layers of viscoelastic polymers, under a history of loading and unloading [12]. Optimization of active and passive

Corresponding author. E-mail address: [email protected] (E.D. Sánchez).

https://doi.org/10.1016/j.compstruct.2019.111794 Received 18 September 2019; Received in revised form 4 December 2019; Accepted 9 December 2019 Available online 16 December 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 1. Multilayered composite beam. Geometry and coordinate system.

damping using a new mixed layerwise finite element model was implement by Araújo et al. [13], where results are compared with an alternative optimization model, based on 3D finite elements using ABAQUS commercial package. Ferreira et al. [14] developed a layerwise finite element model using a unified formulation for the analysis of sandwich laminated plates with a viscoelastic core and laminated anisotropic face layers. Liu et al. [15] formulated a layerwise differential quadrature hierarchical finite element (DQHFE) model for the analysis of sandwich laminated plates with a viscoelastic core and laminated anisotropic face layers. The stiffness and mass matrices in these cases were obtained by Carrera’s Unified Formulation (CUF) [16]. Lei et al. [17] established movement equations for Euler-Bernoulli beams using the Kelvin-Voigt model and the standard three-parameter viscoelastic model with velocity-dependent external damping. Li et al. [18] obtained closed form solutions for stresses and strains using the Laplace transform for simply-supported laminated functionally graded (FG) beams, considering viscoelastic interlayers. Dynamic behavior of sandwich composite beams using different shear deformation theories to formulate diverse layer-wise models (LW), implemented through FEM were studied by Loja et al. [19]. Kpeky et al. [20] proposed the modeling of sandwich structures with a soft core using a finite linearhexahedral solid element, combining Automatic Differentiation (AD) with the Asymptotic Numerical Method (ANM). A multiobjective approach for optimization of passive damping for vibration reduction in sandwich structures is presented by Madeira et al. [21,22], for maximization of modal loss factors and minimization of weight of sandwich beams and plates with elastic laminated constraining layers and a viscoelastic core, using the Direct MultiSearch (DMS) solver. Nguyen et al. [23] implemented a triangular finite element that uses the Laplace transform for laminated viscoelastic composite plates based on an efficient higher order zigzag theory (EHOPT). Huang et al. [24] presented two integral finite elements, with compression and shear mechanisms of damping, to model sandwich type structures with soft core. Li et al. [25] developed a semi-analytical method to investigate the natural frequencies and modal shapes of a system of beams interconnected by a viscoelastic layer. The work of Latifi et al. [26] deals with geometrically non-linear transient analysis of sandwich beams with viscoelastic cores and composite laminated external layers. The non-linear dynamic instability of three-layer composite beams with a viscoelastic core subjected to combined lateral and axial loads was also studied by the same authors in Ref. [27]. Guo et al. [28] analyzed a sandwich beam where monolithic viscoelastic core was replaced by two periodically alternating viscoelastic ones to improve the flexural-wave attenuation performance. This work was later extended by Sheng et al. [29] to consider plate sandwich structures. Mustafa [30] studied laminated Timoshenko beams with two identical external layers bounded by a thin adhesive layer. Zhai et al. [31] analyzed the free vibration of five-layer sandwich plates with two viscoelastic cores using the first order shear deformation theory (FSDT); the motion equations were derived using Hamilton's principle and solved by the closed-form Navier method. It should be noted that he inner-core damping layers undergo strong shear strains, due to the relative motion of the layers. For this reason, it is important to use an appropriate kinematics that considers the effect of

the shearing through advanced structural theories [32]. Asik and Tezcan [33], Bennison and Davies [34] and Ivanov [35], among others, showed that an appropriate consideration of the viscoelastic interlayer shear behavior is essential for an accurate modeling allowing an efficient design. The formulation of an enriched hierarchical one-dimensional finite macro-element suitable to analyze the rheological behavior of thick arbitrarily laminated beams, including soft-core sandwich ones, is presented. The formulation is based on the Equivalent Single Layer (ESL) theory which is well suited in design stages or optimization processes where repetitive computations and a good balance between accuracy and resolution speed are required. The formulation is based on the Equivalent Single Layer (ESL) theory and was developed to allow the use of any high-order beam shear deformation theory (HBST) in a unified approach. A generalized Maxwell model was implemented to analyze the time-dependent behavior of the visco-elastic layers. The finite element employs local Lagrange and Hermitian support functions enriched with orthogonal Gram-Schmidt polynomials. The obtained finite element is free of shear locking and thin beams can be studied with the same formulation without resorting to the use of reduced integration [32,36]. The formulation has been validated with numerical examples of symmetric and non-symmetric laminated beams using the three-dimensional finite element program PLCD. 2. Formulation of the mechanical problem 2.1. Kinematics of deformation A laminated beam of length L , width b and total thickness h , as shown in Fig. 1 is considered. An orthogonal Cartesian axis system (x , y, z ) is used, with the axis x oriented along the longitudinal axis of the beam, the plane x − y coincides with the middle plane and the z axis is perpendicular to the middle plane, resulting in a three dimensional domain where 0 ⩽ x ⩽ L,− b 2 ⩽ y ⩽ b 2,− h 2 ⩽ z ⩽ h 2 . The laminated beam is composed of N layers of different elastic and viscoelastic materials. The index k denotes the layer number from the bottom to the top of the beam, therefore the layer k − th lies between z k − 1 ⩽ z ⩽ z k and its thickness is h(k ) . The kinematics of the laminated beam is characterized by the displacements of its midline and occurs in the plane x − z . The components of the displacement field are obtained using different shear deformation theories in the frame of ESL theories, adopting the following general form:

u1 (x , z , t )= u (x , t ) − z u2 (x , z , t )= 0 u3 (x , z , t )= w (x , t )

∂w (x , t ) ∂x

+ fi (z ) ϕ (x , t ) (1)

where u1, u2, u3 are the displacements of any material point in the beam domain along the axes x , y, z , respectively; u, w are the longitudinal ( x axis) and transverse (z axis) displacements of generic points located on the beam longitudinal axis; ϕ is the additional rotation of the normal to the midplane; fi (z ) , i = 1, ...,n represents the shape function that determine the distribution of strains in the beam thickness for different 2

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 4. Response of the constitutive model for a constant imposed deformation.

strains rate of the damper, respectively. The equilibrium condition requires compliance with the following relationship:

Fig. 2. Distribution of transverse shear deformation in beam thickness. Table 1 Shear deformation theories employed. Model

σ (t ) = σ ∞ (t ) + σ i (t )

Shape functions

Reddy [43] Touratier et al. [44] Soldatos et al. [45]

f1 (z ) = z −

4 3

( ) f (z ) = h senh ( ) − z cosh ( ) f2 (z ) = 3

h π

sin

HBT 01

In the particular case of a constant strain imposed over time (ε0 ), the response of the model is shown in Fig. 4 and the stress expression is reduced to the following simple form:

HBT 02

σ (t ) = (C∞ + C1 e−t tr ) ε0

Acronyms

z3 h2 πz h

z h

1 2

(3)

where the relaxation time is tr = η1 C1 and C0 = C∞ + C1 is the elastic constant. The viscoelastic constitutive equation of the generalized multi-axial Maxwell model [37] expressed in global coordinates is given by:

HBT 03

higher order shear deformation theories (HBT) (Fig. 2). If this function is considered to be zero or one, the kinematic expression (Eq. (1)) collapses to the classical Bernoulli theory (TBC) or to the Timoshenko theory (TBT) respectively, as particular cases. The formulation proposed allows the use of any HBT theory but, in this work, three are used which are comprehensive of the different existing types, i.e.: trigonometric, hyperbolic and polynomial ones (Table 1).

[σ ]t + Δt = [σ ]t ·e (−Δt

tr )

− C: [ε ]t ·e (−Δt

tr )

⎛1 + Δt·C1 ⎞ + C: [ε ]t + Δt 2tr·C0 ⎠ ⎝

⎜

⎟

Δt·C1 ⎞ ·⎛1 − 2tr·C0 ⎠ ⎝ ⎜

⎟

(5)

where σ , ε and C are the stress tensor, the strain tensor and the elastic constitutive tensor respectively; while Δt denotes each time interval. For a laminated beam the constitutive equation for each layer can be expressed considering a plane stress state, as follows:

2.2. Constitutive equations The viscoelastic Generalized Maxwell model is represented conceptually in Fig. 3. This model collapse to the simplified Maxwell model when C1 → ∞ and to Kelvin model when C∞ → 0 . The stress state is expressed as:

σ ∞ (t ) = C∞ ε (t ) i σ (t ) = C1 [ε (t ) − ε i (t )] = η1 ε i̇ (t )

(4)

(k ) t + Δt [σxx ] (k ) t (−Δt ] ·e = [σxx

tr (k ) )

(k ) t (−Δt ] ·e − E0(k ) [εxx

tr (k ) )

E1(k ) Δt ⎞ ⎛1 + (k ) ⎟ + E0 (k ) E0 ·tr (k ) 2 ⎠ ⎝

⎜

E (k ) Δt ⎞ (k ) t + Δt ⎛ [εxx ] ·⎜1 − (k ) 1 ⎟ (k ) 2 E 0 · tr ⎝ ⎠

(2)

(6)

(k ) t + Δt [τxz ]

where σ ∞ represents the equilibrium stress with elastic constant C∞, in the spring-damper analog model (Fig. 3), while σ i represents the time dependent viscoelastic stress of the set formed by a spring of elastic constant C1 and a damper of viscous constant η1. In addition, ε represents the total strain of the system, while ε i and ε i̇ are the strains and

(k ) t (−Δt ] ·e = [τxz

tr (k ) )

(k ) t (−Δt ] ·e − G0(k ) [γxz

tr (k ) )

G (k ) Δt ⎞ (k ) t + Δt ⎛ [γxz ] ·⎜1 − (k ) 1 ⎟ G0 ·tr (k ) 2 ⎠ ⎝

(k ) Δt ⎞ ⎛1 + G1 (k ) ⎟ + G0 (k ) (k ) 2 · G tr 0 ⎝ ⎠

⎜

(7)

where (k ) (k ) (k ) 0 ⎤ ⎡ εxx ⎤ ⎡ σxx ⎤ ⎡E σ = ⎢ (k ) ⎥, ε = ⎢ (k ) ⎥, C = ⎢ 0 (k ) ⎥ γ τ G 0 0 ⎦ ⎣ ⎣ xz ⎦ ⎣ xz ⎦

E0(k )

(8)

G0(k )

where and are the elastic constants of each layer, longitudinal and transverse respectively. 2.3. First thermodynamics law The dynamic equilibrium equation based on the first law of thermodynamics for any time t + Δt is given by the following expression:

Pdeformative = Pintroduced − P kinetic

(9)

For a quasi-static analysis, we only take into account the deformative power and the introduced power, with which Eq. (9) remains:

Fig. 3. Generalized uni-axial constitutive model of Maxwell. 3

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al. L

L

∫ ∬ [σxx(k) ·εxx(̇ k) + τxz(k) γxż(k) ] dA dx = ∫ q (x ) u3̇ dx 0

0

A

Table 2 Local support polynomials and first polynomials of Gram-Schmidt for the three nodal variables.

(10)

where q (x ) represents the transverse load distributed along the beam, A (k ) (k ) ̇ γ̇xz are the is the area of the cross section, t denotes the time and εxx strain rates of each layer. The components of the strain tensor and the strain rate tensor are given from the displacements by the following relation:

ε=

∇s u, ε ̇

=

∇s u̇

Lagrange polynomials for u, ϕ

N1u = 0.5(−ξ + 1)

First Gram-Schmidt polynomial for u, ϕ

N3u = −1 + ξ 2

Hermite polynomials for w

N1w = 0.25(2 − 3ξ + ξ 3)

N2u = 0.5(−ξ + 1)

N2w = 0.25(2 + 3ξ − ξ 3) N3w = 0.25(1 − ξ − ξ 2 + ξ 3)

(11)

N4w = 0.25(−1 − ξ + ξ 2 + ξ 3)

From Eq. (1) we have,

u1 (x , z , t ) ⎤ ∂u u=⎡ , u̇ = ⎢ ∂t ⎦ ⎣ u3 (x , t ) ⎥

(12)

and, therefore, do not affect the displacements in the edge nodes [41]. The local support polynomials and the enrichment polynomials are shown in Table 2. Given that in Eq. (15) polynomials Niu (ξ ) for i = 1, 2 and Niw (ξ ) for i = 1, ...,4 are local support functions, and enrichment functions have a null contribution to the displacements in the edge nodes, the following equalities result,

Finally, the components of the strain tensor for an arbitrary point of the k − th beam layer sheet are obtained from the kinematics defined in Section 2.1.

∂ϕ ∂u̇ ∂u ∂u ∂ 2w (k ) (k ) ̇ = ⎛ 1⎞ εxx = ⎛ 1⎞ = − z 2 + fi (z ) ; εxx ∂x ∂x ∂x ⎝ ∂x ⎠ ⎝ ∂x ⎠ ∂ϕ ̇ ∂u̇ ∂2ẇ = − z 2 + fi (z ) ∂x ∂x ∂x

df (z ) df (z ) ∂u ∂u ∂u̇ ∂u̇ (k ) ̇(k ) = ⎛ 1 + 3 ⎞ = i = ⎛ 1 + 3⎞ = i γxz ϕ; γxz ϕ̇ ∂x ⎠ ∂x ⎠ dz dz ⎝ ∂z ⎝ ∂z

N5w = 1 − 2ξ 2 + ξ 4

First Gram-Schmidt polynomial for w

(13)

u1 = u|ξ =−1 = c1u,

u2 = u|ξ =+1 = c2u,

ϕ1 = ϕ|ξ =−1 = c1ϕ,

ϕ2 = ϕ|ξ =+1 = c2ϕ,

w1 = w|ξ =−1 =

(14)

∂w ∂ξ ξ =−1

c1w,

w2 = w|ξ =+1 = c2w, ∂w ∂ξ ξ =+1

= c3w,

= c4w

(17)

3. Hierarchical finite beam element where node 1 correspond to the natural coordinate ξ = −1 and node 2 corresponds to the natural coordinate ξ = +1. The first two generalized displacements for u and ϕ have a specific physical meaning (longitudinal displacements and rotations of the extreme nodes), while for w the first four represent the transverse displacements and their corresponding derivatives at the ends of the beam. The remaining are generalized variables used by the polynomials of GS for the internal approximation of the element. The approximations of the three kinematic nodal variables of the velocity field are given by:

For the approximation of the three kinematic variables of the displacement field u , w and ϕ , finite macro-elements based on the hierarchical version of FEM are proposed, nu

T

u (ξ ) = ∑ Niu (ξ ) ciu = Nu ·cu i=1 nw

T

w (ξ ) = ∑ Niw (ξ ) ciw = N w ·cw i=1 nϕ

T

ϕ (ξ ) = ∑ Niu (ξ ) ciϕ = Nu ·c ϕ i=1

(15)

T

T

(18)

4. Dynamic equilibrium equation Replacing the Generalized Maxwell constitutive model given by Eq. (6) and Eq. (7) in the first thermodynamics law, Eq. (10), we obtain:

ϕ T c = [cu, c w , c ϕ ]T = [c1u, c2u, ...,cnuu , c1w, c2w, ...,cnww , c1ϕ, c2ϕ, ...,cnϕ ]

L

L

∫ ∬ H·(σ t ) e(−Δt tr

N = [Nu , N w , Nu ]T = [N (ξ )1u , ...,N (ξ )unu , N (ξ )1w , ...,N (ξ )nww , N (ξ )1u , ...,N (ξ )unϕ ]T

T

u̇ (ξ ) = Nu ·ċ u; ẇ (ξ ) = N w ·ċ w ; ϕ ̇ (ξ ) = Nu ·ċ ϕ

where − 1 ⩽ ξ ⩽ 1 is the natural coordinate along the beam length L given by ξ = 2x L − 1; ciu, ciw, ciϕ are the generalized displacements; nu , n w , nϕ are the number of terms in the approximation functions for the three nodal variables and Niu, Niw are local support and enrichment polynomials that are described in the following paragraph. In this way, the vector of generalized displacements and the vector of the shape functions used are given by,

0

(16)

(k )

) dAdx +

∫ ∬ H·(C: εt ) ⎛1 + ⎜

0

A

⎝

A

L (k ) e (−Δt tr ) dAdx +

The first form functions are local support ones. For the approximations of u and ϕ identical functions are used, and they are designated by N1u (ξ ) and N2u (ξ ) , which are the classic linear Lagrange polynomials. For the approximation of w , Niw (ξ )(i = 1, ...,4) , Hermite polynomials are used because C1 type functions are needed. Orthogonal Gram-Schmidt (GS) polynomials are added to generate finite enriched macro finite elements, i.e. Niu (ξ ) (i = 3. ..nu , nϕ) and Ni w (ξ ) (i = 5. ..n w ) . The degree of the first polynomial of GS is 2 for the longitudinal displacement (u ) and for the rotation (ϕ ), and 4 for the transverse displacement (w ). These polynomials must have zero contribution to the generalized displacements corresponding to the edge nodes. The other elements of the characteristic orthogonal polynomial sets are generated following the GS procedure [38–40]. In this way, the hierarchical modes contribute only to the internal generalized displacements of the element

⎡ =⎢ ⎣

L

∫ 0

∫ ∬ H·(C:

0 t + Δt

⎤ q (x ) u3̇ dx⎥ ⎦

A

ε t + Δt ) ⎛1 ⎜

⎝

−

E1(k ) (k ) E0 ·tr (k )

E1(k ) (k ) E0 ·tr (k )

·

·

Δt ⎞ ⎟ 2⎠

Δt ⎞ ⎟ dAdx 2⎠

(19)

where

∂u̇ ∂u̇ ∂u̇ H = ⎡ 1, 1 + 3⎤ ∂x ⎦ ⎣ ∂x ∂z

(20)

Replacing the approximation of the kinematic variables of Eq. (15) and Eq. (18) into Eq. (19), simplifying the generalized velocities ċ and working in natural coordinates, the following expression is obtained

4

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al. 1

∫ ⎛⎜∬ Q(1) ·σ te(−Δt tr −1

⎝

(k )

A

1

+

E1(k ) (k ) E0 ·tr (k )

∫ ⎛⎜∬ Q(1) ·(C: εt ) ⎛1 + ⎜

⎝

−1

⎝

A

ct + Δt= Δct + Δt + ct

L ⎟ 2 dξ ⎠

) dA⎞

t int F⏟

·

Δt ⎞ (−Δt ⎟e 2⎠

L ⎟ 2 dξ ⎠

tr (k ) ) dA⎞

1

T T L ⎡ ⎤ + ⎢ (D(1) : A: D(1) + NuB1 Nu ) dξ⎥ct + Δt = 2 ⎣ −1 ⎦ ⏟ K

∫

Δt ext Ft +⏟

2

0 4 Nw L2 , ξξ

0

c¨ t + Δt=

( ) Δc

γ Δtβ

1 Δt 2β

t + Δt

( ) ċ + Δt (1 − ) c¨ − ( ) ċ − ( − 1) c¨

+ 1−

t + Δt

γ β

1 Δtβ

γ 2β

t

1 2β

t

t

t

(28)

The spatial problem is solved by eliminating the residual Δf by successive linearization, where time (t + Δt ) remains constant until the system converges:

(21)

where 2 u u 0 ⎤ ⎡ L N, ξ ⎡ L N, ξ ⎢ ⎥ ⎢ 4 w 0 ⎥; D(1) = ⎢ 0 Q(1) = ⎢− z L2 N, ξξ ⎢ ⎥ ⎢ df z ( ⎢ 2 fi (z ) Nu, ξ i ) Nu ⎥ ⎢ 0 dz ⎣ ⎣L ⎦ − A A A 11 12 13 ⎡ ⎤ A = ⎢− A12 A22 − A23 ⎥ ⎢ A13 − A23 A33 ⎥ ⎣ ⎦

( ) Δc

5.2. Newton-Raphson method [37]

1

∫ N w q L2 dξ −1

ċ t + Δt=

i [F ]t int

0 ⎤ ⎥ 0 ⎥; ⎥ 2 u N ⎥ L ,ξ ⎦

+ K i [c]t + Δt − i [Fext ]t + Δt = i [Δf]t + Δt

(29)

The current equilibrium in the state (i + 1) is forced and expressed through a serial expansion of Taylor truncated in its first variation: i [Δf]t + Δt i (F t int

∂Δf t + Δt

+ i ⎡ ∂c ⎤ ⎣ ⎦

+ Δt )+ + K ct + Δt − Ftext

(22)

· i + 1[Δc]t + Δt ≅

t i ⎛ ∂Fint

⎝

∂c

+ KT −

i + 1 [Δf]t + Δt + Δt ∂Ftext ⎞ i+1

i [J⏟ ]t + Δt

∂c

⎠

·

=0

[Δc]t + Δt = 0

(30)

(A11 , A12 , A22 , A13 , A23 , A33 ) = ∬ E0 (k ) ⎛1 − ⎝ A ⎜

B1 = ∬ G0 (k ) ⎛1 − ⎝ A ⎜

E1(k ) Δt · ⎟⎞ (1, E0(k ) ·tr (k ) 2 ⎠ E1(k ) E0(k ) ·tr (k )

where J is the tangent Jacobian operator. From the previous equation, the displacement increment necessary for the resolution of the system is obtained as follows,

z , z 2 , f (z ), zf (z ), [f (z )]2 ) dA

Δt

· 2 ⎞ [f (z ), z ]2 dA ⎠ ⎟

i + 1 [Δc]t + Δt

(23)

= −( i [J]t + Δt )−1· i [Δf]t + Δt

(31)

The non-linear dynamic equilibrium equation is given by:

(Ftint

5.3. Algorithm

+ Δt + Kct + Δt) − Ftext =0

(24) 1. Initialization of displacements, speeds, accelerations and internal forces.

where K is the stiffness matrix, Fint is the internal force that contain the viscoelastic terms of the non-linearity of the material and Fext is the vector of external forces in the current time step, whose components are detailed below: ϕ T w Fext = [Fuext , Fext , F ext ] L u u ϕ T ] = [01 , 02 , ...,0unu , qN1w , qN2w , ...,qNnww , 01ϕ, 02ϕ, ...,0nϕ 2

c 0 = {0}; ċ 0 = {0}; c¨ 0 = {0}; Fint 0 = {0} 2. Calculation of constant Tangent Jacobian operator J = KT 3. Control of temporary iterations number from n = 1 to tmax established

(25) A. If n Δt > tmax ⇒ END B. If n Δt ⩽ tmax

uu kuw kuϕ ⎤ ⎡k K=⎢ kww kwϕ ⎥ ⎢ sym kϕϕ ⎥ ⎦ ⎣

4 Obtaining the initial residual force and initial correction of the displacements

(26)

where

kijuu=

2 A11 L 2

1

∫ −1 1

kijuϕ= A13 L ∫

−1

kijwϕ= − A23

4 L2

dNiu dN uj dξ ; dξ dξ dNiu dN uj dξ ; dξ dξ 1

∫ −1

d2Niw dN uj dξ ; dξ 2 dξ

kijuw

=

4 −A12 2 L

1

∫

2 dNiu d N jw

Δf 0 = Ftint + Kc 0 − Fext → Δc 0 = −[J]−1Δf 0

dξ ;

5 Determination of displacement, speed and acceleration fields by Newmark Method:

dξ dξ 2 −1 1 2 w d2N w d Ni j 8 kijww = A22 3 dξ ; L dξ 2 dξ 2 −1 1 1 dNiu dN uj 2 L kijϕϕ = A33 L dξ + B1 2 Niu N uj dξ ; dξ dξ −1 −1

∫

∫

γ ⎞ γ ⎞ γ⎞ ⎛ ⎛ ¨ 0; cN = Δc 0 + c 0; cṄ = ⎛⎜ ⎟ Δc 0 + ⎜1 − ⎟ ċ 0 + Δt ⎜1 − ⎟c Δ 2 β⎠ tβ β ⎝ ⎠ ⎝ ⎠ ⎝ 1 1 ⎞ ⎛ 1 − 1⎟⎞ c¨ 0 c¨N = ⎜⎛ 2 ⎟⎞ Δc 0 − ⎜⎛ ⎟ ċ 0 − ⎜ ⎝ Δt β ⎠ ⎝ Δtβ ⎠ ⎝ 2β ⎠

∫

(27)

6. Calculation of residual force and checking by Newton Raphson Method:

5. Non-linear dynamic problem solution

Δf = Ftint + KcN − Fext

The classic solution of the dynamic behavior of structures assumes the “uncoupled temporal-spatial separation of its variables”. In other words, the semi-discrete equation representing the spatial equilibrium at that instant t is solved at each instant of time.

A. If |Δf| ⩽ Tol ⇒ Assignmentc 0 ← cN . Update Fint andn = n + 1. Return to 3. B. If|Δf| ⩽ Tol 7. Calculation of the new displacement correction: Δc = −[J]−1 ·Δf 8. Determination of updated displacements, speeds and accelerations

5.1. Newmark method [37] The temporal integration is approximated with a Newmark scheme were the displacements and velocities in time t + Δt are obtained partially from a system already known in the previous time t .

γ 1 ⎞ ¨N cA = Δc + cN ; cȦ = ⎛β ⎞ Δc + cṄ ; c¨A = ⎛⎜ ⎟ Δc + c 2 β Δ ⎝ Δt ⎠ ⎝ t ⎠ 5

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 5. (a) Simple supported beam (S-S), (b) Clampedclamped beam (C-C) and (c ) Applied load.

Table 3 Properties of materials for three-layer laminated beams. Composite

Layer

h(k ) [mm]

ν (k )

E∞(k ) [MPa]

E1(k ) [MPa]

tr (k ) [seg ]

Symmetric

1 2 3

8.0 4.0 8.0

0.22 0.34 0.22

34.50 14.79 34.50

– 23.45 –

– 1.398 –

Non-symmetric

1 2 3

6.0 4.0 10.0

0.22 0.34 0.22

34.50 14.79 34.50

– 23.45 –

– 1.398 –

Symmetric

1 2 3

6.0 8.0 6.0

0.22 0.28 0.22

34.50 1.479 34.50

– 34.51 –

– 0.852 –

Non-symmetric

1 2 3

6.0 2.0 12.0

0.22 0.28 0.22

34.50 1.479 34.50

– 34.51 –

– 0.852 –

9. Determination of residual strength updated Δf and assignmentcN ← cA , Return to 6. 6. Shear locking study Shear locking is due to the inability of shear deformable elements to accurately model the bending within an element under a state of zero transverse shearing deformation. When thin beams are analyzed, the energy due to tangential stresses must disappear. “p” version of FEM is free of blocking effects if the polynomial degree is chosen to be moderately high. To show this aspect the analysis of a simple supported beam (S-S) of length L = 0.20m , with b = 10mm , with a uniform load q0 = 0.10 KN m (Fig. 5(a) ) is presented. Two laminates are considered, whose material properties are designated as composite A and B in Table 3. Also, different length-to-thickness ratios λ (with λ = L h ) from thick to thin beams are analyzed. The deflection in the center of the beam is computed using the following non-dimensional variable:

w |x = L

2

=

bE (2) w qλ3L

Fig. 7. Distribution of shear stress τxz in x = L 2 , beam S-S and Composite B, for times (a) 1.0s y (b) 5.0s .

(32)

Fig. 6 shows the variation of the deflection as a function of lengthto-thickness ratios, for both laminated composites (i.e. A and B) for t = 1.0s and using the HBT 02 from Table 1. The tangential stresses are also determined employing both

procedures, from the constitutive relations (Eq. (5)) and from the equilibrium equation at the post-processing level, i.e.: (k ) τxz (z ) = −

z

∫−h 2

(k ) ∂σxx dz ∂x

(33)

The determination is made for a beam of composite B material, with λ = 10 . Fig. 7 shows the distribution of shear stress τxz in a cross section located at a distance x = L 20 , for two different time (t = 1.0s and t = 5.0s ). The transverse shear stress computed from the constitutive relationships are labeled “Constitutive”, while the shear stress calculated from the equilibrium equations are labeled “Equilibrium”. For comparison purposes, the results obtained using the three-dimensional finite element program PLCD [42] are also shown, with a mesh of 10,000 hexahedral elements in the discretization and 20 divisions in height, which are labeled as “3D Model”. An examination of the numerical results presented in these graphs shows that the finite hierarchical element does not experience shear blocking, behaves uniformly well for thin and thick beams, and the results of the stress obtained from the equilibrium equation are in excellent agreement with those obtained with the three-dimensional

Fig. 6. Vertical deflection w in x = L 2 , beam S-S as a function length-tothickness ratios. 6

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Table 4 Convergence study in laminated beam (C-C) symmetric (composite A). No. of Gram-Schmidt polynomials 1

2

3

4

5

6

7

8

9

3D Model

HBT 01

w [mm] σxx [MPa] τxz [MPa]

1.815 −0.019 0.058

2.022 −0.019 0.057

2.022 −0.019 0.057

2.021 −0.019 0.057

2.021 −0.019 0.057

2.041 −0.019 0.057

2.041 −0.019 0.057

2.040 −0.019 0.058

2.040 −0.019 0.058

2.027 −0.018 0.058

HBT 02

w [mm] σxx [MPa] τxz [MPa]

1.815 −0.019 0.058

2.023 −0.019 0.057

2.023 −0.019 0.057

2.022 −0.019 0.057

2.022 −0.019 0.057

2.042 −0.019 0.057

2.042 −0.019 0.057

2.041 −0.019 0.058

2.041 −0.019 0.058

2.027 −0.018 0.058

HBT 03

w [mm] σxx [MPa] τxz [MPa]

1.815 −0.019 0.058

2.021 −0.019 0.057

2.021 −0.019 0.057

2.021 −0.019 0.057

2.021 −0.019 0.057

2.041 −0.019 0.057

2.041 −0.019 0.057

2.040 −0.019 0.058

2.040 −0.019 0.058

2.027 −0.018 0.058

8. Validation

model.

In order to show the accuracy and applicability of the laminated beams free shear blocking formulation, different cases are shown in this section. Non-linear time-dependent constitutive model is implemented with 9 GS orthogonal polynomials in all cases. The analysis is carried out considering time intervals of Δt = 0.10s and a distributed load that remains constant in time (Fig. 5(c)). A single finite element with two nodes coinciding with the ends of the analyzed beams is implemented. The internal or hierarchical nodes depend on the number of GS polynomials and are obtained automatically. As previously mentioned, three-dimensional finite element program PLCD is used to validate the results, this code was developed in Fortran language and use the same viscoelastic model implemented in the hierarchical FEM formulation developed. A mesh of 10,000 hexahedral elements is used in the 3D discretization. In the first case, the analysis of a simple supported beam of composite A material, length L = 0.20m and b = 0.01m , under a uniform distributed load q0 = 0.10 KN m is presented. Fig. 9 shows the vertical displacement w in the middle of the beam (x = L 2) , while Fig. 10 shows the normal stresses σxx (x = L 2, z = −h 20) and shear stresses τxz (x = L 10, z = 0) as a function of time, for the three shear deformation theories depicted in Table 1 and the comparison values obtained by the 3D model. Distribution of normal stress σxx (x = L 2) and tangent stress τxz (x = L 10) obtained by equilibrium (Eq. (33)) as function of height are shown in Fig. 11, for a time t = 20.0s using the same comparison cases. In the second case, the same geometry and load setup are used but the boundary conditions is (C–C) and the laminate is asymmetric (Composite B).

7. Convergence analysis With the use of the hierarchical version of the FEM, the number of finite elements is small and is determined by the geometry. In addition, to verify the results, the order of the approach can be selectively increased. This operation is carried out very efficiently because it is not necessary to generate a new mesh, saving computational cost and operator time. For the convergence study, the analysis of a bi-clamped beam (C–C) of length L = 0.20 m with a uniformly distributed load of value q0 = 0.10 KN m (Fig. 5(b) ) and a symmetric laminate (Composite A) is presented. Vertical displacement w in the half span of the beam(x = L 2) , σxx (x = L 2, z = h 20) normal stresses and shear stresses τxz (x = L 10, z = 0) are shown in Table 4, for a number of orthogonal polynomials GS ranging from 1 to 9. Three named variables were determined by employing the three shear deformation theories depicted in Table 1, at timet = 1.0s . In order to compare the results and analyze the convergence, the values obtained using the three-dimensional finite element program PLCD are also shown, with a mesh of 10,000 hexahedral elements in the discretization. The variation of the deflection at the center of the beam as a function of the number of orthogonal polynomials GS is shown in Fig. 8 for composite A at time t = 1.0s , for the HBT 01 theory and the two edge conditions. The edge condition simple supported is labeled as S-S, while the beam clamped at both edges is labeled as C–C. It can be observed that a very good convergence was achieved in both cases, the simple supported condition being more stable and showing faster convergence. This study allows to determinate the number of polynomials necessary to obtain accurate results.

Fig. 8. Vertical deflection w in x = L 2 for beams S-S and C-C, for different numbers of GS polynomials.

Fig. 9. Deflection w in x = L 2 for symmetric S-S beam (Composite A) as function of time. 7

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 10. Stress variation (a) σxx y (b) τxz in S-S symmetric beam (Composite A) as function of time.

Vertical displacement in the center of the beam w (x = L 2) is shown in Fig. 12, while Fig. 13 shows the normal stress σxx (x = L 2, z = −h 20) and the tangent stress τxz (x = L 10, z = 0) as a function of time. Fig. 14 shows the distribution of normal stress σxx (x = L 2) and the tangential stress τxz (x = L 10) obtained also by equilibrium, for a time t = 20.0s using the same comparison cases.

Fig. 11. Stress distribution (a) σxx in x = L 2 and (b) τxz in x = L 10 in S-S symmetric beam (Composite A).

9. Numerical examples The analysis of a asymmetric S-S beam (b = 0.01 m ) with the properties of Composite D (Table 3) with a uniform load q0 = 0.10 KN m , is presented in this first example. Three levels of length-to-thickness ratios are adopted, keeping the total height of the laminate constant and using HBT 01 from Table 1. Normal stresses σxx at point x = L 4, z = −3h 20 , corresponding to the center of the viscoelastic intermediate layer, as a function of time are shown in Fig. 15. Fig. 16 shows the shear stresses τxz (x = L 10, z = −3h 20) obtained by the constitutive equations (Eq. (5)) and by equilibrium equations (Eq. (33)), as a function of time. The second example presents the analysis of a C–C symmetric beam of length L = 0.20m and width b = 0.01m , with the properties of Composite C (Table 3) and summited to a uniform load q0 = 0.10 KN m . Normal stresses σxx (x = L 4) and shear stresses τxz (x = L 10) distribution obtained by the constitutive and equilibrium equations are shown in Fig. 17 and Fig. 18 respectively, for three instants of time and using HBT 02.

Fig. 12. Deflection w in x = L 2 for asymmetric C-C beam (Composite B) as function of time.

theory (HBTs) for an automatic compliance with null shear stress in the free surfaces requirement. The spatial and temporal problems are considered to be independent and solved using a Newton-Raphson and Newmark scheme respectively. The computational implementation of the aforementioned element is also detailed. The implementation of the generalized viscoelastic constitutive model allows a time history analysis of the structure response.

10. Conclusions In this work the formulation of an enriched hierarchical one-dimensional finite macro-element for the analysis of laminated composites with elastic and viscoelastic layers is presented. The formulation proposed allows the use of any high-order beam shear deformation 8

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 13. Stress variation (a) σxx y (b) τxz in beam C-C asymmetric (Composite B) as function of time.

Validation of results for quasi-static analysis of viscoelastic sandwich beams was achieved comparing the results of a finely discretized 3D finite element model with a single hierarchical beam finite macroelement. Enriched macro-elements allow implementing coarse meshes and controlling accuracy without generating a new mesh and the advantages of the p-version are not limited to the greater convergence rate. In fact, with 3D elements, the accuracy of the solution is determined by executing several analyses with different meshes, an expensive and time-consuming process, both because of the computational cost and because of the operator time required to define the different models, including de meshing process. The convergence study has demonstrated that the element proposed requires a moderate number of degrees of freedom for a very good accuracy. Moreover, it can be applied to study thick and thin beams because no shear locking effects were found. Symmetric and asymmetric laminated beams were considered with two edge conditions. A very good agreement was found for both deflection and stresses evolution over time. The stress distribution on the section was very similar for normal stress while shear stress estimation shows a greater divergence (around 3%). The proposed formulation allows a very good estimation of the stress field in the different layers, provided that the values are obtained from the equilibrium condition. Finally, several numerical examples considering quasi-static analysis of viscoelastic sandwich laminated beams were performed to show the temporal evolution of deflections, normal stresses and shear stresses, together with the stresses distribution along the element thickness.

Fig. 14. Stress distribution (a) σxx in x = L 2 and (b) τxz in x = L 10 in beam CC asymmetric (Composite B).

Fig. 15. Normal stresses in x = L 4, z = −3h 20 for S-S beam (Composite D).

CRediT authorship contribution statement

11. Data availability

Ezequiel D. Sánchez: Investigation, Software, Writing - original draft, Writing - review & editing. Liz G. Nallim: Conceptualization, Supervision, Writing - review & editing. Facundo J. Bellomo: Conceptualization, Investigation, Writing - review & editing. Sergio H.

The raw/processed data required to reproduce these findings cannot be shared at this time due to legal or ethical reasons.

9

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Fig. 16. Shear stress (a) equilibrium and x = L 10, z = −3h 20 for S-S beam (Composite D).

(b )

constitutive

in

Fig. 18. Shear stress distribution (a) equilibrium and (b) constitutive in C-C symmetric beam (Composite C).

Acknowledgements Authors acknowledge the financial support of CONICET (PIP Nº 11220150100132-CO), CIUNSa and MINCyT (PICT–Raíces Nº 20152230). References [1] Arikoglu A, Ozkol I. Vibration analysis of composite sandwich beams with viscoelastic core by using differential transform method. Compos Struct 2010;92:3031–9. [2] Moreira RAS, Rodrigues JD. A layerwise model for thin soft core sandwich plates. Comput Struct 2006;84:1256–63. [3] Moreira RAS, Rodrigues JD, Ferreira AJM. A generalized layerwise finite element for multi-layer damping treatments. Comput Mech 2006;37:426–44. [4] Mahmoodi S, Khadem S, Kokabi M. Non -linear free vibrations of Kelvin-Voigt viscoelastic beams. Int J Mech Sci 2007;49:722–32. [5] Araújo AL, Mota Soares CM, Mota Soares CA, Herskovits J. Optimal design and parameter estimation of frequency dependent viscoelastic laminated sandwich composite plates. Compos Struct 2010;92:2321–7. [6] Araujo AL, Mota Soares CM, Mota Soares CA. A viscoelastic sandwich finite element model for the analysis of passive, active and hybrid structures. Appl Compos Mater 2010;17:529–42. [7] Araujo AL, Mota Soares CM, Mota Soares CA. Finite element model for hybrid active-passive damping analysis of anisotropic laminated sandwich structures. J Sandwich Struct Mater 2010;12:397–419. [8] Arvin H, Sadighi M, Ohadi A. A numerical study of free and forced vibration of composite sandwich beam with viscoelastic core. Compos Struct 2010;92:996–1008. [9] Moita JS, Araújo AL, Martins P, Mota Soares CM, Mota Soares CA. A finite element model for the analysis of viscoelastic sandwich structures. Comput Struct

Fig. 17. Normal stress distribution in x = L 4 in beam C-C symmetric (Composite C).

Oller: Conceptualization, Supervision, Writing - review & editing.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 10

Composite Structures 235 (2020) 111794

E.D. Sánchez, et al.

Int J Mech Sci 2017;123:141–50. [27] Latifi M, Kharazi M, Reza Ovesy H. Nonlinear dynamic instability analysis of sandwich beams with integral viscoelastic core using different criteria. Compos Struct 2018;191:89–99. [28] Guo Z, Sheng M, Pan J. Flexural wave attenuation in a sandwich beam with viscoelastic periodic cores. J Sound Vib 2017;400:227–47. [29] Sheng M, Guo Z, Qin Q, He Y. Vibration characteristics of a sandwich plate with viscoelastic periodic cores. Compos Struct 2018;206:54–69. [30] Mustafa M. Laminated Timoshenko beams with viscoelastic damping. J Mathematical Anal Applications 2018;466:619–41. [31] Zhai Y, Li Y, Liang S. Free vibration analysis of five-layered composite sandwich plates with two layered viscoelastic cores. Compos Struct 2018;200:346–57. [32] Nallim L, Oller S, Oñate E, Flores F. A hierarchical finite element for composite laminated beams using a refined zigzag theory. Compos Struct 2017;163:168–84. [33] Asik M, Tezcan S. A mathematical model for the behavior of laminated glass beams. Compos Struct 2005;83:1742–53. [34] Bennison S, Davies P. High-performance laminated glass for structurally efficient glazing, Innovative Light-weight Structures and Sustainable Facades, pp. 1–12, 2008. [35] Ivanov I. Analysis, modelling, and optimization of laminated glasses as plane beam. Int J Solids Struct 2006;43:6887–907. [36] Rango R, Nallim L, Oller S. Formulation of enriched macro elements using trigonometric shear deformation theory for free vibration analysis of symmetric laminated composite plate assemblies. Compos Struct 2015;119:38–49. [37] Oller S. Nonlinear dynamics of structures. CIMNE Springer; 2014. [38] Nallim L, Oller S, Grossi O. Statical and dynamical behaviour of thin fibre reinforced composite laminates with different shapes. Comput Methods Appl Mech Eng 2005;194:1797–822. [39] Nallim L, Oller S. An analytical–numerical approach to simulate the dynamic behaviour of arbitrarily laminated composite plates. Compos Struct 2008;85:311–25. [40] Rango R, Nallim L, Oller S. Static and dynamic analysis of thick laminated plates using enriched macroelements. Compos Struct 2013;101:94–103. [41] Bardell N, Dunsdon J, Langley R. Free vibration analysis of thin coplanar rectangular plate assemblies – Part I: theory, and initial results for especially orthotropic plates. Compos Struct 1996;34:129–43. [42] PLCd research group, PLCd: Non-linear thermo-mechanic finite element code for research-oriented applications, Free access code developed at CIMNE. Available from: http://www.cimne.com/PLCd, 1991-to present. [43] Reddy J. A simple higher-order theory for laminated composite plates. J Appl Mech 1984;51:745–52. [44] Touratier M. A refined theory of laminated shallow shells. Int J Solids Struct 1992;29:1401–15. [45] Soldatos K. A transverse shear deformationtheory for homogeneous monoclinic plates, 1992.

2011;89:1874–81. [10] Moita JS, Araújo AL, Martins PG, Mota Soares CM, Mota Soares CA. Analysis of active-passive plate structures using a simple and efficient finite element model. Mech Adv Mater Struct 2011;18:159–69. [11] Galuppi L, Royer-Carfagni G. Laminated beams with viscoelastic interlayer. Int J Solids Struct 2012;49:2637–45. [12] Galuppi L, Royer-Carfagni G. The design of laminated glass under time-dependent loading. Int J Mech Sci 2013;68:67–75. [13] Araújo AL, Martins P, Mota Soares CM, Mota Soares CA, Herskovits J. Damping optimisation of hybrid active–passive sandwich composite structures. Adv Eng Softw 2012;46:69–74. [14] Ferreira AJM, Araújo AL, Neves AMA, Rodrigues JD, Carrera E, Cinefra M, et al. A finite element model using a unified formulation for the analysis of viscoelastic sandwich laminates. Compos B 2013;45:1258–64. [15] Liu B, Zhao L, Ferreira AJM, Xing YF, Neves AMA, Wang J. Analysis of viscoelastic sandwich laminates using a unified formulation and a differential quadrature hierarchical finite element method. Compos B 2017;110:185–92. [16] Carrera E. Theories and finite elements for multilayered plates and shells: a unified compact formulation with numerical assessment and benchmarking. Arch Comput Methods Eng 2003;10:215–96. [17] Lei Y, Murmu T, Adhikari S, Friswell M. Dynamic characteristics of damped viscoelastic nonlocal Euler-Bernoulli beams. Eur J Mech A/Solids 2013;42:125–36. [18] Li J, Zheng B, Yang Q, Hu X. Analysis on time-dependent behavior of laminated functionally graded beams with viscoelastic interlayer. Compos Struct 2014;107:30–5. [19] Loja M, Barbosa J, Mota Soares C. Dynamic behaviour of soft core sandwich beam structures using kriging-based layerwise models. Compos Struct 2015;134:883–94. [20] Kpeky F, Boudaoud H, Abed-Meraim F, Daya E. Modeling of viscoelastic sandwich beams using solid–shell finite elements. Compos Struct 2015;133:105–16. [21] Madeira JFA, Araujo AL, Mota Soares CM, Mota Soares CA, Ferreira AJM. Multiobjective design of viscoelastic laminated composite sandwich panels. Compos. Part B 2015;77:391–401. [22] Madeira JFA, Araujo AL, Mota Soares CM, Mota Soares CA. Multiobjective optimization of viscoelastic laminated sandwich structures using the Direct MultiSearch method. Comput Struct 2015;147:229–35. [23] Nguyen S, Lee J, Cho M. A triangular finite element using Laplace transform for viscoelastic laminated composite plates based on efficient higher-order zigzag theory. Compos Struct 2016;155:223–44. [24] Huang Z, Qin Z, Chu F. Damping mechanism of elastic–viscoelastic–elastic sandwich structures. Compos Struct 2016;153:96–107. [25] Li Y, Hu Z, Sun L. Dynamical behavior of a double-beam system interconnected by a viscoelastic layer. Int J Mech Sci 2016;105:291–303. [26] Latifi M, Kharazi M, Ovesy H. Effect of integral viscoelastic core on the nonlinear dynamic behaviour of composite sandwich beams with rectangular cross sections.

11