# Model reduction methods for viscoelastic sandwich structures in frequency and time domains

Souhir Zghal, Mohamed Lamjed Bouazizi, Noureddine Bouhaddi, Rachid Nasri

► To cite this version:

Souhir Zghal, Mohamed Lamjed Bouazizi, Noureddine Bouhaddi, Rachid Nasri. Model reduction methods for viscoelastic sandwich structures in frequency and time domains. Finite Elements in Analysis and Design, 2015, 93, pp.12 - 29. 10.1016/j.finel.2014.08.003. hal-01511270

HAL Id: hal-01511270

https://hal.science/hal-01511270

Submitted on 20 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Model reduction methods for viscoelastic sandwich structures
in frequency and time domains

Souhir Zghal a,b,c, Mohamed Lamjed Bouazizi b, Noureddine Bouhaddi c, Rachid Nasri a

<sup>a</sup> National School of Engineers of Tunis (ENIT), Applied Mechanics and Engineering Laboratory, University of Tunis El Manar, BP 37,1002 Belvédère, Tunis, Tunisia

<sup>b</sup> Preparatory Engineering Institute of Nabeul (IPEIN), Research Unit of Structural Dynamics, Modeling and Engineering of Multi-Physics Systems, University of Carthage, 8000 M'rezgua, Nabeul, Tunisia

<sup>c</sup> FEMTO-ST Institute UMR 6174, Department of Applied Mechanics, University of Franche-Comté, 24 Chemin de L'Epitaphe, 25000 Besançon, France

This paper deals with modeling and model reduction methods intended to sandwich structures with viscoelastic materials. The modeling step is carried out by combining the First order shear deformation theory (FSDT) with the Golla-Hughes-McTavish (GHM) model. The GHM model introduces auxiliary coordinates to take into account the frequency dependence of viscoelastic materials which, combined with the finite element method (FEM), leads to large order models. This paper focuses on the use of model reduction methods. The reduced models compared to the full model are illustrated by three numerical examples in order to outline the performance, the practical interest of these methods and their validity domains.

Keywords:

Sandwich

GHM viscoelastic model

FEM

Model reduction methods

Frequency analysis

Time analysis
# 1. Introduction

The use of viscoelastic [1,2] sandwich structures [3] has been regarded as a convenient strategy for many industries such as aeronautics, marines and automotives. In fact, these structures present a high way of vibration control in terms of lightweight and high specific stiffness, especially when they incorporated viscoelastic materials.

Several theories [4–7] were developed in order to approximate the displacement and the mechanical deformation of such structures. One of the well-known and useful theories is the classical theory of plates (CPT), which assume that a plane section initially normal to the midsurface before deformation remains plane and normal to that surface after deformation. Hence, this theory neglects the effect of shear deformations and leads to inaccurate results for laminated plates. So, it is obvious that transverse shear deformations have to be taken into account in the analysis. Thus, the first order shear deformation theory (FSDT) introduced by Reissner and Mindlin [4,7] takes into account this effect and assumes a linear variation of the midplane displacements through the thickness of the structure. This method has a significant advantage due to its simple implementation and low computational cost. Another laminated theory based on Reddy's refined [8] high order shear deformations theory (HSDT) which includes both bending and shear effects was carried out by Ferreira et al. [9], and by Chugal and Shimpi [10]. Unfortunately, this method requires a prohibitive computational time which is undesirable for such applications. Some other researchers [11,12] have used the layerwise theory for modeling the sandwich structures. Indeed, this theory assumes a displacement field in the form of zigzag along the thickness of the structure, allowing a kinematic description of each layer as a piecewise linear function. In addition, this theory is applicable to both thin and thick structures. Nevertheless, when the study is intended for thin structures, the first order shear deformation theory (FSDT) presents a suitable choice for the modeling of sandwich structures favored by its simple implementation in most finite element codes.

However, these structures exhibit viscoelastic damping, which combines viscous and elastic character. Hence, this dual character leads to a complicated behavior which requires a correct modeling approach. More recently, Golla, Hughes and Mc Tavish [13,14] have proposed the so-called GHM model. This model provides an effective method which includes viscoelastic damping through the addition of auxiliary coordinates called dissipation coordinates as a sum of elementary mini-oscillators.

Furthermore, the GHM model, combined with the finite element method (FEM) [15], allows the introduction of viscoelastic material properties through element mass, stiffness, and damping matrices. The addition of internal mini-oscillators for each viscoelastic finite element allows a general description of frequency-dependent viscoelastic materials properties behavior. The main advantage of this method consists in its efficient modeling of viscoelastic material behavior; but its major drawback is the largely finite element dimension system which requires a prohibitive computing time. Consequently, a model reduction should be applied to the augmented GHM model.

The present paper proposes an alternative of model reduction such as dynamic [16,17], Guyan [18,19], modal and modal in physical space (SEREP) [20–23] reduction methods for this problem. The first one is based on the elimination of unwanted

variables; it partitions the full degree-of-freedom (dofs) into master and slave dofs and uses the modal properties of the slave part of the structure when the master dofs are grounded. Hence, the derived slave modes are operated to enrich the dynamic basis leading to a drastic reduction method. The simplest yet very useful model reduction method is the well-known Guyan reduction method. It is a particular case of dynamic reduction method according to which the inertia associated with the slave coordinates is neglected; only master dofs are retained. Thereby, the unwanted variables are removed, leading to a reduced model which is a subset of the original system in a restricted range of frequency. However this method is limited by its validity domain [24,25]. Another reduction method is the frequently used modal reduction method according to which the derived modes associated with the undamped structure are incorporated in the GHM damped model, yielding an exact transformation basis. This basis restitutes correctly the undamped modes of the original system leading to a drastic reduction. The modal reduction method can expand the projection from generalized coordinates system to the physical coordinates system, leading to another strategy of reduction called modal reduction in physical space method. This method restitutes also the first modes of the undamped structure and partitions the modal basis into master and slave dofs. This leads to several cases which will be tested examining both the number of retained modes and the number of master dofs.

On the other hand, the modeling of viscoelastic sandwich structures has attracted many researchers, but only a few papers have dealt with the GHM model [26,27]. However, these papers remain limited mostly to frequency domain analysis with major uses of the space state modal reduction method for model reduction. In fact, Trindade et al. [28] and De Lima and Rade [29] used the modal reduction in their studies frequently. It consists in transforming the second-order equation of motion into an equivalent first-order form (space-state model). Unfortunately, this method leads generally to a space state model of dimension at least the double of the total dimension of the GHM model (2 N) and the quadruple dimension of the structural dofs which requires a prohibitive time of calculations.

Therefore, the application of the proposed reduction methods, which are often used with the undamped structures, combined with the GHM model allows one to add the effects of viscoelastic components to the sandwich structures without increasing the order of the finite element models. Furthermore, these reduction methods can be applied to sandwich structures described kinematically by the other mentioned theories.

In this paper, both the theory related to the implementation of the FSDT theory combined with the GHM method and the theory related to its reduction methods are presented. Numerical simulations applied to beam, plate and non-linear assembled beams in both frequency and time domains are also illustrated. These examples will highlight the performance of reduction methods and its practical interest in the dynamic analysis of viscoelastically damped sandwich structures.

# 2. Three-layer viscoelastic finite element model

Multilayer structures are typically used for its light-weight, high specific stiffness and strength values in many engineering
fields. In fact, there are attempts to replace components with classical materials (steel, concrete) by laminated materials, notably sandwich structures. Hence, the modeling of such structures has been of particular interest of many studies [7,9–11]. In this paper, the considered sandwich structure is constituted by three laminated materials: a core generally formed by viscoelastic material of thickness $h_c$, incorporated between two elastic faces of thickness $h_{f1}$ and $h_{f2}$, respectively. The studied sandwich panel is assumed to have length L, width b and total thickness h as shown in Fig. 1.

The modeling process of the three-layered sandwich structure is based on the following assumptions:

* The sandwich is a laminate made of a stack of permanently combined layers. There is no slip or delamination between the layers, they are perfectly bonded. Consequently, the continuity of displacements along the interfaces between the layers is considered.

* The sandwich is a homogeneous material on a macro-scale level, but its properties depend in turn on the properties of each layer.

The lamina are macroscopically homogeneous, isotropic and linearly elastic.

* Both extensional and bending deformations are considered.

It should be noted that when the lamina core is made of viscoelastic material, an appropriate model will be used to model such behavior.

The kinematic model of the sandwich structure is based on the first shear deformation theory (FSDT) of Reissner–Mindlin [4,7], which assumes that a plane section and perpendicular to the midplane of the structure before deformation remains plane but not necessary perpendicular to the midplane after deformation.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Sandwich structure geometry.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Serendip finite element in (a) elementary coordinates and (b) local coordinates.</figcaption></figure>

This theory takes into account the effect of transverse shear deformations in both faces and core. Hence, the displacement field of a sandwich laminate structure can be expressed as

