Transverse vibrations of viscoelastic sandwich beams via a Galerkin-based state-space approach

AbstractA new state-space model is formulated for the dynamic analysis of sandwich beams that are made of two thin elastic layers continuously joined by a shear-type viscoelastic (VE) core. The model can accommodate different boundary conditions for each outer layer and accounts for the rate-dependent constitutive law of the core through additional state variables. The mathematical derivation is presented with the Standard Linear Solid (SLS) model (i.e., a primary elastic spring in parallel with a single Maxwell element) and then extended to the generalized Maxwell (GM) model. The kinematics equations are developed by means of Galerkin-type approximations for the fields of axial and transverse displacements in the outer layers, and imposing the pertinent compatibility conditions at the interface with the core. Numerical examples demonstrate the accuracy and versatility of the proposed approach, which endeavors to represent the effects of the VE memory on the vibration of composite beams.


Introduction
Beams and panels with viscoelastic (VE) damping are widely used in state-of-the-art engineering structures, allowing to achieve improved safety, increased durability, and noise and vibration control, as well as reduced costs during the through-life manufacturing, operation, and maintenance of the structural systems (e.g., Lockett 1972;Soong and Dargush 1997;Zhang and Soong 1992;Lee 1997;Rao 2003). This is the case with sandwich beams, where stiff (elastic) external layers are combined with a soft (VE) internal core that has been specifically designed to offer significant damping capabilities in a variety of engineering situations (e.g., Palmeri 2010; Hohe and Librescu 2004;Arikoglu and Ozkol 2010). This in turn has driven the continuous development of new VE materials and laminate configurations and improved and cost-effective manufacturing processes, as well as the formulation of more efficient modeling and design methods. Composite steel-concrete beams (e.g., Adam et al. 1997;Berczyński and Wróblewski 2005;Shen et al. 2011) are another well-known example of multilayered structural members, in which each layer is optimized to achieve an overall better performance and the mechanical properties of the interface (typically represented by shear studs) play a fundamental role in their dynamics (e.g., Kasinos et al. 2015). Different model assumptions for the free undamped vibration of sandwich beams were recently studied by Lenci and Clementi (2012), showing how shear deformations, axial and rotatory inertia, and interface stiffness affect modal frequencies and modal shapes, while Lenci and Warminski (2012) used the multiple-scale method to determine the effects of a nonlinear interface in combination with a pure viscous damping.
The dynamics of three-layered sandwich beams has been addressed in several studies involving the finite element (FE) method (e.g., Soni 1980;Sainsbury and Zhang 1999;Zhang and Sainsbury 2000;Daya and Potier-Ferry 2001;Galucio et al. 2004;Hu et al. 2008;Abdoun et al. 2009;Amichi and Atalla 2009;Won et al. 2013), implementing both linear and nonlinear VE models. For the linear case, the dependance of stiffness and dissipation of the structural elements on the vibration frequency is usually accounted for via complex-valued, frequency-dependent material coefficients. This approach generally involves specialized numerical methods for the solution of nonlinear eigenvalue problems (e.g., Daya and Potier-Ferry 2001), with an increased computational cost for the frequency-domain response functions, as the complex-valued dynamic stiffness matrix has to be formed and decomposed at many different frequencies. A different computational strategy consists of defining the dynamic stiffness matrix for the composite beam (Howson and Zare 2005), but this requires the solution of a transcendental eigenvalue problem with the so-called Wittrick-Williams algorithm (Williams and Wittrick 1983).
In addition to FEs, some researchers (e.g., Vu et al. 2000;Oniszczuk 2000;Oniszczuk 2003; and further works cited in Palmeri and Adhikari 2011) have derived analytical solutions for this class of dynamic problems by modeling the sandwich beam as two parallel Euler-Bernoulli beams continuously connected by a layer of distributed springs. This approach, although very appealing from a theoretical point of view and useful to obtain benchmark solutions, can hardly be applied to cases with different boundary conditions (BCs), functionally graded materials, or both. Such difficulties have been tackled by Palmeri and Adhikari (2011), who proposed a Galerkin-type state-space formulation for the vibration analysis of double-beam systems, in which the computational cost is reduced by (1) adopting simple harmonic functions to discretize the equations of motion and (2) using modal relaxation functions (Palmeri et al. 2004;Muscolino and Palmeri 2007) to represent the VE behavior of the inner layer. In their formulation, however, the inner springs are assumed to provide stiffness in the transverse direction only, so that their contribution is associated with extension or compression of the inner layer, not to any shear strain. On the contrary, such a mode of deformation may play a fundamental role in the dynamics of VE sandwich beams.
Motivated by these considerations, this paper presents a novel Galerkin-type approximation for the fields of both axial and transverse displacements in the outer beams, which allows to represent their relative shear slip. The proposed technique also allows considering (1) any BCs, which in general are different for the two outer layers; and (2) rate-dependent constitutive laws for the inner layer through an enlarged state-space model conveniently built in a reduced modal space. Numerical applications to composite beams with different configurations demonstrate the accuracy and versatility of the proposed approach. Fig. 1(a) shows the composite system under investigation, which consists of two parallel elastic beams of the same length L, continuously connected by an inner layer of shear-type VE springs. In the right-handed Cartesian reference system, z is the abscissa along the beam (which is pointing to the right) and φ is the in-plane rotation (which is positive anticlockwise).

