Forced harmonic response of viscoelastic sandwich beams by a reduction method

ABSTRACT The aim of this article is to develop a reduction method to determine the forced harmonic response of viscoelastic sandwich structures at a reasonable computational cost. The numerical resolution is based on the asymptotic numerical method. This type of problem is complex, and its number of degrees of freedom is double the number of the undamped structure, leading to a high computational time. To address the problem, three reduction methods are evaluated, which differ from the projection operator. Numerical tests have been performed in the case of cantilever sandwich beams. The comparison of the results obtained by the reduction order resolution with those given by the full order resolution shows both a good agreement and a significant reduction in computational cost.


Introduction
Viscoelastic sandwich materials are often used as a passive solution to reduce vibration and noise in many domains, such as automobile, aerospace, civil construction, etc. The dynamic analysis of this type of structure subjected to a harmonic load excitation is an old topic of research. It involves complex stiffness matrices. The complex term is principally due to the complex Young's modulus [1] of the viscoelastic core. Hence, a great deal of literature exists: Lu et al. [2] presented and discussed theoretical and experimental results for the transverse driving point mechanical impedances of damped thin viscoelastic sandwich plates. The dependence on frequency and temperature of the dynamic properties of the viscoelastic materials is taken into consideration. Babert et al. [3] presented a finite element model for the harmonic response of sandwich beams with thin or moderately thick viscoelastic cores. Nonlinear variation of displacements through the thickness of the core is assumed. Hazard and Bouillard [4] studied the passive damping of structural vibrations using viscoelastic layers. The sandwich modeling technique is based on the use of an interface element and on the first-order shear deformation theory. This element couples the lower and upper layers without additional degrees of freedom. Abdoun et al. [5] presented an asymptotic numerical method for forced harmonic vibration analyses of viscoelastic structures. A mathematical formulation that may account for various viscoelastic models is presented. Power series expansions and Padé approximants of the displacement and frequency are developed and the finite element method is used for numerical solutions. The method is applied to harmonic responses of viscoelastic structures with constant and frequency dependent coeffici ts. Lin [6] studied the vibration characteristic for a sandwich beam with a silica/polymer blend as the principal material, and a pure polymer matrix as the surface laminate. A numerical model has been developed and validated experimentally for the vibration of a symmetric elastic-viscoelastic sandwich beam. Moita et al. [7] developed a finite element model for vibration analysis of active-passive damped multilayer sandwich plates, with a viscoelastic core sandwiched between elastic layers, including piezoelectric layers. The elastic layers are modeled using the classical plate theory and the core is modeled using the Reissener-Mindlin theory. Meftah et al. [8] presented an alternative technological solution for passive earthquake control of shear building structures. Alvelid [9] developed a sixthorder differential equation for the dynamic analysis of the defl ction of a three-layer sandwich beam with a viscoelastic middle layer. Transverse shear deformations, as well as rotational inertia eff cts of the covering layers are taken into account. Iurlaro et al. [10] presented the derivation of the nonlinear equations of motion and consistent boundary conditions of the refined zigzag theory (RZT) for multilayered plates.
As previously stated, the study of forced viscoelastic sandwich structures leads to resolving complex equations. This issue has remained relevant to researchers until now, because of its diffi ulty. The problem is that almost all resolution methods need a triangulation of stiff ess matrix, which is of a complex type. This could lead to a high computational cost in the case of large-scale structures. Only a few studies have focused on cost reduction. Park et al. [11] used a condensation method to remove only the internal variables of viscoelastic properties. Chen et al. [12] proposed an iterative reduction method for viscoelastic structures with a constant complex modulus. Some authors have used a one-mode Galerkin's procedure to analyze linear and nonlinear vibrations in viscoelastic sandwich structures [13][14][15]. Lima et al. [16] suggested a robust condensation procedure combined with a sub-structuring technique, intended to be used for dealing with large-scale viscoelastically damped structures, but only in the case of forced frequency response. Chazot et al. [17] used a method based on a Padé approximant to reduce the computation time in the dynamic response computation model of multilayered viscoelastic structures. This acceleration technique leads to fast frequency sweep computations, as compared to a standard direct method. More recently, Bilasse and Oguamanam [18] used the usual modal analysis with real or complex eigenvectors to reduce the forced harmonic equation of large-scale sandwich plates with a viscoelastic core. However, the complex basis requires a nonlinear complex eigenvalue problem resolution.
In this article, a reduction method is combined with an asymptotic numerical method [5] to compute the viscoelastic sandwich plate response. The excitation is assumed to be harmonic. For numerical computation, structures are discretized using the fin te elements method (FEM). This article is presented as follows: in the first section, the governing equation for sandwich structures subjected to harmonic load is specifie . The second section presents the use of the asymptotic numerical method to solve the forced harmonic problem when the complex Young's modulus of the viscoelastic core is constant or frequency dependent. A reduction method is proposed under various basis choices. Finally, the validity and effectiveness of the present method are illustrated in several numerical examples of viscoelastic sandwich beams.