$$ \left\{ \begin{array}{l} u(x, y, z, t) = u_0(x, y, t) + z\psi_x(x, y, t) \\ v(x, y, z, t) = v_0(x, y, t) + z\psi_y(x, y, t) \\ w(x, y, z, t) = w_0(x, y, t) \end{array} \right. \qquad (1) $$

where u, v and w represent the displacements along the axes x, y and z, respectively; $u_0(x, y, t)$, $v_0(x, y, t)$ and $w_0(x, y, t)$ denote the corresponding midplane displacements in the (x,y,z) directions. z is generally the thickness of the structure along the axis (z). $\psi_x(x, y, t)$ and $\psi_y(x, y, t)$ are the rotations of normals to midplane about the y and x axes.

This theory is well applicable for thin and moderately thick plates and allows the compromise: good capacity of prediction/ moderate computational time for large manufactories investigations. Besides, it offers the feasibility of easy implementation in many finite elements codes.

Thus, the finite element formulation uses an eight-node shell finite element with five dofs per node called the Serendip element [15]. The choice of this element is based on the investigations realized by Chee [30], which proved that this element provides an excellent performance for the modeling of composites structures, notably sandwich structures. Furthermore, this element is adapted for the majority of laminated theories, especially first order shear deformation theory (FSDT).

It is a quadratic element belonging to the isoparametric elements family and it uses a bilinear shape function whose coordinates in elementary and local system are presented in Fig. 2.

Therefore, the displacement field can be discretized in local coordinates as follows:

$$ \left\{ \begin{array}{l} u(\xi, \eta, z, t) \\ v(\xi, \eta, z, t) \\ w(\xi, \eta, z, t) \end{array} \right\} = [A(z)]_{(3 \times 5)}[N(\xi, \eta)]_{(5 \times 40)}\{u_e(t)\}_{(40 \times 1)} \qquad (2) $$

where $[A(z)]_{(3 \times 5)}$ is the matrix of z coordinates along the thickness axis; $[N(\xi, \eta)]_{(5 \times 40)}$ is the shape functions matrix; $\{u_e(t)\}_{(40 \times 1)}$ is the elementary nodal displacement vector.

Based on the hypotheses of the stress-states assumed for each layer, the stress–strain relations can be obtained and kinetic energies of the sandwich plate finite element are formulated [31]. Then, the variational Hamilton principle is used considering the nodal displacements and rotations as generalized coordinates, leading one to derive the element's stiffness and mass matrices as
follows:

$$[M_e] = \int_{V_e} \rho^{(k)}[N]^T [A]^T [A][N] dV_e \quad (3a)$$

$$[K_e] = \sum_{k=1}^{N_c} \int_{\xi=-1}^{1} \int_{\eta=-1}^{1} \int_{z=z_k}^{z_{k+1}} ([B_b]^T [C_b^e]^{(k)} [B_b] + [B_s]^T [C_s^e]^{(k)} [B_s]) dz d\eta d\xi \quad (3b)$$

where the subscript k denotes the kth layer of the laminated structure; $\rho^{(k)}$ is the corresponding lamina density; and $N_c$ denotes the number of laminated layers. Herein, $N_c$ is considered equal to 3. $z_k$ denotes the thickness of the kth sandwich layer; $[B_b]$ and $[B_s]$ refers to the strain-displacement matrices where the extensional, bending, plane shear and transverse shear effects are uncoupled separately.

$[C_b^e]^{(k)}$ and $[C_s^e]^{(k)}$ refers to the strain-stress matrices associated to each layer k where the extensional, bending, plane shear effect (subscript b) and transverse shear effect (subscript s). $dV_e$ indicates the elementary variation volume and J is the Jacobian matrix.

After derivation of the element mass and stiffness matrices using the Gaussian quadrature integration, the corresponding global matrices are assembled accounting for the connectivity using the standard assembling procedure and the equation of motion of undamped structure is established as follows:

$$[M][\dot{q}] + [K][q] = \{F\} \quad (4)$$

where $[M] \in R^{N \times N}$ is the mass matrix (symmetric and positive definite), $[K] \in R^{N \times N}$ is the stiffness matrix (symmetric and non-negative definite), $[q]$ is the displacements vector and $\{F\} \in R^{N \times 1}$ is the external load vector.

Nevertheless, when the sandwich structure is made of viscoelastic material, this equation of motion is unable to describe the frequency-dependence of such materials. Indeed, it omitted the damping effect. Hence, the use of a consistent model across a broad range of frequencies should be considered.

Several approaches are presented in the literature to describe this behavior such as the anelastic displacement fields model proposed by Leisuture [32], Fractional derivatives models proposed by Bagley and Torvik [33] and especially the Golla-Hughes-Mc Tavish (GHM) model [13,14]. Hence, the GHM model can be developed for direct incorporation into the finite element method.

### 3. GHM viscoelastic approach

For a sandwich structure incorporating viscoelastic materials, the stiffness matrix can be decomposed as follows:

$$[K] = [K_e] + [K_v(s)] \quad (5)$$

where $[K_e]$ is the stiffness matrix corresponding to the purely elastic layers and $[K_v(s)]$ is the stiffness matrix associated with the viscoelastic layer. The inclusion of the frequency-dependent behavior of the viscoelastic material can be made by generating $[K_v(s)]$ for specific types of elements (beams, plates, etc.) considering initially constant moduli ($E(s)$ or $G(s)$). Then, using the elastic-viscoelastic correspondence principle [34,35], these moduli are factored out of the stiffness matrix reflecting the frequency dependence of viscoelastic materials.

Hence, the viscoelastic stiffness can be written as

$$[K_v(s)] = G(s)[\bar{K}_v] \quad (6)$$

Golla-Hughes and Mc Tavish [13,14] introduced the so-called GHM model to describe the shear modulus of viscoelastic structure as a series of mini-oscillator terms:

$$G(s) = G_0 + \sum_{i=1}^{N_G} \alpha_i \frac{s^2 + 2\zeta_i \omega_i s}{s^2 + 2\zeta_i \omega_i s + \omega_i^2} \quad (7)$$

where $G_0$ is the static modulus; $s$ is the Laplace complex variable; $(\alpha_i, \zeta_i, \omega_i)$ are the parameters of the ith mini-oscillator, and $N_G$ is the number of mini-oscillators. The parameters $(\alpha_i, \zeta_i, \omega_i)$ are identified from the experimental fit curves of the corresponding viscoelastic material [26,34]. In fact, different viscoelastic materials have different frequency dependence and so have a different number of terms $N_G$ of the GHM fit.

Substituting Eq. (6) into Eq. (5) and replacing $G(s)$ by its expression, the equation of motion can be written as follows:

$$\begin{gathered} (s^2[M] + s[D] + [K_e] + [K_V^0])\{q(s)\} + [K_V^0] \sum_{i=1}^{N_G} \frac{s^2 + 2\zeta_i \omega_i s}{s^2 + 2\zeta_i \omega_i s + \omega_i^2} \\ \{q(s)\} = \{F(s)\} \end{gathered} \quad (8)$$

Now, by adding extra-coordinates $\{z_i\}(1, ..., N_G)$ called dissipation coordinates as

$$\{z_i(s)\} = \left\{ \frac{\omega_i^2}{s^2 + 2\zeta_i \omega_i s + \omega_i^2} \right\} \{q(s)\} \quad (9)$$

The equation of motion may be rewritten, in the Laplace domain, as two coupled second order equations:

$$\{s^2[M] + s[D] + [K_e] + [K_V^\infty]\}\{q\} - \alpha[K_V^0]\{z\} = \{F\} \quad (10a)$$

$$s^2\{z\} + 2s\zeta\omega(z) - \omega^2\{q\} + \omega^2\{z\} = \{0\} \quad (10b)$$

where $[K_V^0] = G_0[\bar{K}_V]$ and $[K_V^\infty] = [K_V^0](1 + \sum_{i=1}^{N_G} \alpha_i)$ are, respectively, the static or low-frequency stiffness matrix and the dynamic or high-frequency stiffness matrix corresponding to the viscoelastic layer.

The matrix $[D]$ represents the structural damping of the structure without the viscoelastic effect.

After some manipulations and back to time domain, the following equation of motion in the Laplace domain is obtained:

$$\left\{ s^2 \begin{bmatrix} [M] & 0 \\ 0 & [M_z] \end{bmatrix} + s \begin{bmatrix} [D] & 0 \\ 0 & [D_z] \end{bmatrix} + \begin{bmatrix} [K_q] & [K_{qz}] \\ [K_{qz}^T] & [K_z] \end{bmatrix} \right\} \begin{Bmatrix} q(s) \\ z(s) \end{Bmatrix} = \begin{Bmatrix} F(s) \\ 0 \end{Bmatrix} \quad (11)$$

Or in compact form:

$$\{s^2[M_G] + s[D_G] + [K_G]\}\{q_G\} = \{F_G(s)\} \quad (12)$$

The derived second-order time domain equation of motion is expressed as

$$[M_G]\{\ddot{q}_G\} + [D_G]\{\dot{q}_G\} + [K_G]\{q_G\} = \{F_G\} \quad (13)$$

where $[M_G]$; $[D_G]$ and $[K_G] \in R^{n_G \times n_G}$, with $n_G = N(1 + N_G)$, are, respectively, the mass, damping and stiffness matrices of the global viscoelastic GHM model expressed as follows:

$$\begin{align} [M_G] &= \begin{bmatrix} [M] & 0 & \cdots & 0 \\ 0 & \frac{\alpha_i}{\omega_i}[K_V^0] & 0 & \cdots \\ \vdots & 0 & \ddots & 0 \\ 0 & \cdots & 0 & \frac{\alpha_{N_G}K_V^0}{\omega_i} \end{bmatrix}, & [D_G] &= \begin{bmatrix} [D] & 0 & \cdots & 0 \\ 0 & \frac{2\alpha_i}{\omega_i}[K_V^0] & 0 & \cdots \\ \vdots & 0 & \ddots & 0 \\ 0 & \cdots & 0 & \frac{2\alpha_{N_G}K_V^0}{\omega_i} \end{bmatrix} \nonumber \\ [K_G] &= \begin{bmatrix} [K_e] + [K_V^\infty] & -\alpha_1[K_V^0] & \cdots & -\alpha_{N_G}[K_V^0] \\ -\alpha_1[K_V^0]^T & \alpha_1[K_V^0] & 0 & \cdots \\ \vdots & 0 & \ddots & 0 \\ -\alpha_{N_G}[K_V^0]^T & \cdots & 0 & \alpha_{N_G}[K_V^0] \end{bmatrix}, & \{F_G\} &= \begin{Bmatrix} F \\ 0 \\ \vdots \\ 0 \end{Bmatrix}, & \{q_G\} &= \begin{Bmatrix} q \\ z_1 \\ \vdots \\ z_{N_G} \end{Bmatrix} = \begin{Bmatrix} q \\ z \end{Bmatrix} \quad (14) \end{align}$$

The dissipative coordinate $\{z\}$ appears as an augmenting state variable which increases considerably the order of the differential equation of motion. In fact, the dimension of the system is at least doubled and the computational time is notably increased. This motivates the use of model reduction methods, as an alternative solution for this problem, which will be presented in the following section.
## 4. Model reduction methods

Model reduction is necessary to reduce the high order finite element models to a smaller size so that direct dynamic analysis can be performed. Several model reduction methods have commonly been used including dynamic [16,17], Guyan [18,19], modal and modal in physical space [20–23] reduction methods. They can be classified either by type of reduction approach or by type of reduced space coordinates.

The former can be further divided into two categories as the elimination dofs approach and the modal projection approach. Indeed, the elimination dofs approach is based on the partition of the full dofs of the structure into ‘m’ master and ‘s’ slave dofs. In the reduction process, master dofs are retained and slave dofs are removed. This category includes notably dynamic and Guyan reduction methods. For the modal reduction approach, a partition of the modal matrix for the associated undamped model is established. The lowest modes are retained and all else are removed. Both modal reduction and modal reduction in physical space belong to this category. The concept of physical, generalized and hybrid coordinates will be clarified in the latter type of reduction.

In fact, based on the type of coordinates retained as the reduced order coordinates, the existing model order reduction methods fall into three basic categories:

* Physical coordinates reduction.

* Generalized coordinates reduction.

* Hybrid coordinates reduction.

In the physical coordinates model method, the reduced model is obtained by removing part of the physical coordinates of the full model. Thus, the coordinates of the reduced model are a subset of the full model. This is the most straightforward model reduction among the three categories. Guyan and modal in physical space reduction methods belong to this type of coordinates.

In the generalized coordinates model reduction, all the coordinates that are not physical coordinates are generally referred to as generalized coordinates. The modal reduction method is one of the frequently used generalized coordinates.

The hybrid coordinates model reduction uses a combination of physical and generalized coordinates. Thus, this technique provides a good representation of the dynamic behavior of the sandwich structures. One very useful method belonging to this type is the dynamic reduction method.

Each method is attached to the definition of a transformation matrix $[T] \in R^{n_G \times n_c}$, relating the $n_G$ full dofs of the viscoelastic sandwich structure to the $n_c$ reduced dofs where $n_c \ll n_G$. Thereby, the displacement vector $\{q_G\}$ can be written as a linear combination of the subspace elements presented by the columns of $[T]$ as

$$ \{q_G\} = [T]\{q_c\} \qquad (15) $$

where $\{q_G\} \in R^{n_G \times 1}$ is the displacement vector of the full GHM model. $\{q_c\} \in R^{n_c \times 1}$ is the vector of reduced coordinates through the projection on $[T]$. $[T] \in R^{n_c \times n_c}$ is the transformation matrix.

This transformation takes various forms depending on the used reduction technique.

The equation of motion in full space Eq. (12) is then written in reduced space as follows:

$$ ([T]^T[M_G][T])\{\ddot{q}_c\} + ([T]^T[D_G][T])\{\dot{q}_c\} + ([T]^T[K_G][T])\{q_c\} = [T]^T\{F_c\} \qquad (16) $$

Hence, the reduced mass, stiffness and damping matrices, as well as the reduced load vector can be written as

$$ [K_c] = [T]^T[K_G][T] $$

$$ \begin{aligned} [M_c] &= [T]^T[M_G][T] \\ [D_c] &= [T]^T[D_G][T] \\ \{F_c\} &= [T]^T\{F_G\} \end{aligned} \qquad (17) $$

For each type of reduction process, the transformation from the full space to the reduced space is established through the partition of structural displacement vector $\{q\}$ into two subvectors as follows:

$$ \{q\} = \begin{cases} q^m \\ q^s \end{cases} \qquad (18a) $$

where the subscript $m$ is related to the master dofs and the subscript $s$ is related to the slave dofs.

Following the master and slave dofs partitioning and assuming that only the master dofs are loaded, the external load vector $\{F\}$ can be written as

$$ \{F\} = \begin{cases} F^m \\ 0 \end{cases} \qquad (18b) $$

Consequently, the dynamic equilibrium of the associated undamped system can be expressed as

$$ \begin{bmatrix} K^{mm} & K^{ms} \\ K^{sm} & K^{ss} \end{bmatrix} - \lambda \begin{bmatrix} M^{mm} & M^{ms} \\ M^{sm} & M^{ss} \end{bmatrix}) \begin{Bmatrix} q^m \\ q^s \end{Bmatrix} = \begin{Bmatrix} F^m \\ 0 \end{Bmatrix} \qquad (19) $$