Definition of the Problem
The two elastic beams are prismatic and slender enough for the classical Euler-Bernoulli theory to be valid (i.e., shear deformations and rotatory inertia are negligible). Different material and sectional properties are allowed for the beams, and each of them is fully characterized by modulus of elasticity E r , mass density ρ r , frequency-independent viscous damping ratios ζ r , cross-sectional area A r , and second moment I r , where the subscript r denotes the top (r ¼ 1) and bottom (r ¼ 2) beam, respectively.
The inner layer is homogeneous and, within the limits of the linear viscoelasticity theory, its behavior in shearing is fully characterized by the complex-valued stiffness per unit length, k inn ðωÞ, which depends on the vibration frequency ω; m inn ¼ ρ inn A inn is the mass density per unit length; and d 0 is the distance between their centroids.
Without loss of generality, the Standard Linear Solid (SLS) model is adopted as constitutive law for the inner layer. As illustrated within Fig. 1(b), the model consists of a primary elastic spring (called equilibrium modulus), K 0 ≡ k inn ð0Þ, acting in parallel with a Maxwell element, which in turn is given by a secondary elastic spring, K 1 , in series with a viscous dashpot, C 1 ¼ K 1 τ 1 , where τ 1 is the single relaxation time of the VE material.
For the SLS model, the function k inn ðωÞ can be expressed as where i ¼ ffiffiffiffiffiffi −1 p = the imaginary unit. The shear reaction force, SðtÞ, experienced by the VE inner layer can be related to the associated displacement, sðtÞ (i.e., the internal slippage between top and bottom beam), as where F −1 = the inverse Fourier transform operator; * stands for the convolution operator; and the over-dot means derivative with respect to time t, so thatṡðtÞ is the pertinent velocity (i.e., the rate of variation of the shear strain within the inner layer). The reaction force SðtÞ can be expressed as (Palmeri et al. 2003) where K 0 sðtÞ = the elastic part in the VE constitutive law; K 1 λ 1 ðtÞ = the contribution of the Maxwell element; λ 1 ðtÞ = the additional internal variable, which measures the elongation of the spring, K 1 , and is ruled bẏ The extension to more sophisticated linear VE models such as the generalized Maxwell (GM) model is straightforward (Palmeri et al. 2003) and requires considering more internal variables λ k ðtÞ on the right side of Eq. (3). The associated parameters K i and τ i can be determined from the storage modulus and loss modulus of the VE material, which are measured experimentally for different frequencies of vibration.
It is worth mentioning here that Eqs.
(1)-(4) are formally similar to those used for double-beam systems by Palmeri and Adhikari (2011), with the key difference that shear forces and shear deformations are considered in the current work instead of extensions and compressions of the core.

Kinematics of the System
By resorting to a Galerkin-type approximation, the in-plane vibration of the sandwich beam of Fig. 1(a) can be fully characterized by three independent time-varying fields of displacement: (1) two axial displacements, u 1 ðz; tÞ and u 2 ðz; tÞ, for the top and bottom beams, respectively; and (2) a single transverse displacement, vðz; tÞ, as long as the inner VE layer can be assumed to be incompressible in the transverse direction. It follows that the top and bottom beams will have the same transverse displacements, rotations φðz; tÞ ¼ −v 0 ðz; tÞ, and curvatures χðz; tÞ ¼ −v 0 0 ðz; tÞ, in which the prime means a derivative with respect to the spatial coordinate z.