Problem formulation
The problem of interest here is the forced vibrations of elasticviscoelastic-elastic sandwich beam (Figure 1a). The constitutive materials of the beam are homogeneous and isotropic. The Young modulus of the core is complex frequency dependent, but the Poisson ratio is assumed to be constant. Because of the contrast between the moduli of the faces and the core, the viscoelastic material is subject to a large shear strain. This cannot be accounted for within the classical kinematic models. Hence, the use of zig-zag models is the most appropriate solution (Figure 1b). An essential review of applied theories and models for sandwich structures can be found in [19]. Assuming a common transverse displacement and the continuity of the displacement at the interfaces [20,21], one can built a sandwich fin te element with only eight degrees of freedom (d.o.f) per node (the longitudinal displacements of the faces, the common defl ction, and three rotations) as proposed in [22]. In this article, the sandwich fin te element developed in [23] is used. This element is based on shell elements from Büchter et al. [24].
Using the fin te element method, the governing equation for any viscoelastic structure excited by a transverse harmonic force can be written in the following form: where ω is the vibration circular frequency; U is the complex displacement vector of dimension [ND]; K is the complex stiff ness matrix of dimension [ND × ND] and depends nonlinearly on the excitation frequency ω; M is the global mass matrix of dimension [ND × ND]; F is the vector of the amplitude excitation of dimension [ND]. In many studies, the materials constituting the viscoelastic sandwich structures are assumed to be homogeneous and isotropic [5,18]. This assumption enables the breakdown of the stiffness matrix K as follows: where K 0 and K v are real constant stiffness matrices of dimension [ND × ND]. E(ω) is the complex Young's modulus of the viscoelastic core. Two models are studied in this article: a constant Young's modulus and a Young's modulus varying with respect to the circular frequency following the generalized Maxwell model [23]. Then the equation (Eq. 1) can be written in the following form: To simplify the development, Eq. (3) can be written in a single form for the two cases as follows: where As stated in the introduction, many methods exist to solve Eq. (4). The vibration problem given in Eq. (4) can be solved by the available commercial fin te element code, such as the Abaqus code. However, the determination of the response curves can require a high computational cost because of the assembly and the inversion of the dynamic stiffness matrix at each frequency. Abdoun et al. [5] presented an asymptotic numerical method for forced harmonic vibration analyses of viscoelastic structures. A mathematical formulation that may account for various viscoelastic models is presented. The solutions were obtained over a wide frequency band with a few inversions of the dynamic stiff ness matrix. More recently, Bilasse and Oguamanam [18] combined the asymptotic numerical method (ANM) and the classical modal analysis in order to reduce the size of the problem (Eq. 1). Two projection bases have been tested: the real and complex eigenvectors. However, the complex basis requires a nonlinear complex eigenvalue problem resolution.
In this article, the combination of the use of the ANM and an appropriate reduction method to obtain the forced harmonic response of viscoelastic sandwich structures is proposed.