Thus, the definition of transformation matrices for each type of reduction method is based on the use of the equilibrium equation (Eq. (19)) as will be shown in the following sections.

### 4.1. Dynamic reduction method

This method was proposed by Leung [16] and then by Petersmann [17]. It uses jointly the modal synthesis method and the dynamic reduction method, similar to the substructuring technique proposed by Craig and Bampton [36] in the component modal synthesis context.

Using the second part of rows of Eq. (19), the sub-vector of slave dofs $\{q^s\}$ can be expressed in terms of master dofs $\{q^m\}$ as

$$ \{q^s\} = -([K^{ss}] - \lambda [M^{ss}])^{-1}([K^{sm}] - \lambda [M^{sm}]) \{q^m\} \qquad (20) $$

It should be noted that this expression is defined when the slave dynamic stiffness $Z(\lambda) = ([K^{ss}] - \lambda [M^{ss}])$ is nonsingular. Indeed, this condition is satisfactory for each frequency, else the eigenfrequencies $(\lambda \neq \sigma_v, v = 1, ..., s)$ of the slave problem are defined as follows:

$$ ([K^{ss}] - \sigma_v[M^{ss}])\{\varphi_v\} = 0 \qquad (21) $$

where $\Sigma = \text{diag}(\sigma_v)$ and $[\Phi] = [\cdots \{\varphi_v\} \cdots]$; $v = 1, ..., s$ are the spectral and modal matrices, respectively.

This leads to defining the dynamic contribution of the slave dofs as

$$ [T_d(\lambda)] = -([K^{ss}] - \lambda [M^{ss}])^{-1}([K^{sm}] - \lambda [M^{sm}]) \qquad (22) $$

As can be seen, this relation (22) is an exact dynamic relation which depends strongly on the value of the unknown eigenvalue $\lambda$. This leads to resolving a non-linear eigenvalue problem. Nevertheless, the viscoelastic behavior of the structure is linearized [35], so it is necessary to approximate this relation to be adequate for the linear viscoelastic problem.

Hence, according to Leung and Petersmann method [16,17], the hybrid projection coordinates can be expressed as follows:

$$ \{q^s\} = -[K^{ss}]^{-1}[K^{sm}][q^m] + [\Phi_p]\{c\} \qquad (23) $$

where $[\Phi_p]$ is the $p$-truncated modal basis of the slave structure with $p \ll s$.

This hybrid formulation is similar to the Craig-Bampton method applied in the case of the sub-structuring procedure.
Thus, the master dofs are the junction dofs and the slave dofs are the interior dofs.

While this dynamic reduction formulation was developed for undamped systems, a straightforward application of the above developments to the viscoelastic damped sandwich structures yields the following expression of the slave displacement vector:

$$ \{q^s\} = -[K_q^{ss}]^{-1}[K_q^{sm}]\{q^m\} - [K_q^{ss}]^{-1}[K_{qz}^{sm} \quad K_{qz}^{ss}]\{z\} + [\Phi_p]\{c\} \quad (24) $$

Hence, the reduction of the full dofs to the reduced dofs is achieved as follows:

$$ \begin{Bmatrix} q^m \\ q^s \\ z \end{Bmatrix} = \underbrace{\begin{bmatrix} I_1 & 0 & 0 \\ t_1 & t_2 & \Phi_p \\ 0 & I_2 & 0 \end{bmatrix}}_{[T_{Dyn}]} \begin{Bmatrix} q^m \\ z \\ c \end{Bmatrix} \quad (25) $$

where $[T_{Dyn}]$ is the dynamic transformation matrix; $[I_1]$ and $[I_2]$ are the identity matrices of appropriate size; $[t_1] = -[K_q^{ss}]^{-1}[K_q^{sm}]$ and $[t_2] = -[K_q^{ss}]^{-1}[K_{qz}^{sm} \quad K_{qz}^{ss}]$ represent the static contribution of the structure and $\Phi_p$ basis represents the dynamic contribution of the structure.

For the damped GHM model, the slave problem can be written as follows:

$$ ([K_q^{ss}] - \sigma_v[M^{ss}])\{\phi_v\} = 0 \quad (26) $$

The reduced mass, stiffness and damping matrices can then be written in the form of Eq. (17) using Eq. (25).

This method can predict well the dynamic behavior of viscoelastic sandwich structures combining static and dynamic contributions through hybrid reduced coordinates. Nevertheless, it requires the computation of $p$ truncated modes which increase the size order of the system and leads to less additional CPU time.

## 4.2. Guyan reduction method

The simplest, yet very useful, model reduction method is the Guyan reduction method, which was introduced by Guyan [18] and Iron [19] in 1965. This method is an approximation of the dynamic reduction method, according to which the inertia associated with the slave coordinates is neglected. Thus, applying this static reduction procedure for the damped sandwich structures, the relationship between the full dofs and the reduced dofs can be expressed as follows:

$$ \begin{Bmatrix} q^m \\ q^s \\ z \end{Bmatrix} = \underbrace{\begin{bmatrix} I_1 & 0 \\ t_1 & t_2 \\ 0 & I_2 \end{bmatrix}}_{[T_{St}]} \begin{Bmatrix} q^m \\ z \end{Bmatrix} \quad (27) $$

where $[T_{St}]$ is the static (or Guyan) transformation matrix $[I_1]$; $[I_2]$; $[t_1]$ and $[t_2]$ are the same quantities as defined for Eq. (25).

Under this form, it appears that Guyan reduction is a particular case of dynamic reduction when no slave modes are taken into account.

The reduced mass, stiffness and damping matrices can then be written in the form of Eq. (17) using Eq. (27).

### 4.2.1. Validity domain

This method is valid and useful in an accurate domain. This domain is limited by the cutoff frequency [24,25]. It is the smallest frequency determined by the resolution of the eigenvalue problem (26) defined as $f_1 = f_c$. Thereby, in the practice of applications, the validity domain of the Guyan reduction method is $[0 : f_c/3]$, which reflects the “effective” frequency band.

Consequently, the quality of Guyan approximation depends on the good selection of master dofs which defines its validity domain. In practice, an optimal selection of master dofs must be based on the maximization of the cutoff frequency $f_c$. Out from the validity domain of this method, the accuracy of the obtained results is not well controlled.

### 4.3. Modal reduction method

Modal reduction consists in the derivation of $[\Lambda]$ spectral matrix and $[Q]$ modal basis corresponding to the eigenvalues and eigenvectors of the associated undamped system Eq. (4). Then, these matrices are divided into two parts as follows:

$$ [Q] = [Q_1 \quad Q_2]; \quad [\Lambda] = \begin{bmatrix} \Lambda_1 & 0 \\ 0 & \Lambda_2 \end{bmatrix} \quad (28) $$

The displacement vector is presented as a combination of the $p$ first eigenvectors contained in $[Q_1] \in R^{N \times p}$:

$$ \{q\} = [Q_1]\{c\} \quad (29) $$

Applying this procedure to damped viscoelastic sandwich structures, the full model can be reduced through the projection in the generalized space as follows:

$$ \begin{Bmatrix} q \\ z \end{Bmatrix} = \underbrace{\begin{bmatrix} Q_1 & 0 \\ 0 & I \end{bmatrix}}_{[T_{Mod}]} \begin{Bmatrix} c \\ z \end{Bmatrix} \quad (30) $$

where $[T_{Mod}]$ is the modal transformation matrix.

The reduced mass, stiffness and damping matrices can then be written in the form of Eq. (17) using Eq. (30).

This method uses non-physical coordinates and the truncation can induce errors in the evaluation of dynamic responses. In practice, the modal reduction method leads to good accuracy results when the $p$ first modes are chosen typically from 1.5 to 3 times the frequency band of interest.

### 4.4. Modal reduction in physical space

This method was proposed by O' Callahan [20]. It is based on the modal projection in the physical coordinates. O' Challahan [21,22] marked that this technique allows, after expansion, to return from the reduced model $p$ exact solutions of the full model.

Indeed, the base $[Q_1]$ is partitioned into $m$ master dofs and $s$ slave dofs as follows:

$$ \{q\} = \begin{Bmatrix} q^m \\ q^s \end{Bmatrix} = \begin{bmatrix} Q_{1m} \\ Q_{1s} \end{bmatrix} \{c\} \quad (31) $$

The first line of Eq. (31) leads to

$$ \{q^m\} = [Q_{1m}]\{c\} \quad \text{with} \quad [Q_{1m}] \in R^{N \times p} \quad (32) $$

According to the used $m$ master dofs and $p$ retained undamped modes, three cases can be highlighted:

(a) $m = p$ and $[Q_{1m}]$ is nonsingular: Eq. (32) can be solved exactly as

$$ \{c\} = [Q_{1m}]^{-1}\{q^m\} \quad (33a) $$

(b) $m < p$:

This case leads to infinity of solutions of $\{c\}$ which is not accurate.
<table><caption>**Table 1**
Mechanical and geometrical characteristics of the viscoelastic sandwich structures.</caption><thead><tr><th>Elastic layer (1)</th><th>Viscoelastic core</th><th>Elastic layer (2)</th></tr></thead><tbody><tr><td>L<sub>beam</sub> = L<sub>plate</sub> = 500 mm</td><td>L<sub>beam</sub> = L<sub>plate</sub> = 500 mm</td><td>L<sub>beam</sub> = L<sub>plate</sub> = 500 mm</td></tr><tr><td>b<sub>beam</sub> = 38 mm/b<sub>plate</sub> = 400 mm</td><td>b<sub>beam</sub> = 38 mm/b<sub>plate</sub> = 400 mm</td><td>b<sub>beam</sub> = 38 mm/b<sub>plate</sub> = 400 mm</td></tr><tr><td>h<sub>f<sub>1</sub></sub> = 4.5 mm</td><td>h<sub>c</sub> = 0.2 mm</td><td>h<sub>f<sub>2</sub></sub> = 0.5 mm</td></tr><tr><td>G<sub>f<sub>1</sub></sub> = 70.3 × 10<sup>9</sup> N/m<sup>2</sup></td><td>G<sub>c</sub> (GHM modulus)</td><td>G<sub>f<sub>2</sub></sub> = 70.3 × 10<sup>9</sup> N/m<sup>2</sup></td></tr><tr><td>ρ<sub>f<sub>1</sub></sub> = 2750 kg/m<sup>3</sup></td><td>ρ<sub>c</sub> = 1099.5 kg/m<sup>3</sup></td><td>ρ<sub>f<sub>2</sub></sub> = 2750 kg/m<sup>3</sup></td></tr><tr><td>ν<sub>f<sub>1</sub></sub> = 0.3</td><td>ν<sub>c</sub> = 0.5</td><td>ν<sub>f<sub>2</sub></sub> = 0.3</td></tr></tbody></table>

<table><caption>**Table 2**
Parameters of the GHM viscoelastic model identified for material 242F01 3M™ for one mini-oscillator.</caption><thead><tr><th>Mini-oscillator (i=1)</th><th>Value</th></tr></thead><tbody><tr><td>α<sub>i</sub></td><td>1.047</td></tr><tr><td>ψ<sub>i</sub></td><td>3911.89</td></tr><tr><td>ω<sub>i</sub> [rad/s]</td><td>4943.06</td></tr><tr><td>G<sub>0</sub> [MPa]</td><td>0.079</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Clamped-free (C-F) sandwich beam finite element.</figcaption></figure>

(c) $m > p$ and $(rank([Q_{1m}]) = p)$ is maximal:
Eq. (32) can be solved in the sense of linear least squares as
$$ \{c\} = [Q_{1m}]^+ \{q^m\} \qquad (33b) $$
where $[Q_{1m}]^+ = ([Q_{1m}]^T [Q_{1m}])^{-1} [Q_{1m}]^T$ is the Moore Penrose pseudo inverse.

Thus, by substituting Eq. (33b) into Eq. (31), the structural vector dofs can be written as

$$ \{q\} = \begin{bmatrix} Q_{1m} & Q_{1m}^{+} \\ Q_{1s} & Q_{1m}^{+} \end{bmatrix} \{q^m\} \qquad (34) $$

Hence, the relationship between full and reduced dofs through the projection in physical coordinates for viscoelastic sandwich structures can be expressed generally when $m \ge p$ as follows:

$$ \left\{ \begin{matrix} q^m \\ q^s \\ z \end{matrix} \right\} = \underbrace{\begin{bmatrix} Q_{1m} & Q_{1m}^+ & 0 \\ Q_{1s} & Q_{1m}^+ & 0 \\ 0 & 0 & I \end{bmatrix}}_{[T_{SP}]} \left\{ \begin{matrix} q^m \\ z \end{matrix} \right\} \qquad (35) $$

where $[T_{SP}]$ is the modal transformation matrix defined in physical space.

Consequently, the reduced mass, stiffness and damping matrices can then be written in the form of Eq. (17) using Eq. (35). The definition of a transformation matrix using a number of master dofs either equal to or up the number of undamped modes (as mentioned in the case a or c) leads to the use of maximum rank sub-basis. In fact, the rank of a matrix is defined by the number of rows or columns linearly independent. Numerically, this linear independence is evaluated by the conditioning number.

O' Challahan [20–22] and Friswell [23] shows in their previous studies, that using a sub-basis $[Q_{1m}]$ with high conditioning

number ($> 10^7$; $10^8$) can affect the accuracy of results and can generate erroneous responses. So, the process of selection of master dofs is achieved such that $[Q_{1m}]$ has the minimum conditioning number.

The modal reduction in physical space method allows the derivation of a reduced model which is a subset of the original model expressed in physical coordinates. Furthermore, this technique provides an expanded choice of master dofs but it remains limited by the minimum conditioning condition.

## 5. Numerical applications

In this section, numerical applications are presented in order to illustrate the finite element procedure used for viscoelastic sandwich beam and plate models and outline the practice interest of the proposed reduction strategies. Hence, we consider one mini-oscillator ($N_G=1$) of viscoelastic sandwich beam and plate which are constituted by two elastic layers (faces) in aluminum and a viscoelastic layer (core) of the nuance 242F01. The material and geometrical characteristics of the used sandwich structures are shown in Table 1.

The values of the parameters of the viscoelastic commercially available 242F01, manufactured by 3M™ used at 25 °C for one mini-oscillator, are presented in Table 2 [29].

### 5.1. Viscoelastic sandwich beam

The used FE mesh for the viscoelastic sandwich beam involves 2 elements through the width and 20 elements along the length, having a total number of 1600 dofs. The excitation point is selected in the extremity of the beam (point P, dof of translation $u_z$) and the responses are depicted in two different points P and K (dof of translation $u_z$) as shown in Fig. 3.

In the remainder of this section, the results derived from the implementation of the GHM model, as well as the responses of reduced models used the reduction techniques described as above will be presented both in frequency and in time domains.

#### 5.1.1. Frequency domain evaluation

The interest herein is focused on the frequency-domain responses for both full and reduced GHM models.

The full GHM model response can be derived directly using Eq. (11) in order to calculate the frequency response function (FRF) matrix as follows:

$$ H(\omega) = \{a\}[Z(\omega)]^{-1}\{b\}^T \qquad (36) $$

where $[Z(\omega)] = -\omega^2[M_G] + j\omega[D_G] + [K_G]$ is the dynamic stiffness matrix associated with the damped viscoelastic structure; $\{b\}^T$ is a column vector which defines among all discretized dofs of the structure the position of the selected excitation degree of freedom; $\{a\}$ is a row vector containing the coordinates where the responses are taken into account.
<figure><img src="image_5.png" /><figcaption>Fig. 4. FRFs of the full GHM model plotted in (a) point P and (b) point K.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. FRFs of full model derived from two mathematical tools.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. FRFs of the damped and the undamped GHM model for the sandwich beam.</figcaption></figure>

<table><caption>**Table 3**
Undamped and damped eigenfrequencies of the sandwich beam.</caption><thead><tr><th>Frequency</th><th>Undamped eigenfrequencies [Hz]</th><th>Damped eigenfrequencies [Hz]</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>18.17</td><td>14.65</td></tr><tr><td>f<sub>2</sub></td><td>113.83</td><td>91.75</td></tr><tr><td>f<sub>3</sub></td><td>318.62</td><td>256.92</td></tr><tr><td>f<sub>4</sub></td><td>420.80</td><td>343.94</td></tr><tr><td>f<sub>5</sub></td><td>624.35</td><td>503.72</td></tr></tbody></table>

Fig. 4 depicts the frequency responses of the full GHM model plotted in the points P and K in the frequency band of interest [0–700] Hz.

The frequency responses represented by Fig. 4(a) and (b) are considered as the reference for the full GHM model. Indeed, the FRF amplitudes in [dB] have been computed using a convenient reference factor through the relation Amplitude[dB] = 20 log₁₀(|H(ω)|/1e⁻⁶).

These responses can be determined by the resolution of Eq. (8), which describes the shear modulus as a rational function as well as by Eq. (11) which derives the problem as a second order differential equation. In fact, the mathematical development established in order to derive a second order equation of motion (Eq. (12)) has interest in time domain analysis while the frequency analysis can be carried out directly by the resolution of Eq. (8).

Fig. 5 represents the FRFs derived from the resolution of Eqs. (8) and (11) plotted for the point P.

It can be observed that the two frequency responses derived from the resolution of Eqs. (8) and (11) are perfectly identical. This ensures the equivalence of the two equations.

The next step consists in the determination of damped and undamped frequencies of the sandwich beam which is carried out by the resolution of the eigenvalue problem associated with the damped Eq. (11) and the undamped Eq. (4) systems, respectively. Table 3 represents the five first undamped and damped frequencies of the viscoelastic sandwich beam.

Fig. 6 shows the frequency responses corresponding to the damped and undamped systems plotted for point P. This will illustrate the effect of viscoelastic damping.

5.1.1.1. *Elimination dofs reduction approach.* Now, we will compare the dynamic and Guyan reduction methods belonging to the elimination dofs approach with the full model for the viscoelastic sandwich beam.

At the beginning, we choose $m=30$ master dofs for both reduction methods such that we obtain two transformation matrices having the same dimension: $Dim([T_{Dyn}]) = Dim([T_{St}])$ and the size of each transformation is equal to (1600 × 830). For the Guyan reduction method, the choice of $m=30$ master dofs is constituted by the $u_z$ dofs which are the translation dofs normal to the midplane of the sandwich beam. Furthermore, this choice is carried out maximizing the cutoff frequency. Then, the reduction process is applied and the reduced mass, stiffness and damping
<figure><img src="image_9.png" /><figcaption>Fig. 7. FRFs for full and reduced models: dynamic/Guyan reduction methods: (a) point P and (b) point K of the viscoelastic sandwich beam.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Table 4**
Full and Guyan reduced eigenfrequencies for the viscoelastic sandwich beam.</figcaption></figure>

<table><thead><tr><th>Frequency</th><th>Full<br>frequencies [Hz]</th><th>Guyan reduced<br>frequencies [Hz]</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>14.65</td><td>14.65</td></tr><tr><td>f<sub>2</sub></td><td>91.75</td><td>91.73</td></tr><tr><td>f<sub>3</sub></td><td>256.92</td><td>283.94</td></tr><tr><td>f<sub>4</sub></td><td>343.94</td><td>354.48</td></tr><tr><td>f<sub>5</sub></td><td>503.72</td><td>563.75</td></tr></tbody></table>

matrices as well as the external load vector corresponding to each reduction method are calculated. The reduced dynamic stiffness is evaluated as follows: $[Z_c(\omega)] = -\omega^2[M_c] + j\omega[D_c] + [K_c]$ and the reduced model corresponding to each reduction method compared to the full model is presented. Fig. 7(a) and (b) shows the corresponding FRFs plotted for points P and K, respectively.

From both Fig. 7(a) and (b), the frequency response for the dynamic reduced model is identical to that of the full model while the Guyan reduced model agreed well with the full model for the two first modes of vibration but deviated after the cutoff frequency ($f_c = 165$ Hz) which defined the validity domain of the Guyan reduction method. After this cutoff frequency, the reduced model does not agree with the full model but follows its shape curve. This leads to the conclusion that the Guyan reduction method has the capacity to reproduce the original system but it remains limited by its validity domain. Nevertheless, the dynamic reduction method enriched by the first slave modes ($s=10$) provides a very satisfactory agreement with the full model, making it a suitable method for prediction of the dynamic behavior of viscoelastic sandwich structures.

Table 4 shows the frequency values corresponding to the full and the Guyan reduced model. This affirms that the reduced model derived from the Guyan reduction method is able to reproduce only the two first modes of the full model and outlines that the validity domain of this method is limited by the cutoff frequency.

### 5.1.1.2. Modal reduction approach.
Here, the frequency responses derived from modal and modal in physical space reduction methods are compared to the full model. Indeed, we determine first the number of modes associated with the undamped structure, which covers 1.5 the frequency band of interest ($1.5f_u = 1100$ Hz) and ($p=17$ modes). Then a projection on the physical coordinates is achieved by the partition of the modal basis into master and slave dofs where

$m=p=17$ (case a in Section 4.4) and we inverse directly the modal basis corresponding to the master dofs contribution such that $[Q_{1m}]$ has the minimum conditioning number. Let $cond([Q_{1m}]) = 50$. Thus, we obtain two transformation matrices having the same size $Dim([T_{Mod}]) = Dim([T_{SP}])$ with the dimension of each basis being equal to $(1600 \times 817)$.

Fig. 8(a) and (b) depicts the frequency responses for full and reduced models.

As expected, the frequency responses curves for reduced model and full model are in good agreement for the points P and K as shown in Fig. 8(a) and (b). This leads to the conclusion that modal reduction projected in generalized or in physical coordinates is a viable method for the prediction of the dynamic behavior of structures incorporating viscoelastic materials.

### 5.1.1.3. Physical coordinates approach.
An overlap of elimination dofs reduction approach and modal reduction approach is realized through the projection of the physical coordinates leading to a comparison of the Guyan reduction method to modal reduction in physical space with the full model. Hence, the two transformation matrices must have the same size to compare them such that the master dofs of the Guyan reduction method is equal to the master dofs of modal reduction in the physical space method. Let the number of master dofs $m=30$ which is higher than the number of modes $p=17$. Consequently, we test the case where $m > p$ (case c in Section 4.4) for modal reduction in physical space method. This leads to two transformation matrices as $Dim([T_{SP}]) = Dim([T_{St}])$ with size $(1600 \times 830)$.

Fig. 9(a) and (b) represents the corresponding frequency responses comparison.

As can be seen, the frequency response derived from the modal reduction method is in good agreement with the reference while the frequency response of the Guyan reduction method deviates after the cutoff frequency in both Fig. 9(a) and (b). This affirms that the Guyan reduction method is limited by its validity domain but it can generally predict the dynamic behavior of viscoelastic structures with less accuracy than modal reduction in physical space method, which needs an additional time of evaluation relative to the Guyan method.

Consequently, through the projection on physical coordinates both Guyan and modal reduction in physical space are viable methods which can reproduce the original model. However, for the Guyan reduction method the optimal choice of master dofs is conditioned by the maximum of cutoff frequency.
<figure><img src="image_11.png" /><figcaption>Fig. 8. FRFs for full and reduced models: Modal/Modal in physical space reduction methods: (a) point P and (b) point K of the viscoelastic sandwich beam.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. FRFs for full and reduced models: Guyan/Modal in physical space reduction methods: (a) point P and (b) point K of the viscoelastic sandwich beam.</figcaption></figure>

<table><caption>**Table 5**
Performance of the proposed reduction methods in frequency domain.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Total CPU time [min]</th></tr><tr><th>Full</th><th>Guyan</th><th>Dynamic</th><th>Modal</th><th>Modal in physical space</th></tr></thead><tbody><tr><th>Reduction ratio (%)</th><td>263 -</td><td>35 87</td><td>52 80</td><td>65 75</td><td>78 70</td></tr></tbody></table>

The performance of all proposed reduction methods in terms of CPU time is shown in Table 5.

Table 5 shows the total computing time for full and reduced models. This time, evaluated for each reduction method, includes the calculations of the transformation matrix and the FRF response which is obtained by solving linear equations at each frequency point. For the clarity of comparison, it should be mentioned that all reduced models have the same size 830. Hence, as can be remarked, while the reduction ratio in terms of models dimension is about 50%, it is so advantageous in terms of CPU time. In fact, this ratio reaches 87% with the Guyan reduction method and 80% with the dynamic reduction method while it does not exceed 75% with the modal reduction method and 70% with modal reduction in physical space. The additional CPU time for the modal reduction method in generalized or physical space can be explained by the calculation requirements of undamped modes and the verification of minimum conditioning number in the case of projection on

physical space. Hence, these reduction methods allow generally a drastic reduction, making them a suitable choice to handle both the prohibitive computational effort and the viscoelasticity, especially for complex structures with a large finite element model or in optimization procedure when the dynamic calculations of such models become more complicated. Consequently, the application of these direct reduction methods in the frequency domain can save time considerably, leading to the applicability and the efficiency of these methods in time domain.

### 5.1.2. Time domain evaluation

The interest here is intended to time domain analysis for the viscoelastic sandwich structures. In fact, the prediction of the dynamic behavior of such structures remains until now focused on the frequency analysis more than time analysis. Here both steady state and transient analysis are carried out.

The resolution of temporal equation of motion Eq. (13) is performed using Newmark's integration technique [37] with an unconditionally stable scheme. This technique is used in order to derive the time responses for both full and reduced models which will be compared for each reduction method. These comparisons are performed through static tools called temporal prediction indicators.

#### 5.1.2.1. Temporal prediction indicators
Results comparison tools are based on the statistic indicator calculations associated with the full
<figure><img src="image_13.png" /><figcaption>Fig. 10. Steady state responses for full and reduced models: dynamic/Guyan reduction methods of viscoelastic sandwich beam: (a) ($f < f_c$) and (b) ($f > f_c$).</figcaption></figure>

<table><caption>**Table 6.1**
Temporal moments for the steady state responses of the viscoelastic sandwich beam ($f < f_c$).</caption><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>9.3002E-6</td><td>1.3804</td><td>0.8220</td></tr><tr><td>Dynamic</td><td>9.3000E-6</td><td>1.3804</td><td>0.8220</td></tr><tr><td>Guyan</td><td>9.3586E-6</td><td>1.3808</td><td>0.8218</td></tr></tbody></table>

<table><caption>**Table 6.2**
Temporal moments for the steady state responses of the viscoelastic sandwich beam ($f > f_c$).</caption><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>4.65001E-6</td><td>0.6902</td><td>0.4110</td></tr><tr><td>Dynamic</td><td>4.65000E-6</td><td>0.6902</td><td>0.4110</td></tr><tr><td>Guyan</td><td>4.58005E-6</td><td>0.6700</td><td>0.3889</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Steady state responses for full and reduced models: Modal/Modal in physical space reduction methods of the viscoelastic sandwich beam.</figcaption></figure>

and reduced responses. In fact temporal moments are usually used to quantify a temporal signal in order to compare several transient responses [38]. Hence, the ith order of the temporal moment of a

<figure><img src="image_15.png" /><figcaption>**Table 7**
Temporal moments for the steady state responses of modal reduction approach of the viscoelastic sandwich beam.</figcaption></figure>

<table><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>9.3002E-6</td><td>1.3804</td><td>0.8220</td></tr><tr><td>Modal</td><td>9.2904E-6</td><td>1.3804</td><td>0.8220</td></tr><tr><td>Modal in physical space</td><td>9.2904E-6</td><td>1.3804</td><td>0.8220</td></tr></tbody></table>

response $y(t)$ is defined as [39]

$$M_i = \int_{-\infty}^{+\infty} (t - t_s)^i (y(t))^2 dt \quad (37)$$

where $t_s$ represents the temporal shift and i the moment index order.

In this case, the temporal moment $M_i$ is defined for $t_s = 0$ and normalized as follows:

$$ \begin{cases} E = M_0, & \text{Energy (ms}^2) \\ T = \frac{M_1}{M_0}, & \text{Central time (centroid) (s)} \\ D^2 = \frac{M_2}{M_0} - \left(\frac{M_1}{M_0}\right)^2, & \text{Root means square duration (s)} \end{cases} \quad (38) $$

Thereby, this triplet of indicators ($E$, $T$ and $D$) enables one to determine the error generated both in the amplitude and in the time scales. Indeed, $E$ is used to identify the error in the amplitude of the response; $T$ and $D$ are used to identify the error in the periodicity of the response.

Gerges [40] proved that a relative error of order of ±4% in energy $E$, ±2% in central time $T$ and ±4% in root mean square $D$ are admissible in order to validate the reduced model compared to the full model.

In the remainder and for good clearance, it should be noticed that all time responses will be plotted only on point P.

5.1.2.2. Steady state analysis. The sandwich beam is subjected to an harmonic load of the form ($F(t)$) = $F_0 \sin(\omega t)$ where $F_0 = 1N$ and $\omega = 50\pi$ rad/s; ($f = 25$ Hz). The steady state response is established after 1 s of transient response and the oscillations are well stable over a period of time from 2 to 3 s. Therefore, in the following, the time response derived from a harmonic excitation will be presented in the interval of time [2–3] s.

The same procedure of comparison for reduced and full models presented in Section 5.1.1 is also carried out in this section with the time domain analysis. For each type, the size of the reduced models is kept the same as mentioned in the previous sections.
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Transient responses for full and reduced models: dynamic/Guyan reduction methods of the viscoelastic sandwich beam.</figcaption></figure>

<table><caption>**Table 8**
Temporal moments for the transient responses of eliminated dofs reduction approach of the viscoelastic sandwich beam.</caption><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>1.0387E-6</td><td>0.2055</td><td>0.0498</td></tr><tr><td>Dynamic</td><td>1.0387E-6</td><td>0.2055</td><td>0.0498</td></tr><tr><td>Guyan</td><td>1.0403E-6</td><td>0.2054</td><td>0.0498</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Transient responses for full and reduced models: Modal/Modal in physical space reduction methods of the viscoelastic beam.</figcaption></figure>

<table><caption>**Table 9**
Temporal moments for the transient responses of modal reduction approach of the viscoelastic sandwich beam.</caption><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>1.0387E-6</td><td>0.2055</td><td>0.0498</td></tr><tr><td>Modal</td><td>1.0369E-6</td><td>0.2055</td><td>0.0498</td></tr><tr><td>Modal in physical space</td><td>1.0369E-6</td><td>0.2055</td><td>0.0498</td></tr></tbody></table>

5.1.2.3. Steady state responses for elimination dofs reduction approach. The time responses for full and reduced models are presented in Fig. 10(a) and (b).

Fig. 10(a) shows that the time responses curves for the full model and those of reduced models are in good correlation over a

<table><caption>**Table 10**
Performance of the proposed reduction methods in time domain.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Total CPU time [min]</th></tr><tr><th>Full</th><th>Guyan</th><th>Dynamic</th><th>Modal</th><th>Modal in physical space</th></tr></thead><tbody><tr><th>Reduction ratio (%)</th><td>365</td><td>40</td><td>65</td><td>85</td><td>92</td></tr><tr><th></th><td>-</td><td>88</td><td>82</td><td>77</td><td>75</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 14.** FE model for the viscoelastic sandwich plate with the position of the optimal master dofs (●).</figcaption></figure>

period of time of 1 s. Here, the Guyan reduction method agrees well with the full model because the excitation frequency in which the sandwich beam is subjected is less than the cutoff frequency ($f < f_c$). So, the excitation covers the validity domain of this method, leading to good agreement with the reference.

Furthermore, Table 6.1 shows that the relative error in energy $E$ between the full model and the Guyan reduced model is of order of -0.62% while this error does not exceed 0.002% with the dynamic reduced model. This confirms the visual impression in amplitudes.

In addition, the relative error in central moment $T$ and means square root $D$ does not exceed 0.04% for the Guyan reduced model and it is practically equal to zero for the dynamic reduced model. Consequently, the obtained results in the case where ($f < f_c$) presents a satisfactory accuracy compared to the full model, thus enabling a validation of the reduced models.

However, when the excitation frequency is higher than the cutoff frequency ($f > f_c$), the obtained results for full and reduced models of viscoelastic sandwich beam subjected to harmonic load under a frequency excitation $f = 300$ Hz are depicted in Fig. 10(b).

As can be seen, in this case, the results start to lose their accuracy. Indeed, the reduced Guyan response presents an apparent deviation in both amplitudes and time scales (Table 6.2). The deviation in amplitudes scale is indicated by a relative error which reaches 1.5% in energy $E$. In the time scale, the relative error has the order of 5% in $D$ and 3% in $T$. These values are significant compared to the case where ($f < f_c$) and leads one to conclude that the Guyan reduction method is limited by its validity domain. Hence, beyond the cutoff frequency, the Guyan reduction method is less accurate. Nevertheless, the dynamic reduction method preserves its capacity to reproduce the full steady state response in both cases, leading to affirm the performance of this reduction method in the prediction of the dynamic behavior of viscoelastic sandwich structures.

5.1.2.4. Steady state results for modal reduction approach. For modal reduction methods the obtained results are shown in Fig. 11.
<table><caption>**Table 11**
Undamped and damped frequencies for the viscoelastic sandwich plate.</caption><thead><tr><th>Frequency</th><th>Undamped frequencies [Hz]</th><th>Damped frequencies [Hz]</th></tr></thead><tbody><tr><td>f<sub>1</sub></td><td>264.26</td><td>213.07</td></tr><tr><td>f<sub>2</sub></td><td>463.78</td><td>374.13</td></tr><tr><td>f<sub>3</sub></td><td>604.83</td><td>488.05</td></tr><tr><td>f<sub>4</sub></td><td>788.08</td><td>636.37</td></tr><tr><td>f<sub>5</sub></td><td>1095.50</td><td>885.47</td></tr><tr><td>f<sub>6</sub></td><td>1125.20</td><td>908.92</td></tr><tr><td>f<sub>7</sub></td><td>1226.53</td><td>991.01</td></tr><tr><td>f<sub>8</sub></td><td>1300.83</td><td>1051.70</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** FRFs for full and reduced models: Guyan/dynamic reduction methods for the viscoelastic sandwich plate.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** FRFs for full and reduced models: Modal/Modal reduction in physical space methods for the viscoelastic sandwich plate.</figcaption></figure>

It can be observed in Fig. 11 that the time responses of steady state motion for the reduced models are in good agreement with the full model. This is confirmed by the values of the three central temporal moments. In fact, Table 7 shows that the relative error in *E* does not exceed 0.1% while the relative error in *T* and *D* is practically zeros. This leads to validation of the reduced models which allows a perfect reproduction of the original model in amplitude and time scales.

The comparison of the Guyan reduction method with the modal reduction method in physical space for steady state responses is also carried out. The obtained results indicate a relative error in energy *E* of 0.1% for modal reduction in physical space which reaches 0.6% for the Guyan reduction method. This implies that through the projection on physical space, the modal reduction method in physical space has the capacity to reproduce the original model better than the Guyan reduction method.

5.1.2.5. *Transient analysis.* In this section, the viscoelastic beam is subjected to an impulse load of duration $T_{impulse} = 2$ ms and amplitude equal to 10 N. The same strategy of comparison between the different reduction methods is carried out. Furthermore, as mentioned in the previous sections, for each reduction approach the equality of reduction basis is provided.

5.1.2.6. *Transient results for elimination dofs reduction approach.* The comparison between reduced models and full model is shown in Fig. 12.

The time responses to an impulse excitation at the point P of the viscoelastic sandwich beam are well correlated before and after reduction. Besides, the three central moments reflect that the dynamic reduction method is a viable method that reproduces entirely the original model. Since the frequency spectrum of the impulse excitation covers the validity domain of the Guyan reduction method, the reduced response derived from this method agrees well with the original while it represents a relative error in *E* of the order of 0.15% and a relative error in *T* of the order of 0.04% (Table 8). Thus, these values can validate the Guyan reduced model. As a result, this method presents a suitable choice in terms of simplicity, feasibility of implementation and also satisfactory accurate results.

5.1.2.7. *Transient results for modal reduction approach.* The derived transient results for both reduced and full models are presented in Fig. 13.

Fig. 13 shows the transient responses for the reduced models derived from modal and modal in physical space reduction methods compared to the full model. It can be observed that these responses are identical. In fact, the modal reduction method returns the *p* first exact modes of the associated undamped model, allowing a reproduction of the original model through a generalized coordinate projection while modal reduction in physical space method allows a reproduction of the full model through a projection in physical coordinates. This is affirmed by the three central moments (Table 9), which indicates that both reduced models preserve the periodicity of the full response with a relative error in the energy *E*, which does not exceed 0.17%, leading to the validation of these two reduction methods in temporal domain.

Table 10 presents the performance of the proposed reduced models compared to the full model in time domain. There is a significant reduction ratio in total CPU time required for the evaluation of reduced basis and temporal responses at each iteration which justifies the efficiency of these reduction methods in time domain.

## 5.2. Viscoelastic sandwich plate

In this example, the interest is focused on the validity domain extension of the Guyan reduction method. In fact, after meshing the plate into 20 × 15 finite elements as shown in Fig. 14, the master dofs, which are translation dofs $u_z$, are chosen such that the cutoff frequency is maximal. Hence, the distribution of the chosen master dofs ($m=40$) is illustrated in Fig. 14.

The plate is clamped on the four sides (C-C-C-C). The FE discretization scheme leads to 8310 dofs in total number. The excitation and the responses are depicted in point *E* as presented in Fig. 14.

### 5.2.1. Frequency domain analysis

The frequency analyses for the sandwich plate are carried out in the same manner as mentioned in the previous sections. In fact, for each type of reduction, the equality of bases is assured.
<figure><img src="image_21.png" /><figcaption>Fig. 17. Steady state response for the full and reduced models: Guyan/dynamic methods of the viscoelastic sandwich plate: (a) ($f < f_c$) and (b) ($f > f_c$).</figcaption></figure>

<caption>**Table 12.1**</caption>

Temporal moments for the steady state responses ($f < f_c$) of eliminations dofs reduction approach of the viscoelastic sandwich plate.

<table><caption>**Table 12.2**</caption><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>1.8377E-8</td><td>0.2428</td><td>0.0218</td></tr><tr><td>Dynamic</td><td>1.8377E-8</td><td>0.2428</td><td>0.0218</td></tr><tr><td>Guyan</td><td>1.8516E-8</td><td>0.2428</td><td>0.0218</td></tr></tbody></table>

Temporal moments for the steady state responses ($f > f_c$) of eliminations dofs reduction approach of the viscoelastic sandwich plate.

<table><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>2.1123E-8</td><td>0.3578</td><td>0.0412</td></tr><tr><td>Dynamic</td><td>2.1123E-8</td><td>0.3578</td><td>0.0412</td></tr><tr><td>Guyan</td><td>2.1128E-8</td><td>0.3578</td><td>0.0412</td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Steady state response for the full and reduced models: Modal/Modal in physical space methods of the viscoelastic sandwich plate.</figcaption></figure>

<caption>**Table 13**</caption>

Temporal moments for the steady state responses of modal reduction approach of the viscoelastic sandwich plate.

<table><thead><tr><th></th><th>E</th><th>T</th><th>D</th></tr></thead><tbody><tr><td>Full</td><td>1.8377E-8</td><td>0.2428</td><td>0.0218</td></tr><tr><td>Modal</td><td>1.8375E-8</td><td>0.2428</td><td>0.0218</td></tr><tr><td>Modal in physical space</td><td>1.8375E-8</td><td>0.2428</td><td>0.0218</td></tr></tbody></table>

<figure><img src="image_23.png" /></figure>

For the elimination dofs reduction approach, Guyan and dynamic basis are constructed such that the two bases have the same size (8310 × 4195). Then, modal and modal in physical space bases belonging to the modal reduction approach are also constructed such that the size of each basis is equal to (8310 × 4180), with $p=25$ modes which covers 1.5 the frequency band of interest [0–1200] Hz ($1.5f_u = 1800$ Hz).

Table 11 presents the eight first damped and undamped frequencies of the viscoelastic sandwich plate. Indeed, the difference between the damped and undamped frequencies values indicates the effect of the viscoelastic damping.

### 5.2.1.1. Elimination dofs reduction approach.
The obtained results for the viscoelastic sandwich plate are presented in Fig. 15.

As can be seen, the frequency response derived from the Guyan reduction method has the capacity here to reproduce the frequency response of the full model for the first four modes. Furthermore, beyond the cutoff frequency which is equal to 806 Hz, the Guyan reduced response follows the shape curve of the full model with small difference. This implies that the Guyan reduction method is a viable method for the prediction of the dynamic behavior of viscoelastic sandwich structures, when the choice of master dofs is optimal. So, more the choice is optimal, more the results are accurate. For the dynamic reduction method, its frequency response is in good agreement with the response of the full model. This affirms the efficiency of this method in the reproduction of the full model dynamics.

### 5.2.1.2. Modal reduction approach.
Fig. 16 shows the frequency responses for the reduced models derived from modal and modal reduction in physical space methods compared to the full model.

It can be observed that the two reduced frequency responses derived from modal and modal reduction in physical space are identical to the frequency response of the full model. This leads to confirm that modal and modal reduction method in physical space have the capacity to reproduce the original coordinates of the sandwich structures through a projection on the generalized coordinates as well as on the physical coordinates with good accuracy.

### 5.2.2. Time domain analysis

In this section, the time responses are focused on the steady state analysis for the viscoelastic sandwich plate in order to show the performance of the Guyan reduction method. Indeed, the viscoelastic sandwich plate is subjected to a harmonic load of amplitude equal to 1 N under an excitation frequency equal to
400 Hz (around the second mode of vibration for the sandwich plate). The steady state response is reached after 0.15 s of transient oscillations and where the oscillation becomes more stable, the time responses are plotted.

For the Guyan reduction method, two cases are tested: the first one examines the steady state responses of the sandwich plate where it is subjected to excitation frequency less than the cutoff frequency and the second shows the other case where the plate is subjected to an excitation frequency higher than the cutoff frequency.

5.2.2.1. Case 1: $f < f_c$. Fig. 17(a) shows the steady state responses of the Guyan and dynamic reduced models compared to the full model. It can be seen that the dynamic response presents a satisfactory agreement with the full model. This is clarified by the values of the three central moments (E, T, D) (Table 12.1), which is identical to those of the full model. Thus dynamic reduction method remains a good choice of reduction methods. For the Guyan reduction method, the steady state response reproduces the original response with a relative error in energy E which does not exceed 0.007% as shown in Table 12.2 while the central moments T and D, indicators of error in periodicity, are identical to those of the full model. Hence, the Guyan reduction method is validated for each frequency excitation less than the cutoff frequency.

5.2.2.2. Case 2: $f > f_c$. When the sandwich plate is subjected to a harmonic excitation frequency higher than the cutoff frequency ($f = 1100$ Hz), the steady state response of the Guyan reduced model presents a little shift relative to the full model while dynamic reduced response preserves its capacity to reproduce the full response (Fig. 17(b)). In fact, the reduced response derived from the Guyan reduction method presents a few relative errors in energy E of order of 0.02% while it preserves the periodicity of the full model as shown in Fig. 17(b). Consequently, the Guyan reduction method can predict with good accuracy the viscoelastic behavior of sandwich structures when the choice of master dofs is optimal. So, compared to the case of sandwich beam where the excitation frequency is high than the cutoff frequency, the Guyan reduction method presents in the case of the plate more satisfactory results.

For the modal reduction approach, the obtained results are presented in Fig. 18.

The reduced responses obtained from modal and modal reduction in physical space methods are identical to the full model. This is affirmed by the values of the three central moments (E, T, D) presented in Table 13. In fact, the relative error for each moment for the three compared responses is practically equal to zero. This leads to the conclusion that modal and modal reduction in physical space are viable methods for predicting the dynamic behavior of viscoelastic sandwich plate.

The transient analysis for the viscoelastic sandwich plate subjected to an impulse excitation is also established. Indeed, the transient reduced responses present a good agreement with the full model for each type of reduction. This can be explained for the Guyan reduction method by the frequency spectrum of the impulse excitation which covers the validity domain of this method. Hence, the optimal choice of master dofs is an important step in all reduction procedures, notably for the Guyan reduction method in order to predict well the dynamic behavior of viscoelastically damped structures.

For this example and for the sake of brevity, only the CPU time evaluated in time domain is illustrated (Table 14).

The dynamic potential of the proposed reduction methods is more highlighted with the viscoelastic plate example. In fact, the

saved time required for calculating the full and reduced models increase by increasing the degrees of freedom. Furthermore, these calculations take into account the evaluation of reduced basis and the iterative procedure generated by the use of the Newmark scheme in time domain for each applied reduction method. Hence, these reduction methods constitute an efficient solution to gain time and to handle large finite elements models with viscoelastic components. On the other hand, these methods are used in the direct reduction context and they improve their efficiency notably in terms of CPU time, leading to performing both frequency and temporal analysis. So, when more than one structure is used and takes into account the non-linear behavior of the most structures, the use of model reduction method in the substructuring context or component mode analysis [36] for viscoelastic sandwich structures appears quite attractive.

## 5.3. Temporal analysis with localized nonlinearities in the substructuring context

In this section, attention is focused on assembled viscoelastic sandwich structures. Indeed, the bolted joints are usually modeled by non-linear elements in the junctions of such structures. Therefore, local nonlinearities are introduced to take into account this effect. However, this is done at the price of generation firstly a large systems dimension induced by viscoelastic components and secondly time consuming due to the resolution scheme which become more complicated with the introduction of local nonlinearities. So, it remains challenging to develop an efficient reduction strategy that can overcome this problem. For that, we propose combining the Guyan reduction method with the modal synthesis method for local non-linear viscoelastic structures in the substructuring context. This is done by the addition of a non-linear term in the equation of motion Eq. (13). In fact, the form of this equation as a standard temporal second-order equation leads to introduction of local nonlinearities in a simple and soft way. Thereby, the obtained temporal non-linear equation of motion can be written as follows:

$$[M_G](\dot{q}_G) + [D_G]\{\dot{q}_G\} + [K_G]\{q_G\} + \{f_{nl}(q_G)\} = \{F_G\} \quad (39)$$

where $\{f_{nl}(q_G)\}$ indicates the added non-linear load, where its ith component can be expressed by the Duffing oscillator as follows:

$$\{f_{nl}(q_G)\}_i = \sum_{j=1}^{m} \mu_j [(q_G)_i - (q_G)_j]^3 = [K_{nl}(q_G)]\{q_G\} \quad (40)$$

<table><caption>**Table 14**
CPU time of the viscoelastic sandwich plate.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Total CPU time [min]</th></tr><tr><th>Full</th><th>Guyan</th><th>Dynamic</th><th>Modal</th><th>Modal in physical space</th></tr></thead><tbody><tr><td>Reduction ratio (%)</td><td>1440</td><td>168</td><td>258</td><td>324</td><td>356</td></tr><tr><td></td><td>-</td><td>88</td><td>82</td><td>77</td><td>75</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>**Fig. 19.** FE model of the global non-linear assembled viscoelastic sandwich beam.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 20. Temporal responses of full and reduced models of the non-linear assembled viscoelastic sandwich beam: (a) $f<f_c$ and (b) $f>f_c$.</figcaption></figure>

<table><caption>**Table 15**
Performance of the non-linear assembled viscoelastic sandwich beam.</caption><thead><tr><th rowspan="2"></th><th colspan="2">CPU time [min]</th></tr><tr><th>Full</th><th>Reduced</th></tr></thead><tbody><tr><td></td><td>684</td><td>61</td></tr><tr><td>Reduction ratio (%)</td><td></td><td>92</td></tr></tbody></table>

$m$ represents the number of attached non-linear springs depending on the $i$th dof; $\mu_j$ represents the non-linear stiffness factor for each non-linear spring and $[K_{nl}(q_G)]$ is the non-linear stiffness matrix contribution.

The application of the proposed reduction strategy leads to the following non-linear reduced model:

$$[M_c]\{\ddot{q}_c\} + [D_c]\{\dot{q}_c\} + [K_c]\{q_c\} + [K_{nlc}]\{q_c\} = \{F_c\} \quad (41)$$

where $[M_c]$; $[D_c]$; $[K_c]$; $[K_{nlc}]$ and $\{F_c\}$ represent respectively the reduced mass, damping, linear and non-linear stiffness matrices and the reduced load vector obtained by the application of the Guyan transformation matrix which is described in Section 4.2 with master ($m$) and slave ($s$) dofs expressed for the direct method being replaced, respectively, by junction ($j$) and interior ($i$) dofs for the substructuring procedure.

The FE model of the global viscoelastic sandwich beam is illustrated in Fig. 19.

The used FE model of the global viscoelastic sandwich beam involves 80 elements with 320 nodes and 5 dofs per node, leading to 3200 total dofs. This beam is clamped at its two edges and the mechanical and geometrical properties for each substructure (SS1)

or (SS2) are the same as described for the viscoelastic sandwich beam in Table 1. The value of each used non-linear spring coefficient is $\mu = 10^9$ N/m$^3$.

First, we start from the knowledge of the dynamic behavior of each substructure (SS1) and (SS2) which are reduced separately by the application of the Guyan reduction method.

### 5.3.1. Guyan reduction of substructure (SS1)

The displacement vector $\{q_G\}_{(SS1)}$ of the viscoelastic substructure (SS1) is decomposed according to the junction ($j$) and interior ($i$) dofs partition as follows:

$$\{q_G\}_{(SS1)} = \begin{Bmatrix} q^i \\ z \end{Bmatrix}_{(SS1)} = [T_{St}] \begin{Bmatrix} q^i \\ z \end{Bmatrix} \quad (42)$$

$[T_{St}]$ is the Guyan transformation matrix as defined in Section 4.2. Then, the reduced system is obtained by substituting Eq. (42) into Eq. (17). Thus, the size of the reduced model for the substructure (SS1) is 805 for the studied example. The choice of junction dofs ($j=5$), which are translation dofs $u_z$, is carried out on maximizing the cutoff frequency of the viscoelastic substructure (SS1), which is equal to 165 Hz.

### 5.3.2. Guyan reduction of substructure (SS2)

In the same manner, the displacement vector of the second viscoelastic substructure (SS2) is partitioned in terms of junction
(j) and interior (i) dofs as follows:

$$ \{q_G\}_{(SS2)} = \begin{cases} q^i \\ z \end{cases}_{(SS2)} = [T_{st}] \begin{cases} q^i \\ z \end{cases} \quad (43) $$

The reduced model is obtained in the form of Eq. (17) using Eq. (43). Thereby, its dimension is equal to 805 with $j=5$ dofs. Furthermore, the cutoff frequency of the viscoelastic substructure (SS2) is equal to 165 Hz.

After that, the reduced matrices are assembled taking into account the localized nonlinearities in the junctions between the two viscoelastic substructures (SS1) and (SS2) leading to a global reduced system of order 1610. The obtained temporal results of the global viscoelastic sandwich beam, which is subjected to a harmonic load in point A of amplitude 50 N to arise effectively the non-linear behavior, in terms of displacement and velocity, are presented in Fig. 20(a) and (b).

As can be seen, Fig. 20(a) shows that the temporal responses in terms of displacement and velocity are in good agreement in the case of a harmonic excitation ($f=25$ Hz) less than the cutoff frequency of the non-linear assembled viscoelastic ($f_c=165$ Hz). Furthermore, the evaluation of the three temporal moments ($E, T, D$) provided identical values for both full and reduced models, leading to validate the visual correlation. For the case of high excitation ($f=300$ Hz) relative to the cutoff frequency ($f_c=165$ Hz), the full and reduced models present a shift in amplitude and time scales. This shift is about 3% in energy $E$, 0.1% in $T$ and 0.2% in $D$ for the displacement responses and 5% in $E$, 0.3% in $T$ and 0.1% in $D$ for the velocity responses. This leads to validate the applicability of the proposed method for non-linear viscoelastic structures in time domain. On the other hand, while the reduction ratio in terms of systems order is around 50% for such as a non-linear example, the accuracy of the obtained results in terms of displacement and velocity is satisfactory. Hence, this reduction method presents an efficient tool to handle non-linear structures with viscoelastic materials in time domain, which enables one to perform the frequency analysis with more specific techniques such as the harmonic balance method.

The performance of the proposed method in terms of CPU time is shown in Table 15.

There is a significant CPU reduction ratio of 92%, leading to the conclusion that the proposed reduction method for non-linear viscoelastic sandwich structures enables one to bring two levels: viscoelasticity and nonlinearity for the compromise good accuracy and time efficiency.

It should be mentioned that from the studied examples of viscoelastic sandwich (beam, plate, assembled beams) which are academic structures, the reduction ratio in terms of systems dimension does not exceed 50% but it can be possible to increase this ratio further with more complex structures.

## 6. Conclusions

In this paper, finite element procedures are combined to first-order shear deformation theory (FSDT) and to the GHM model for the modeling of viscoelastic sandwich structures. The introduction of internal variables or dissipation coordinates through a series of mini-oscillators to take into account the viscoelastic damping is achieved. Unfortunately, this was done at the expense of increasing the model order. Consequently, model reduction methods have been proposed as a convenient alternative for this problem. First, the dynamic reduction method based on the elimination of slave dofs and enrichment of the transformation basis with first slave modes is developed. As a result, the reduced model reproduces well the original model with good accuracy and less CPU time,

making it the best choice of model reduction methods for the compromise accuracy-time gain in the direct reduction procedure. Next, the Guyan reduction method is expressed by a static basis, neglecting the inertia associated with the slave coordinates. This method allows a simple implementation in the most finite elements codes with a significant reduction ratio in terms of CPU time and a good capacity of prediction of the original model, especially in the substructuring context where the necessity of an efficient reduction method becomes twice reinforced first by the large systems dimensions induced by viscoelastic components and second by the consuming time generated by the introduction of local nonlinearities. Then, the modal reduction method based on the derivation of the first modes associated with the undamped structure is established. This method constitutes a good representation of the original model with reduced CPU time, making it a suitable choice for the reduction of sandwich structures incorporating viscoelastic materials. Finally, modal reduction in the physical space method is outlined as a projection of modal basis in physical coordinates system. Thereby, the projection on physical space is realized, leading to good results. However, this method needs an additional time compared to other reduction methods and requires verifying the minimum conditioning number condition.

In all reduction procedures, the proposed methods provided good accuracy results and a satisfactory agreement with the full model in both frequency and time domains. On the other hand, even the reduction ratio in terms of systems size did not exceed 50%; it reached 90% in terms of CPU time, which makes these methods a suitable choice in handling viscoelastic sandwich structures in an accurate and efficient way. Furthermore, the kernel of the idea to use model reduction methods in time domain can be explained by the temporal interest of the GHM model which allows transformation from a frequency rational shear modulus function to a temporal resolved second-order equation which improves its importance notably in the substructuring context for structures with local nonlinearities.

## References

[1] J. Salençon, Viscoélasticité, Presse des ponts et chaussés, Paris, 1983.

[2] J.D. Ferry, Viscoelastic Properties of Polymers, John Wiley & Sons, New York, 1980.

[3] F.J. Plantema, Sandwich Construction, John Wiley & Sons, New York, 1966.

[4] E. Reissner, Finite deflections of sandwich plates, AIAA J. 15 (7) (1948) 435-440.

[5] E. Reissner, The effect of transverse shear deformation on the bending elastic plates, J. Appl. Mech. 12 (1945) 69-77.

[6] R.D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic, elastic plates, J. Appl. Mech. 18 (1951) 31-38.

[7] E. Reissner, On the theory of bending of elastic plates, J. Maths Phys. 12 (1944) 184-191.

[8] J.N. Reddy, A simple higher order theory for laminated composites plates, J. Appl. Mech. 51 (1948) 745-752.

[9] A.J.M Ferreira, C.M.C. Roque, M. PALS, Analysis of composite plates using high order shear deformation theory and a finite point formulation based on the multiquadric radial basis function method, Compos. Part B 34 (2003) 627-636.

[10] Y.M. Chugal, R.P. Shimpi, A review of Refined shear deformation theories of isotropic and anisotropic laminated plates, J. Reinf. Plast. Compos. 21 (9) (2002).

[11] J.N. Reddy, Mechanics of Laminated Composites Plates, Theory and Analysis, 2nd ed., CRC Press, Florida, 1997.

[12] M. Di Sciuva, Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: an evaluation of a new displacement model, J. Sound Vib. 105 (1986) 425-442.

[13] D.J. Mc-Tavish, P.C. Hughes, Dynamics of viscoelastic structure—a time domain, finite element formulation, J. Appl. Mech. 52 (1985) 897-906.

[14] D.J. Mc-Tavish, P.C. Hughes, Modeling of linear viscoelastic space structures, J. Appl. Mech. 52 (1993) 103-113.

[15] O.C. Zienkiewicz, R.L. Taylor, TheFinite Element Method, fifth edition, vol. 1: The Basis, 2000, Butterworth-Heinemann, Oxford, Woburn, MA.

[16] Y.T. Leung, An accurate method of dynamic condensation in structural analysis, Int. J. Numer. Methods Eng. 14 (1979) 1241-1256.
[17] N. Petersmann, Calculation of eigenvalues using substructuring and dynamic condensation, in: Proceedings of the Recent International Conference in Recent Advance in Structures Dynamics, Southampton, 1984, pp. 211–219.

[18] R.J. Guyan, Reduction of stiffness and mass matrices, AIAA J. 3 (1965) 380.

[19] B. Iron, Structural eigenvalues problem: elimination of unwanted variables, AIAA J. 3 (1965) 961–962.

[20] J. O'Callahan, A procedure for an improved reduced system (IRS) model, in: Proceedings of the 7th International Modal Analysis Conference, Las Vegas, 1989, pp. 17–21.

[21] J. O'Callahan, System equivalent reduction expansion process (SEREP), in: Proceedings of the 7th International Modal Analysis Conference, Las Vegas, 1989, pp. 27–29.

[22] J. O'Callahan, A non smoothing SEREP process for modal expansion, in: Proceedings of the 12th International Modal Analysis Conference, Honolulu (Hawaii), 1994, pp. 232–238.

[23] M.I. Friswell, S.D. Garvey, J.E.T. Penny, The convergence of the iterated IRS method, J. Sound Vib. 211 (1998) 123–132.

[24] N. Bouhaddi, R. Fillod, Model reduction by a simplified variant of dynamic condensation, J. Sound Vib. 19 (1996) 233–250.

[25] N. Bouhaddi, R. Fillod, A method for selecting master dof in dynamic substructuring using the Guyan condensation method, Comput. Struct. 45 (5/6) (1992) 941–946.

[26] C.H. Park, D.J. Inman, M.J. Lam, Model reduction of viscoelastic finite element models, J. Sound Vib. 219 (1999) 619–637.

[27] M. Meunier, R.A. Shenoi, Forced response of FRP sandwich panels with viscoelastic materials, J. Sound Vib. 263 (2003) 131–151.

[28] M.A. Trindade, A. Benjeddou, R. Ohayon, Modeling of frequency-dependent viscoelastic materials for active-passive vibration damping, J. Vib. Acoust. 122 (2) (2000) 169–174.

[29] A.M.G. De Lima, D.A. Rade, Model Reduction Methods As Applied to Viscoelastically Damped Finite Element Models, Catalao.ufg.br, Brasil, 2010.

[30] C.Y.K. Chee, Static shape control of laminated composite plate smart structure using piezoelectric actuators (Ph.D. Thesis), Sydney University, Department of Aeronautical Engineering, Australia, 2000.

[31] A.W. Faria, Finite element modelling of composites plates: contribution to the damping, damage and incertitudes (Ph.D. thesis), Uberlandia University, Brasil, 2010.

[32] G.A. Lesieutre, Finite element for dynamic modeling of uniaxial rods with frequency dependent material properties, Int. J. Solids Struct. 29 (1992) 1567–1579.

[33] R.L. Bagley, P.J. Torvik, Fractional calculus—a different approach to the analysis of viscoelastically damped structures, AIAA J. 21 (5) (1983) 741–748.

[34] A.M.G. De Lima, A.R. Da Silva, D.A. Rade, N. Bouhaddi, Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped systems, Eng. Struct. 32 (2010) 1479–1488.

[35] R.M. Christensen, The Theory of Linear Viscoelasticity, An Introduction, 2nd ed., Academic Press, New York, 1982.

[36] Jr. R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analysis, AIAA J. 6 (1968) 1913–1919.

[37] N. Newmark, A method of computation for structural dynamics, J. Eng. Mech. Div. 85 (7) (1959) 67–94.

[38] F.M. Hemez, S.M. Doebling, From shock response spectrum to temporal moments and vice-versa, in: Proceedings of the 21st SEM International Modal Analysis Conference, Kissimmee, 2003.

[39] G. Masson, Robust modal synthesis adapted for the optimization of high order models (Ph.D. thesis), Franche Comté University, France, 2003.

[40] Y. Gerges, Model reduction methods for non-linear vibroacoustics (Ph.D. thesis), Franche Comté University, France, 2013.