Axial Assumed Modes
For the axial displacements of the outer beams, the Galerkin shape functions (or assumed modes) can be defined by taking their first n axial modes of vibration, which are the nontrivial solutions of the eigenproblem (Den Hartog 1984;Clough and Penzien 2010): where ϕ r;j ðzÞ and α r;j are the jth pair of eigenfunction and eigenvalue for the rth bar [the eigenpairs satisfying Eq. (5) are offered in Table 1 for different BCs in the axial direction; i.e., Constrained- Once the arrays ϕ r ðzÞ ¼ fϕ r;1 ðzÞ · · · ϕ r;n r ðzÞg ⊤ are defined for the top (r ¼ 1) and bottom (r ¼ 2) beams, with ⊤ denoting the transpose operator, the displacements in the rth outer beam can be expressed as u r ðz; tÞ ¼ ϕ ⊤ r ðzÞ · q r ðtÞ ¼ X n r j¼1 ϕ r;j ðzÞq r;j ðtÞ ð 6Þ in which the n r -dimensional array q r ðtÞ ¼ fq r;1 ðtÞ · · · q r;n r ðtÞg ⊤ collects the Lagrangian coordinates associated with the axial shape functions; and the central dot · = matrix product. Importantly, the BCs in the axial direction can be different for the two outer beams, allowing for maximum versatility.

Transverse Assumed Modes
A convenient array of assumed modes for the transverse deflections can be obtained by taking the first n buckling modes of a homogeneous slender beam, which are the nontrivial solutions of the eigenproblem (e.g., Timoshenko and Gere 2009): where, similar to the previous case, ϕ 0;j ðzÞ and α 0;j are the jth pair of eigenfunction and eigenvalue. Table 2 displays the nontrivial roots of Eq. (7) for different BCs at z ¼ 0 and z ¼ L; i.e., Pinned-Pinned (P-P), Clamped-Free (C-F), Clamped-Pinned (C-P), and Clamped-Clamped (C-C). It is worth noting here that the same assumed modes were exploited by Palmeri and Adhikari (2011) for the dynamic analysis of VE double-beam systems as a way to avoid numerical difficulties due to the presence of hyperbolic functions in the exact modes of vibration of homogeneous Euler-Bernoulli beams. This is advantageous for two reasons: first, it allows avoiding complicated numerical integrations when evaluating the mass and stiffness coefficients (see the appendix); second, it removes any inaccuracy associated with the unbounded nature of the hyperbolic functions, which can affect the higher modes of vibration (Tang 2003).
It also should be noted that, although in the double-beam system studied by Palmeri and Adhikari (2011) the outer beams are connected by transverse springs (and therefore two fields of transverse displacements are required for them), in the present case the inner layer imposes a rigidity constrain in the transverse direction. Thus, a single field of transverse displacements is introduced. Moreover, the associated BCs must be defined in such a way that deflection vðz; tÞ and rotation φðz; tÞ at z ¼ 0 and z ¼ L are allowed if and only if they are compatible with both restraints at the same end of the top and bottom beams. For example, the BCs in the transverse direction for the sandwich beam sketched in Fig. 1(a) must be C-P since the top beam prevents both deflection and rotation at z ¼ 0, while at z ¼ L, the deflection is not permitted by the roller support in the bottom beam and the rotation is allowed in both beams.

Kinetic and Potential Energy
According to Eqs. (6) and (8), total kinetic energy, TðtÞ, and total potential energy, VðtÞ, can be evaluated as the sum of three terms each: Specifically, for the axial vibration of top and bottom layer (r ¼ 1, 2): where μ r ¼ ρ r A r ; and for the transverse vibration (r ¼ 0): where μ 0 ¼ m inn þ ρ 1 A 1 þ ρ 2 A 2 = the mass per unit length of the whole composite beam; EI 0 ¼ E 1 I 1 þ E 2 I 2 = the flexural stiffness of the sandwich beam without shear interaction (i.e., at the limiting condition when K 0 → 0; and sðz; tÞ = the shear slip experienced by the inner layer, given by which couples the axial and transverse movements of the outer beams. In the presence of full interaction between the top and bottom beam (i.e., at the limiting condition when K 0 → ∞), the shear slip goes to zero, sðz; tÞ → 0. Therefore, the axial displacements u 1 ðz; tÞ and u 2 ðz; tÞ become constrained by the transverse displacement vðz; tÞ; i.e., the sandwich beam behaves like a solid beam, whose overall cross section remains plane in bending.
Substituting Eqs. (6) and (8), respectively, into Eqs. (10) and (11) leads to the following expressions for the axial vibration of the outer beams (r ¼ 1, 2): and for the synchronous transverse vibration of the three layers (r ¼ 0): j;kq0;j ðtÞq 0;k ðtÞ ð 14aÞ The sets of mass and stiffness coefficients M ðr;rÞ j;k , K ðr;sÞ j;k , and ΔK ðr;rÞ j;k in Eqs. (13) and (14) couple the jth assumed mode of the rth field of displacements with the kth assumed mode of the rth (M and ΔK) or sth (K) field. The mathematical expressions of these coefficients, which are different with respect to the analogous coefficients presented by Palmeri and Adhikari (2011) for double-beam systems, are provided in the Appendix. In particular, the shear deformations within the core determine the full coupling between the three layers in the sandwich beams considered in this study [Eq. (19b)], while there was no direct coupling between the top and bottom beams in the model investigated by Palmeri and Adhikari (2011).

Generalized Forces
In order to derive the forcing terms in the equations of motion, the external dynamic load fðz; tÞ, transversally applied to the composite beam, is projected onto the assumed modes in the transverse direction. The jth generalized force is then given by Since u 1 ðz; tÞ and u 2 ðz; tÞ are orthogonal to fðz; tÞ, the corresponding generalized forces are identically zero [i.e., Q 1;j ðtÞ ¼ Q 2;j ðtÞ ¼ 0]. Accordingly, the n-dimensional array Q 0 ðtÞ ¼ fQ 0;1 ðtÞ · · · Q 0;n 0 ðtÞg ⊤ fully describes the dynamic loading, while Q 1 ðtÞ ¼ 0 n 1 ×1 and Q 2 ðtÞ ¼ 0 n 2 ×1 , with 0 r×s being a zero matrix with r rows and s columns [Eq. (18b)].

Lagrangian and Modal Equations of Motion
Having defined kinetic energy, potential energy and generalized forces [Eqs. (13)-(15)] in the previous subsections of this paper, the Lagrange equations for the undamped dynamic system can be written as (for r ¼ 0, 1, 2 and j ¼ 1; : : : ; n r ) d dt ∂ ∂q r;j ðtÞ LðtÞ − ∂ ∂q r;j ðtÞ LðtÞ ¼ Q r;j ðtÞ ð 16Þ where LðtÞ ¼ TðtÞ − VðtÞ is the Lagrangian function of the system, and is then cast in the following matrix form: where the arrays uðtÞ and FðtÞ, of size N ¼ n 0 þ n 1 þ n 2 , collect Lagrangian coordinates and generalized forces for the three subsystems, respectively: FðtÞ Importantly, matrix assembly procedures are not required in the proposed formulation, as the mass and stiffness coefficients can be allocated directly to their pertinent positions within the elementary square blocks M ðr;rÞ and K ðr;rÞ (of size n r ) and the (generally) rectangular block ΔK ðr;sÞ , of dimensions n r × n s .
For the purposes of the modal analysis, the matrices M and K can be simultaneously diagonalized through the real-valued eigenproblemω with the orthonormal conditioñ where δ j;k = Kronecker delta symbol, so that δ j;k ¼ 0 when j ≠ k and δ j;k ¼ 1 when j ¼ k;ω j is the approximate jth natural circular frequency of the undamped sandwich beam; and the corresponding modal shape is given by the vector where ΓðzÞ = the 3 × N transformation matrix, defined as follows: Eq. (17) can be further reduced to the modal (decoupled) form where θðtÞ ¼ fθ 1 ðtÞ : : : θ m ðtÞg ⊤ = the array listing the first m modal coordinates of the sandwich beam under investigation, with m ≤ N; Ω ¼ diagfω 1 ; : : : ;ω m g = the m × m spectral matrix; and X ¼ ½x 1 · · ·x 1 = the N × m modal matrix, whose jth column is the jth eigenvectorx j satisfying Eq. (20). After some algebra, one can prove that when the first m modes of the sandwich beam are retained in the analysis, the time-domain dynamic response can be expressed as while the frequency-domain dynamic response at a given abscissa z ¼z can be evaluated as where F = the Fourier transform operator; andZðωÞ ¼ fZ 0 ðωÞ;Z 1 ðωÞ;Z 2 ðωÞg ⊤ = three-dimensional (3D) array collecting the frequency response functions (FRFs) at the selected abscissa z ¼z, given byZ where H 0 ðωÞ is the N × N complex-valued frequency-dependent matrix that collects the FRFs of the Lagrangian coordinates of the system when the first m modes of the system are retained in the analysis. For the undamped case, the matrix H 0 ðωÞ at a given circular frequency ω can be evaluated as where D 0 = 2m × 2m state-space matrix of the dynamic system without damping; and G 0 = the corresponding 2m × N input matrix: and the symbol I s stands for the identity matrix of size s.

Damped Equations of Motion
Eq. (17) in the Lagrangian space and Eq. (23) in the modal subspace do not include energy dissipation. Following a similar approach to Palmeri and Adhikari (2011), both viscous and VE damping can be introduced. Assuming that the former is proportional to the mass density ρ r and modulus of elasticity E r of the rth outer layer (r ¼ 1, 2), and adopting the Rayleigh model (Clough and Penzien 2010; Chopra 2007), the following expression can be devised for the N × N viscous damping matrix in the Lagrangian space: in which the elementary n × n blocks are given by where α K ¼ ðω 1 þω m Þ −1 and α M ¼ω 1ωm α K are the two proportionality coefficients for stiffness and mass, respectively. Since the VE damping is provided by the shear deformations within the core, their effects on the sandwich beam are proportional to the shear stiffness of the inner layer, whose N × N influence matrix in the Lagrangian space is Now, exploiting the same modal transformation of variables as in the previous section [i.e., uðtÞ ¼ X · θðtÞ], the matrices of Eqs. (30) and (31) can be projected onto the reduced modal subspace, which then gives the m × m matrices of modal viscous damping, Ξ ¼ X ⊤ · C · X, and rigidity influence of the inner layer on the modal subspace, B inn ¼ X ⊤ · L inn · X. These new matrices allow including viscous and VE damping directly into the modal equations of motion [Eq. (23) for the undamped case], which now take the form where the m-dimensional array λ 1 ðtÞ ¼ ½λ 1;1 ðtÞ · · · λ 1;m ðtÞ ⊤ has an additional time-varying internal variable for each modal coordinate. Since the spectral matrix Ω is diagonal, it appears from Eq. (32) that the modal coordinates collected by the array θðtÞ are only coupled by the modal matrices Ξ and B inn . If these matrices are diagonal (or if their out-of-diagonal terms are negligible), the dynamic system becomes classically damped in the sense that the modes of vibration are decoupled (Palmeri et al. 2004;Adhikari 2006). Eq. (32) can then be arranged in a more convenient state-space form:ẏ where y 1 ðtÞ ¼ f θðtÞ ⊤θ ðtÞ ⊤ λ 1 ðtÞ ⊤ g ⊤ = enlarged state array (which includes modal displacements, modal velocities, and additional internal variables in the modal subspace), while dynamic matrix D 1 and load influence matrix G 1 are and the subscript 1 means that only one relaxation time (τ 1 ), and therefore only one additional state variable [λ 1 ðtÞ], have been used for the VE behavior of the core. If more than one relaxation time is required, the GM model can be used instead. In this case, the shear force SðtÞ within the inner layer becomes where a viscous dashpot of damping constant C 0 and further (l − 1) Maxwell elements are ideally appended in parallel to the spring-dashpot system of Fig. 1(b). The additional internal variable λ i ðtÞ then describes the elongation of the elastic spring K i within the ith Maxwell element, whose evolution in time is ruled by with τ i as the ith relaxation time of the VE core. The introduction of the GM model has the effect of enlarging the size of the dynamic problem. One can prove that if l Maxwell elements are considered in the analysis, the state-space equations of motion can be posed in the forṁ y l ðtÞ ¼ D l · y l ðtÞ þ G l · FðtÞ ð 37Þ where the expanded array y l ðtÞ ¼ f θðtÞ ⊤θ ðtÞ ⊤ λ 1 ðtÞ ⊤ · · · λ l ðtÞ ⊤ g ⊤ collects the state variables of the system; and dynamic matrix and load influence matrix take the generalized expressions D l ¼ 2 6 6 6 6 6 6 6 6 6 6 6 6 4 The time-domain solution of Eq. (33) or (37) can be sought either with standard techniques [e.g., the classical RK4 Runge-Kutta method (2012)] or specifically tailored numerical schemes (Palmeri et al. 2003;Palmeri and Muscolino 2011).
Alternatively, the dynamic response of the system can be evaluated in the frequency domain. In this case, the solution of Eq. (25) can still be used, provided that the FRF matrix of Eq. (27) is generalized to take into account the effect of viscous and VE damping, as modeled previously. If m modes of vibration are retained in the analysis, and l is the number of relaxation times needed to describe the VE behavior of the core, the N × N complex-valued FRF matrix can be evaluated as which at given frequency ω requires the inversion of a square matrix of size ð2 þ lÞm; i.e., the size of the problem increases linearly with both m and l.

Modal Shapes and Modal Frequencies
In the first stage, the eigenproperties of three undamped (ζ 1 ¼ ζ 2 ¼ 0; l ¼ 0) sandwich beams have been determined. The two outer layers are assumed to be simply supported at their ends; the geometrical and mechanical parameters of the top layer are ρ 1 ¼ 2,000 kg=m 3 , E 1 ¼ 10 GPa, A 1 ¼ 500 cm 2 , and I 1 ¼ 40,000 cm 4 ; and for the bottom layer: The length of the beam is L ¼ 10 m and the distance between the centroids of the outer layers is d 0 ¼ 50 cm; the core is assumed to be massless (i.e., m inn ¼ 0). It is worth noting that, although the outer layers are the same as in some examples presented by Oniszczuk (2003) and by Palmeri and Adhikari (2011), the shear stiffness of the core results in a completely different type of interaction between the top and bottom layers (for instance, all the transverse modes of vibration of the sandwich beams are affected by the stiffness of the inner layer, while this does not happen with the double-beam assemblies studied in those previous papers). Table 3 lists the values of the first 12 (m ¼ 12) modal frequencies predicted by the commercial FE software SAP2000 against those computed with the proposed procedure [Eq. (20)], for three different values of the shear stiffness K 0 of the inner layer; namely, Case A, K 0 ¼ 0 (no shear interaction between the outer layers); Case B, K 0 ¼ 2.0 MN=m 2 (medium core); and Case C, K 0 ¼ 200 MN=m 2 (stiff core). The results of the proposed Galerkin-type approximation have been obtained using eight shape functions for each component (n 0 ¼ n 1 ¼ n 2 ¼ 8), and the bold values in Table 3 denote the predominantly axial modes of the system.
For the first two values of K 0 (Cases A and B), the modal frequencies so computed are in excellent agreement with the values predicted with a relatively coarse FE model (21 spring elements for the inner layer and 40 beam elements for the two outer layers), while the effects of a very stiff core (Case C) are not fully captured Note: Bold values denote the predominantly axial modes of the system. by this FE model. Increasing the fidelity of the FE model (161 spring elements for the core and 320 beam elements for the two outer layers), the values of the modal frequencies so predicted become closer to the values computed with the proposed method (keeping n 0 ¼ n 1 ¼ n 2 ¼ 8), as revealed by the last column of Table 3. It follows that, independent of the relative stiffness of the core, the proposed approach is able to deliver accurate results by using only a few assumed modes for each component (8 × 3 ¼ 24 in total), while a large number of degrees of freedom (2 × 320 ¼ 640 for Case C) may be required by the FE model. Figs. 2 and 3 illustrate the relative contribution of transverse and axial deformations to some selected modal shapes of the undamped sandwich beam [Eq. (21)]. Each of the top subplots [(a), (c), (e), and (g)] shows the vertical deflection v of the sandwich beam (inner dashed line) in a given mode. To make visible the deformations experienced by the core, this has been divided into 30 equally spaced intervals, each one delimited by vertical links, whose inclination in the deformed shape depends on the relative slippage be- Case B with a medium core has been selected for this analysis, and the inspection of Fig. 2 reveals that the first four modes (j ≤ 4) are dominated by transverse deflections, as their amplitude is two orders of magnitude larger that the axial displacements. This is also confirmed by the links, which are all substantially vertical and equally spaced (i.e., they simply appear to move vertically within Fig. 2). The lowest four modes of vibration with dominant axial deformations (j ¼ 6, 8, 11 and 12) are depicted within Fig. 3. Interestingly, the 6th and 11th modes show exactly the same axial displacements in the two outer beams and, being perfectly in phase, there are no shear deformations within the inner layer. Once again, this is confirmed by the links, which are all vertical in Figs. 3(a and e). As a result, these modes are not affected by the stiffness K 0 of the core, as is also evident in Table 3, where the modal frequencies atν ¼ω=ð2πÞ ¼ 111.80 Hz andν ¼ 223.61 Hz are the same for all three cases. On the contrary, the axial displacements experienced by the outer beams in the 8th and 12th mode of vibration are in opposition of phase. For this reason, they also involve shear deformations in the inner layer, which appear in Figs. 3(c and g). Accordingly, the associated modal frequencies increase with the stiffness K 0 of the core.
Figs. 4(a and c) display the absolute value of the FRFs jZ r j of displacements v, u 1 , and u 2 at abscissaz ¼ L=6 [Eq. (26)], all normalized with respect to the pseudostatic value of jZ 0 j when the frequency approaches zero (ν → 0). The dynamic load is uniformly distributed, and the undamped sandwich beams of Cases B and C  Fig. 4 reveals that the axial response of the outer beams (bottom lines) increases significantly due to the shear impedance introduced by the stiffness of the inner layer, which also increases the modal frequencies.

Damped Beams
In the second stage, the effects of VE damping have been quantified. In this case, the top and bottom layers are identical, and they have the same geometrical and mechanical parameters as the top layer (r ¼ 1) in the undamped sandwich beams considered in the previous subsection; the length (L ¼ 10 m) and inner depth (d 0 ¼ 50 cm) are also the same. The viscous damping ratios of the outer layers are taken as ζ 1 ¼ ζ 2 ¼ 0.05. For the inner layer, the mass per unit length is m inn ¼ 12 kg=m, while its VE constitutive law consists of the equilibrium modulus K 0 ¼ 300 kN=m 2 and one or two Maxwell elements (l ¼ 1, 2), with varying stiffness coefficients K 1 and K 2 and relaxation times τ 1 and τ 2 . The BCs are those shown in Fig. 1, given as C-F for the axial deformation of the top layer, C-C for the bottom layer, and C-P for the transverse deflections. Fig. 5 shows four selected modal shapesṽ j ðzÞ, namely: (1) the first two transverse modes (j ¼ 1, 2), in which the amount of axial displacements is negligible in comparison to the transverse  Fig. 5. Selected modal shapes of the damped sandwich beam deflections (and for this reason, the links are all substantially vertical); and (2) the first two axial modes, in which either the bottom layer (j ¼ 5) or the top layer (j ¼ 7) experiences relatively large axial displacements. The number of assumed modes considered for this analysis is n 0 ¼ n 1 ¼ n 2 ¼ 8. As expected, due to the different BCs, the axial profiles of the top and bottom layers can be largely different, which confirms the versatility of the proposed approach. Fig. 6 plots the absolute value jZ 0 ðωÞj and phase angle ∠Z 0 ðωÞ (in rad) for the transverse response of the sandwich beam at abscissaz ¼ L=6, considering five different damping conditions, namely: • Viscous damping only (thin solid lines) • Additional VE damping for the core, modeled with a single Maxwell element (l ¼ 1) of parameters K 1 ¼ 5K 0 and τ 1 ¼ 0.1 s (dashed lines) • Additional VE damping, again with a single Maxwell element (l ¼ 1), with the same stiffness K 1 ¼ 5K 0 and a smaller relaxation time τ 1 ¼ 0.01 s (dotted lines) • Additional VE damping with two Maxwell elements (l ¼ 2), having K 1 ¼ K 2 ¼ 2.5K 0 , τ 1 ¼ 0.1 s and τ 2 ¼ 0.01 s (dotdashed lines) • Additional VE damping with larger stiffness coefficients, K 1 ¼ K 2 ¼ 5K 0 , and the same relaxation times as in the previous case, τ 1 ¼ 0.1 s and τ 2 ¼ 0.01 s (thick dot-dot-dashed lines) A comparison among the different curves demonstrates the beneficial effect of the additional VE damping, which allows reducing the dynamic response at the resonant frequencies, particularly for the lower modes of vibration. Moreover, the values of the resonant frequencies reduce with the relaxation time τ 1 , meaning that the stiffness and energy dissipation are intrinsically interdependent in presence of VE damping. It follows that, for a complex structure made of sandwich elements, selecting the most appropriate mechanical properties of the VE inner layer (in this example, K 0 , K 1 , K 2 , τ 1 , and τ 2 ) is an effective way to optimize and control the global behavior of the structure under dynamic forces, and the proposed state-space formulation can be used to assess the performance of different solutions with less computational effort.
Finally, a convergence study is presented to highlight the influence of the number of assumed modes on the accuracy of the proposed model. For the same composite beam considered previously, the modal frequenciesω fn 0 ;n 1 ;n 2 g j for the first six transverse modes of vibration have been calculated, simultaneously increasing the number of assumed modes n ¼ n 0 ¼ n 1 ¼ n 2 . The convergency measure has been defined as and plotted in Fig. 7(a). Monotonic convergence is observed for all the modal frequencies, and in this example, an accuracy of 0.1% (i.e., e j ¼ 10 −3 ) for the jth transverse mode is achieved with n ¼ j þ 4 assumed modes.
A second analysis has been performed by keeping constant the number of assumed modes for the transverse vibration, n 0 ¼ 12, and simultaneously increasing n 1 and n 2 for the axial vibrations. In order to enhance the interaction between transverse and axial deformations, the elastic modulus of the outer layers has been assumed in this case to be E 1 ¼ E 2 ¼ 0.
and the results are reported in Fig. 7(b) for the first four transverse modes and the first two axial modes. This time, the convergence is not always monotonic, as a result of the complex dynamic interaction between the transverse and axial deformations. For this application, at least n 1 ¼ n 2 ¼ 6 assumed modes are required for the axial vibrations to achieve an accuracy of 0.1%.

Conclusions
A new state-space formulation has been developed and numerically validated for the dynamic analysis of sandwich beams consisting of two parallel Euler-Bernoulli elastic beams continuously connected by a shear-type VE layer. The proposed approach can be used with any boundary conditions of the outer layers and any number l ≥ 0 of Maxwell's elements to represent the constitutive law of the core. The derivation of the state-space model requires the following three steps: 1. The fields of displacements are conveniently described by means of Galerkin-type approximations with simple harmonic functions; namely, the first n 1 and n 2 axial modes of vibration for the two outer layers and the first n 0 bucking modes for the whole composite beam, which are all available in closed analytical form (N ¼ n 0 þ n 1 þ n 2 functions are used in total). In this way, the use of hyperbolic functions is avoided. 2. The undamped equations of motion are derived in the Lagrangian space with a convenient matrix form, in which each element of the N × N matrices of mass and stiffness can be evaluated with a simple numerical integration and then allocated to its position. An assembly procedure is not required. 3. The equations of motion are projected onto a reduced modal subspace (of dimensions m ≤ N) and then cast into an enlarged state-space form, where a set of l × m additional internal variables takes into account the VE damping provided by the inner layer, and the inherent viscous damping of the outer layers can be included.

Appendix. Mass and Stiffness Coefficients
This appendix offers the analytical expressions for the mass and stiffness coefficients introduced in Eqs. (13) and (14) and collected in the n × n block matrices M ðr;rÞ , K ðr;sÞ , and ΔK ðr;rÞ appearing in Eq. (19). The generic mass coefficient M ðr;rÞ j;k for the rth field of displacements in the composite sandwich beam is given by where the mass per unit length μ r takes different expressions for axial displacements (r ¼ 1, 2) and transverse displacements (r ¼ 0). The stiffness coefficients associated with flexural rigidity of the outer beams are given by while those due to their axial rigidity (r ¼ 1, 2) take the expression K ðr;rÞ The remaining sets of stiffness coefficients in Eq. (14b) are those associated with deformations within the inner layer of sheartype springs, which in turn couples the three fields of displacements characterizing the kinematics of the system. These coefficients can be evaluated as follows (with r ¼ 1, 2 and s ¼ 1, 2): Importantly, all the integrands appearing in the right side of these expressions are the product of simple trigonometric functions of the abscissa z, and therefore the evaluation of these coefficients can be done in a closed form. m = number of modal coordinates; m inn = mass of the unit length of the core; N = total number of assumed modes (N ¼ n 0 þ n 1 þ n 2 ); n r = number of assumed modes in the rth component; Q, Q = generalized force and array of generalized forces; q, q = Lagrange coordinate and array of Lagrange coordinates; S = shear force within the core; s = internal slippage between top and bottom beam; T = kinetic energy; t = time instant; u = axial displacement; u = superarray of Lagrangian coordinates; V = potential energy; v = transverse displacement; v = array of the modal shape components; X = modal matrix; x = eigenvector; y = array of the state variables; Z = array of FRFs at a given abscissa; z = abscissa along the beam; 0 = zero matrix; α K , α M = proportionality coefficients for the Rayleigh damping; Γ = transformation matrix; δ = Kronecker delta; ζ = viscous damping ratio; θ, θ = modal coordinate and array of the modal coordinates; λ, λ = additional internal variable and array of additional internal variables; μ = mass per unit length; ν = modal frequency (Hz); Ξ = viscous damping matrix in the modal subspace; ρ = mass density; τ = relaxation time; ϕ, ϕ = assumed mode and array of assumed modes; φ = rotation; χ = curvature; Ω = spectral matrix; and ω,ω = frequency of vibration and modal frequency (rad=s).