Asymptotic numerical method
Within asymptotic numerical techniques the residual equation can be solved (Eq. 4) where the unknown X is sought in the form of a truncated power series with respect to a path parameter 'a': where X 0 is a regular solution of Eq. (4) and X p is the unknown at each order p (p = 1, n) to be computed. n is the truncation order. Series (Eq. 5) are inserted in the equations (Eq. 4). Then, from the identific tion of the similar powers of 'a' , a set of recurrent linear problems is obtained at each order p as follows: where K t0 denotes the Jacobian matrix at the initial point X 0 and F nl p is the second member vector at order p. The second member is known since it depends only on known parameters X i computed in previous orders (i = 1, p -1).
The governing system has one additional unknown relative to the number of equations. Then, another equation must be added. To do so, the path parameter 'a' is identified as the projection of the unknown vector increment (X − X 0 ), on the tangent vector X 1 as follows: where Z is a diagonal matrix, its components are equal to 1 or 0 to defin the ANM parameterization [25]. In this study, a pseudo arc-length scheme is used. The implementation of series (Eq. 5) in the latter (Eq. 7) makes it possible to provide the additional equation at each order [26].
In the next two sections, Eqs. (5) and (7) are written for the two models of the Young's modulus studied without giving development details (for more details in each case see Appendices A and B).

Constant Young's modulus
In this case, the Young's modulus of viscoelastic material is complex and constant. Thus, the only unknowns in Eq. (4) are the displacements and the square of the vibration circular frequency (parameter C in the following). They are developed into the power series as follows: where U 0 and C 0 are initial known solutions; n is the troncature order; U p and C p are the new unknowns at each order p (p = 1, n) to be computed; and ( * ) symbolizes the product operator.
The path parameter used in the series (Eq. 8) can be identifi d as the projection of the displacement increment (U − U 0 ) and the square frequency increment (C − C 0 ) on the tangent vector where <.,.> designates the Euclidian scalar product.

Variable Young's modulus
In the second model, the Young's modulus is assumed to be variable with respect to the frequency following a generalized Maxwell model: where k j and η j are material Maxwell constants (given experimentally) and 'i' is the complex number (i 2 = −1).
Then, in addition to Eq. (8), the Young's modulus and the circular frequency should be developed into a power series with respect to a path parameter 'a': where ω 0 and E 0 are initial known solutions and ω p and E p are the new unknowns at each order p (p = 1, n) to be computed. The path parameter used in the series (Eq. 8) and (Eq. 11) can be identifi d as the projection of the displacement increment (U − U 0 ) and the frequency increment (ω − ω 0 ) on the tangent vector (U 1 , ω 1 ): where <.,.> designates the Euclidian scalar product.

Resolution and validity range
The resolution of Eq. (3), using the ANM, is widely explained in different papers [5,22]. The polynomial solutions coincide almost perfectly within the convergence radius, but they diverge out of this zone of validity. This limit can be computed automatically by using the following simple criterion proposed by Cochelin et al. [27]: where ε is a small fi ed number. At this stage, the equation to be solved has a high number of degrees of freedom. In order to reduce it, as well as the computation time, the reduced-order models are introduced in the discrete form (Eq. 1). This procedure is described in the following section.

Reduction model
As previously stated, the resolution of the governing equation (3) entails a considerable computational cost, especially in the case of large-scale structures. This is mainly due to the complex tangent matrix triangulations required. In order to reduce the computational cost, a reduction technique has been developed and applied to Eq. (3). This technique involves projecting the displacement vector onto a small basis as follows: where is the projection matrix of dimension [ND × nd] and u is the reduced vector of dimension [nd]. The reduction model must generate small and limited approximation errors. The system properties, such as stability, must also be preserved. The key point of the reduction technique is the choice of the reduction matrix . This one must adequately characterize the dynamic response of the structure, and it must be able to approximate the solution in a significant time interval. Moreover, its columns must be linearly independent. In this study, three reduction bases are tested.
To obtain the reduced governing equation, the displacement vector U is replaced by its expression (Eq. 14), in Eq. (3). After that, both members of Eq. (3) are multiplied at the left by T . Then, the governing equation to be resolved becomes smaller. Its dimension is now [nd] with (nd ND) and the reduced problem can be written as follows: where k 0 = T K 0 , k v = T K v , m = T M , and f = T F. Then, the reduced problem (Eq. 15) is similar to the initial one (Eq. 3). Its resolution can be achieved by any method. In this study, it is performed by the ANM but only reduced matrices and vectors are considered (Figure 2). This size reduction of the initial problem may considerably decrease the computational cost. This is demonstrated through the numerical examples examined in the last sections. Before that, the three reduction bases are detailed.

Eigenvector basis (EB)
The first matrix contains the linear eigenvectors of the equivalent elastic problem: Then the reduction matrix contains the 'n' first eigenvectors: This basis is denoted in the following EB.

Enriched eigenvector basis (EEB)
The second basis contains the previous linear eigenvectors enriched by vectors containing viscoelastic properties [28]: For this reason, the basis is named the enriched eigenvector basis (EEB).

Asymptotic numerical basis (ANB)
In this last projection matrix, one ANM step without reduction is performed and displacement vectors at each order U i issued from series (Eqs. 8) are saved to build the basis. This kind of technique has been successfully applied to solve nonlinear vibrations [29], to defin an effici t prediction-correction scheme [30], or to compute bifurcations in fluid mechanics [31]: This basis is named the asymptotic numerical basis (ANB).

Numerical examples
Two examples of cantilever beams are tested to validate the proposed reduced order models (ROM). For the first one, the complex Young's modulus is considered constant and for the second one, it depends on frequency. Before verifying the proposed reduced models, validation of the full order model (FOM) is performed in the first example by comparing it to results obtained with the commercial Abaqus code. The finite element developed in [23] is used for the discretization. Using the procedure proposed in [22], this element is obtained by assembling three fin te elements throughout the thickness of the sandwich structure, two fin te elements from Büchter et al. [24] for the elastic layers, and a volume element for the viscoelastic core. Considering the classical assumptions used for modeling sandwich structures [20,21], a shell element is obtained with eight degrees of freedom (d. r For the reduced order model (Eq. 15): where r is the residual value in the reduced order:

Example 1
In this example a cantilever beam is considered, with the dimensions (177.8 * 12.7 mm²), that has been analyzed theoretically and experimentally in the literature (Figure 1). The material characteristics are given in Table 1. The Young's modulus of the core is a constant complex number E = E 0 (1 + iη c ) where (η c = 1.5) represents the loss factor. The beam is discretized into 569 nodes (ND = 4552 d.o.f). The transverse excitation force (10 N) is assumed to be applied to the free side at coordinates (177.8 mm, 6.35 mm). First, to check the validity of ANM, results are compared with those obtained by: • Commercial software (Abaqus) using 528 elements of type C3D20R (a 20-node quadratic brick, reduced integration); • Direct computation by inversing the tangent matrix at each frequency as follows: For ANM, the following parameters are used: the truncated order n = 20 and tolerance ε = 10 −7 . The results from the direct method, the ANM, and Abaqus are superposed in Figure 3. The first curves, superposed in Figure 3a, represent the real part of transverse displacement w r evolution with respect to the frequency at the load position. The second curves, in Figure 3b, represent the global residual value evolution with respect to frequency. These figu es show that the results from the ANM program coincide perfectly with direct computation. In addition, the residual value remains very small in all of the frequency ranges.
Subsequently, results from the full order model (using ANM without reduction) and those from the reduced order model are superposed for comparison ( Figure 4). Results represent the evolution of the transverse displacement modulus at the excitation point with respect to the frequency: where w r and w i are the real and imaginary parts of the transverse displacement. Note that the two last bases give good results in comparison to those from the full-order model. Reduced and full-order residual value evolution is presented as well to view the global result quality (  remark that the global residual value, in full-order model and using the fi st basis, approaches zero in all of the frequency range. In Figure 6, more than one mode is shown to estimate the validity range of the proposed reduction model. The transverse excitation force (of magnitude 100 N) is assumed to be applied at position (x = 57.79 mm, y = 6.35 mm) to show more than one mode. Displacement is measured in the excitation position ( Figure 6). Figures clearly show good quality results.

Example 2
In this example, a cantilever beam is also considered. The dimensions of the previous beam example and excitation force are retained. Only the material properties have changed ( Table 2). In this case, the Young's modulus depends on a frequency using the generalized Maxwell model. The beam is discretized into 569 nodes (ND = 4552 d.o.f). The temperature used for simulation is equal to 20°C. As it was carried out in the  Figure 7). Results from the full-order model and those from the reduced-order model are superposed for comparison (Figure 8). The results of Figure 8 represent the evolution of the transverse displacement at the excitation point with respect to frequency. Note that the basis (EB, based on a linear eigenvalue problem) is not effici t. The two other bases (EEB and ANB) give good results, which coincide perfectly with those obtained from a full calculation around the first mode. In Figure 9, the residual values evolutions using the three reduced matrices are compared. Figure 9a shows the residual values using reduction methods, which are computed in the reduced order equation form (Eq. 21). Figure 9b shows the residual values using reduction methods computed in the full-order size equation form (Eq. 20). To obtain Figure 9b, the displacement vector in the complete-order form is computed from the reduced displacement vector using Eq. (14). Hereafter, Eq. (20) is used to compute the residual value in a full-order model. Note that the residual values (logarithm value), computed in the reduced-order equation (21), are really smaller than (-7) (Figure 9a). However, the residual values computed in the full-order size equation (20) approach zero but stay smaller than zero using the two last bases (Figure 9b). The residual value in full-order form is poor in terms of error accumulation since the resolution is entirely carried out using the reduced order.

Computational time
In this section, the computational time is analyzed using the second example (4.    Table 3. From this table, note that one step with ANM takes 6 s but with the reduction method it takes 0.20 s. However, to obtain solutions for a large band of frequency more than one step is carried out. For the last example (4.2), to fin solutions from 0 to 150 Hz, 20 iterations with ANM are required, 29 iterations using the second basis (ANM + EEB), and 201 using the third basis (ANM + ANB). To estimate the computational time, the time taken to build the basis and to compute the reduced initial vectors and matrices was considered. The gain in computational time is better with ANB and equal to 83%.

Conclusion
In this article, a reduction method has been developed for computing the response curves of viscoelastic sandwich structures submitted to a harmonic force. Three bases are tested to find the displacement response of beams. The first basis, which consists of linear eigenvectors, does not reproduce the same results as a full order model. The second basis, where the viscoelastic properties are taken into account in its construction, gives a good result in the case of the beam. The last basis is built from the first step of ANM in the full order. However, this last basis contains all the problem characteristics, even the excitation force. The results of the reduction method using this last basis are obtained with good precision and signific nt computational time gain. As a continuity of this work, it would be interesting to apply the reduction method using the last basis to analyze plate and shell structures for a short time duration. In addition to material nonlinearity, since thin structures are studied here, geometrical nonlinearity should be taken into account [26].
where U 0 and C 0 are initial known solutions; n is the troncature order; and U p and C p are the new unknowns at each order p (p = 1, n) to be computed. The path parameter used in the series (A.2) can be identifi d as the projection of the displacement increment (U-U 0 ), and the frequency increment (C -C 0 ), on the tangent vector (U 1 , C 1 ): where <.,.> designates the Euclidian scalar product. It is mentioned above that the Young's modulus E and the displacement vector U are complexes and can be broken down into real and imaginary parts, as follows: By introducing expressions (A.2) into Eqs. (A.1) and (A.3) and equating similar powers of 'a' , the following set of linear problems are obtained: order 0

Appendix B: Young's modulus variable
The governing equation, for a structure excited by a force, F, is given with respect to circular frequency and displacement as follows: Young's modulus variation is modeled by the generalized Maxwell model: The Young's modulus should be developed into a truncated power series with respect to a path parameter 'a': where E 0 is the initial known solution and E p is the new unknown at each order p (p = 1, n) to be computed. The Young's modulus is a complex number. First, develop the expression under the sum: By replacing series (B.3) into (B.1) and identifying a similar power of 'a' , the real E R pj and imaginary E I pj parts of Young's modulus are found at each order p, where: The representations of the last values are listed in Table B.1.

The problem to be resolved at each order
Here, in addition to Eqs. (A.2) and (B.3), the circular frequency should be developed into a power series with respect to a path parameter 'a': where ω 0 is the initial known solution and ω p is the new unknown at each order p (p = 1, n) to be computed.
where B is the same as given before (A.6). Preparing a second member for the following order: Preparing a second member for the following order: