# Forced harmonic response of viscoelastic sandwich beams by a reduction method

Faiza Boumediene, El Mostafa Daya, Jean-Marc Cadou, Laëtitia Duigou

► To cite this version:

Faiza Boumediene, El Mostafa Daya, Jean-Marc Cadou, Laëtitia Duigou. Forced harmonic response of viscoelastic sandwich beams by a reduction method. Mechanics of Advanced Materials and Structures, 2016, 23 (11), pp.1290-1299. 10.1080/15376494.2015.1068408. hal-01514672

HAL Id: hal-01514672

https://hal.univ-lorraine.fr/hal-01514672

Submitted on 26 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
# Forced harmonic response of viscoelastic sandwich beams by a reduction method

Faiza Boumediene<sup>a</sup>, El Mostafa Daya<sup>b,c</sup>, Jean-Marc Cadou<sup>d</sup>, and Laetitia Duigou<sup>d</sup>

<sup>a</sup>Laboratoire de Mécanique Avancée, Faculte de Genie Mecanique & Genie des Procedes, Université des Sciences et de la Technologie Houari Boumediene, Alger, Algeria; <sup>b</sup>Laboratoire d'Etude des Microstructures et Mécanique des Matériaux, Université de Lorraine, Metz Cedex, France; <sup>c</sup>Laboratory of Excellence on Design of Alloy Metals for Low-Mass Structures (DAMAS), Universite de Lorraine, Metz Cedex, France; <sup>d</sup>Laboratoire d'Ingénierie des Matériaux de Bretagne, Université de Bretagne Sud, Lorient Cedex, France

## 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.

## KEYWORDS

sandwich; viscoelastic; reduction model; beam; harmonic; asymptotic numerical method

## 1. 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 coefficients. 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 sixth-order differential equation for the dynamic analysis of the deflection of a three-layer sandwich beam with a viscoelastic middle layer. Transverse shear deformations, as well as rotational inertia effects 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 difficulty. The problem is that almost all resolution methods need a triangulation of stiff element 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–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 finite elements method (FEM). This article is presented as follows: in the first section, the governing equation for sandwich structures subjected to harmonic load is specified. 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.

## 2. Problem formulation

The problem of interest here is the forced vibrations of elastic-viscoelastic-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 build a sandwich finite element with only eight degrees of freedom (d.o.f) per node (the longitudinal displacements of the faces, the common deflection, and three rotations) as proposed in [22]. In this article, the sandwich finite element developed in [23] is used. This element is based on shell elements from Büchter et al. [24].

<figure><img src="image_2.png" /><figcaption>Figure 1. The studied beam.</figcaption></figure>

Using the finite element method, the governing equation for any viscoelastic structure excited by a transverse harmonic force can be written in the following form:

$$[\mathbf{K}(\omega) - \omega^2 \mathbf{M}] \mathbf{U} = \mathbf{F}, \quad (1)$$

where $\omega$ is the vibration circular frequency; $\mathbf{U}$ is the complex displacement vector of dimension $[ND]$; $\mathbf{K}$ is the complex stiffness matrix of dimension $[ND \times ND]$ and depends nonlinearly on the excitation frequency $\omega$; $\mathbf{M}$ is the global mass matrix of dimension $[ND \times ND]$; $\mathbf{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 $\mathbf{K}$ as follows:

$$\mathbf{K}(\omega) = \mathbf{K}_0 + E(\omega) \mathbf{K}_v, \quad (2)$$

where $\mathbf{K}_0$ and $\mathbf{K}_v$ are real constant stiffness matrices of dimension $[ND \times ND]$. $E(\omega)$ 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:

$$[\mathbf{K}_0 + E(\omega) \mathbf{K}_v - \omega^2 \mathbf{M}] \mathbf{U} = \mathbf{F}. \quad (3)$$

To simplify the development, Eq. (3) can be written in a single form for the two cases as follows:

$$\mathbf{R}(\mathbf{X}) = 0, \quad (4)$$

where

$$\mathbf{R}(\mathbf{X}) = \mathbf{F} - [\mathbf{K}_0 + E(\omega) \mathbf{K}_v - \omega^2 \mathbf{M}] \mathbf{U} \quad \text{and} \\ \mathbf{X} = (\omega^2, \mathbf{U}) \quad \text{for } E \text{ constant,}$$

$$\mathbf{R}(\mathbf{X}) = \mathbf{F} - [\mathbf{K}_0 + E(\omega) \mathbf{K}_v - \omega^2 \mathbf{M}] \mathbf{U} \quad \text{and} \\ \mathbf{X} = (\omega, \omega^2, \mathbf{U}, E) \quad \text{for } E \text{ variable.}$$

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 finite 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

<figure><img src="image_3.png" /></figure>
a wide frequency band with a few inversions of the dynamic stiff ness matrix. More recently, Bilasse and Oguamanam [18] com- bined the asymptotic numerical method (ANM) and the clas- sical modal analysis in order to reduce the size of the problem (Eq. 1). Two projection bases have been tested: the real and com- plex eigenvectors. However, the complex basis requires a nonlin- ear 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.

### 3. 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':

$$X(a) = X_0 + \sum_{p=1}^{n} a^p X_p, \quad (5)$$

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 identification of the similar powers of 'a', a set of recurrent linear problems is obtained at each order $p$ as follows:

$$K_{40}X_p = F_p^{nl}, \quad (6)$$

where $K_{40}$ denotes the Jacobian matrix at the initial point $X_0$ and $F_p^{nl}$ 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:

$$a = X_1 Z (X - X_0), \quad (7)$$

where Z is a diagonal matrix, its components are equal to 1 or 0 to define 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).

### 3.1. 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:

$$U = U_0 + \sum_{p=1}^{n} a^p * U_p,$$

$$\omega^2 = C = C_0 + \sum_{p=1}^{n} a^p * C_p, \quad (8)$$

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 identified as the projection of the displacement increment ($U - U_0$) and the square frequency increment ($C - C_0$) on the tangent vector ($U_1, C_1$):

$$a = \langle U - U_0, U_1 \rangle + (C - C_0) C_1, \quad (9)$$

where <.,> designates the Euclidian scalar product.

### 3.2. 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:

$$E(\omega) = k_0 + \eta_0 \cdot i\omega + \sum_{j=1}^{N_{\max}} i\omega / \left( \frac{i\omega}{k_j} + \frac{1}{\eta_j} \right), \quad (10)$$

where $k_j$ and $\eta_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':

$$\begin{aligned} E &= E_0 + \sum_{p=1}^{n} a^p * E_p, \\ \omega &= \omega_0 + \sum_{p=1}^{n} a^p * \omega_p, \end{aligned} \quad (11)$$

where $\omega_0$ and $E_0$ are initial known solutions and $\omega_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 identified as the projection of the displacement increment ($U - U_0$) and the frequency increment ($\omega - \omega_0$) on the tangent vector ($U_1, \omega_1$):

$$a = \langle U - U_0, U_1 \rangle + (\omega - \omega_0) \omega_1, \quad (12)$$

where <.,> designates the Euclidian scalar product.

### 3.3. 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]:

$$a_{\max} = \left( \varepsilon \frac{\|U_1\|}{\|U_n\|} \right)^{1/(n-1)}, \quad (13)$$
where $\epsilon$ is a small fixed 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.

### 4. 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:

$$U = \mathfrak{R}u, \qquad (14)$$

where $\mathfrak{R}$ is the projection matrix of dimension $[ND \times 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 $\mathfrak{R}$. 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 $\mathfrak{R}^T$. Then, the governing equation to be resolved becomes smaller. Its dimension is now $[nd]$ with ($nd \ll ND$) and the reduced problem can be written as follows:

$$[\mathbf{k}_0 + E(\omega) \mathbf{k}_v - \omega^2 \mathbf{m}] \mathbf{u} = \mathbf{f}, \qquad (15)$$

where $\mathbf{k}_0 = \mathfrak{R}^T \mathbf{k}_0 \mathfrak{R}$, $\mathbf{k}_v = \mathfrak{R}^T \mathbf{k}_v \mathfrak{R}$, $\mathbf{m} = \mathfrak{R}^T \mathbf{M} \mathfrak{R}$, and $\mathbf{f} = \mathfrak{R}^T \mathbf{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.

### 4.1. Eigenvector basis (EB)

The first matrix contains the linear eigenvectors of the equivalent elastic problem:

$$[\mathbf{K}_0 - \omega_i^2 \mathbf{M}] \mathbf{\Phi}_i = \mathbf{0}. \qquad (16)$$

Then the reduction matrix contains the 'n' first eigenvectors:

$$\mathfrak{R} = \mathbf{\Phi} = [\mathbf{\Phi}_1 \dots \mathbf{\Phi}_n]. \qquad (17)$$

This basis is denoted in the following EB.

<figure><img src="image_4.png" /><figcaption>Figure 2. Resolution chart.</figcaption></figure>

### 4.2. Enriched eigenvector basis (EEB)

The second basis contains the previous linear eigenvectors enriched by vectors containing viscoelastic properties [28]:

$$\mathfrak{R} = [\mathbf{\Phi} \mathbf{K}_0^{-1} \mathbf{K}_v \mathbf{\Phi}]. \qquad (18)$$

For this reason, the basis is named the enriched eigenvector basis (EEB).

### 4.3. 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 define an efficient prediction-correction scheme [30], or to compute bifurcations in fluid mechanics [31]:

$$\mathfrak{R} = [U_1 \dots U_n]. \qquad (19)$$

This basis is named the asymptotic numerical basis (ANB).
## 5. 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 finite elements throughout the thickness of the sandwich structure, two finite 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.o.f) by node. For all examples, the following parameters are used: ANM truncation order $n = 20$ and 20 vectors in the projection basis ($nd = 20$). The convergence test is performed by computing the residual value at each step as follows:

* For the full order model (Eq. 3):

$$ RES = LOG_{10} (\|\mathbf{R}(\mathbf{X})\| / \|\mathbf{K}_0 \cdot \mathbf{U}\|). \quad (20) $$

* For the reduced order model (Eq. 15):

$$ res = LOG_{10} (\|\mathbf{r}(\mathbf{x})\| / \|\mathbf{k}_0 \cdot \mathbf{u}\|), \quad (21) $$

where $\mathbf{r}$ is the residual value in the reduced order: $\mathbf{r} = \mathbf{f} - [\mathbf{k}_0 + E(\omega)\mathbf{k}_v - \omega^2\mathbf{k}]\mathbf{u}$.

### 5.1. 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\eta_c)$ where ($\eta_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);

<table><caption>Table 1. Material properties, beam 1.</caption><thead><tr><th>Material properties</th><th>Elastic layers</th><th>Viscoelastic layer</th></tr></thead><tbody><tr><td>Young's modulus</td><td>E<sub>a</sub> = 6.9 10<sup>10</sup> N/m<sup>2</sup></td><td>E<sub>0</sub> = 1794 10<sup>3</sup> N/m<sup>2</sup></td></tr><tr><td>Poisson's ratio</td><td>ν<sub>a</sub> = 0.3</td><td>ν<sub>p</sub> = 0.3</td></tr><tr><td>Density</td><td>ρ<sub>a</sub> = 2766 kg/m<sup>3</sup></td><td>ρ<sub>p</sub> = 968.1 kg/m<sup>3</sup></td></tr><tr><td>Thickness</td><td>h<sub>a</sub> = 1.524 mm</td><td>h<sub>p</sub> = 0.127 mm</td></tr></tbody></table>

* Direct computation by inverting the tangent matrix at each frequency as follows:

$$ \mathbf{U}(\omega) = [\mathbf{K}_0 + E(\omega)\mathbf{K}_v - \omega^2\mathbf{M}]^{-1} \cdot \mathbf{F}. $$

For ANM, the following parameters are used: the truncated order $n = 20$ and tolerance $\epsilon = 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 figures 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:

$$ \|\mathbf{w}\| = \sqrt{w_r^2 + w_i^2}, \quad (22) $$

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 (Figure 5). The first figure (Figure 5a) shows the residual values, using reduction methods, which are computed in the reduced order equation form (Eq. 21). The second one (Figure 5b) shows the residual values using the reduction methods computed in the full-order size equation form (Eq. 20). To obtain this last figure (Figure 5b), 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 the full-order model. One can

<figure><img src="image_5.png" /><figcaption>Figure 3. Full order model results—constant Young's modulus.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. Transverse displacement modulus (comparison between full order model—reduced order model).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Residual value evolution vs. circular frequency.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Transverse displacement modulus for beam 1, F = 100 N at position (x = 57.79 mm, y = 6.35 mm).</figcaption></figure>

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.

5.2. 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 fre-
quency using the generalized Maxwell model. The beam is dis-
cretized 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
<table><caption>**Table 2.** Material properties, beam 2.</caption><thead><tr><th>Material properties</th><th>Elastic layers</th><th>Viscoelastic layer</th></tr></thead><tbody><tr><td>Young's modulus (N/m²)</td><td>2.1 × 10<sup>11</sup></td><td>27.216 × 10<sup>6</sup></td></tr><tr><td>Poisson's ratio</td><td>0.3</td><td>0.44</td></tr><tr><td>Density (kg/m³)</td><td>7800</td><td>1200</td></tr><tr><td>Thickness (mm)</td><td>1.524</td><td>0.127</td></tr></tbody></table>

previous example, the ANM results are validated by the direct method (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 efficient. 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.

### 5.3. Computational time

In this section, the computational time is analyzed using the second example (4.2). The different time costs are presented in

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Full order model results—variable Young's modulus.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution of transverse displacement modulus with respect to frequency (comparison: full order model—reduced order model).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Residual value in the reduced form [Eq. (21)] (b) Residual value in the full order form [Eq. (20)]</figcaption></figure>

<caption>**Figure 9.** Residual value variation (comparison: full order model—reduced order model).</caption>

<table><caption>**Table 3.** Computational CPU time (ND = 4552 d.o.f., nd = 20).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">One step time</th><th rowspan="2">Basis building time</th><th rowspan="2">Reduction of initial vectors and matrices</th><th colspan="2">0 Hz &lt; f &lt; 150 Hz</th><th rowspan="2">Gain (%)</th></tr><tr><th>Steps number</th><th>Computational time</th></tr></thead><tbody><tr><td>ANM</td><td>6 s</td><td></td><td></td><td>20</td><td>120 s</td><td></td></tr><tr><td>ANM + EEB (basis 2)</td><td>0.20 s</td><td>53 s</td><td>10 s</td><td>29</td><td>69 s</td><td><strong>42.50</strong></td></tr><tr><td>ANM + ANB (basis 3)</td><td>0.20 s</td><td>6 s</td><td>10 s</td><td>21</td><td>20 s</td><td><strong>83.33</strong></td></tr></tbody></table>

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 find 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%.

## 6. 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 significant 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].

## References

[1] E.M. Daya and M. Potier-Ferry, A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures, Comput. Struct., vol. 79, pp. 533–541, 2001.

[2] Y.P. Lu, J.W. Killian, and G.C. Everstine, Vibrations of three layered damped sandwich plate composites, J. Sound Vib., vol. 64, no. 1, pp. 63–71, 1979.

[3] T.T. Babert, R.A. Maddox, and C.E. Orozco, A fin te element model for harmonically excited viscoelastic sandwich beams, Comput. Struct., vol. 65, no. 1, pp. 105–113, 1998.

[4] L. Hazard and P. Bouillard, Structural dynamics of viscoelastic sandwich plates by the partition of unity finite element method, Comput. Methods Appl. Mech. Eng., vol. 196, pp. 4101–4116, 2007.

[5] F. Abdoun, L. Azrar, E.M. Daya, and M. Potier-Ferry, Forced harmonic response of viscoelastic structures by an asymptotic numerical method, Comput. Struct., vol. 87, pp. 91–100, 2009.

[6] J.C. Lin, Characterization and numerical evaluation of vibration on elastic-viscoelastic sandwich structures, Compos. Struct., vol. 92, pp. 669–675, 2010.

[7] J.S. Moita, A.L. Araújo, P. Martins, C.M. Mota Soares, and C.A. Mota Soares, A fin te element model for the analysis of viscoelastic sandwich structures, Comput. Struct., vol. 89, pp. 1874–1881, 2011.

[8] S.A. Meftah, E.M. Daya, and A. Tounsi, Finite element modelling of sandwich box column with viscoelastic layer for passive vibrations control under seismic loading, Thin-Walled Struct., vol. 51, pp. 174–185, 2012.

[9] M. Alvelid, Sixth order differential equation for sandwich beam deflection including transverse shear, Compos. Struct., vol. 102, pp. 29–37, 2013.

[10] L. Iurlaro, M. Gherlone, M. Di Sciuva, and A. Tessler, Assessment of the refined zigzag theory for bending, vibration, and buckling of sandwich plates: A comparative study of different theories, Compos. Struct., vol. 106, pp. 777–792, 2013.

[11] C.H. Park, D.J. Inman, and M.J. Lam, Model reduction of viscoelastic finite element models, J. Sound Vib., vol. 219, pp. 619–637, 1999.

[12] X. Chen, H.L. Chen, and X.L. Hu, Damping prediction of sandwich structures by order-reduction-iteration approach, J. Sound Vib., vol. 222, pp. 803–812, 1999.

[13] E.M. Daya, L. Azrar, and M. Potier-Ferry, An amplitude equation for the nonlinear vibration of viscoelastically damped sandwich beams, J. Sound Vib., vol. 271, pp. 789–813, 2004.

[14] M. Bilasse, E.M. Daya, and L. Azrar, Linear and nonlinear vibrations analysis of viscoelastic sandwich beams, J. Sound Vib., vol. 329, pp. 4950–4969, 2010.

[15] M. Bilasse, L. Azrar, and E.M. Daya, Complex modes based numerical analysis of viscoelastic sandwich plates vibrations, Comput. Struct., vol. 89, pp. 539–555, 2011.

[16] A.M.G. de Lima, A.R. da Silva, D.A. Rade, and N. Bouhaddi, Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures, Eng. Struct., vol. 32, pp. 1479–1488, 2010.

[17] J.D. Chazot, B. Nennig, and A. Chettah, Harmonic response computation of viscoelastic multilayered structures using a ZPST shell element, Comput. Struct., vol. 89, pp. 2522–2530, 2011.

[18] M. Bilasse and D.C.D. Oguamanam, Forced harmonic response of sandwich plates with viscoelastic core using reduced-order model, Compos. Struct., vol. 105, pp. 311–318, 2013.
[19] H. Hu, S. Belouettar, M. Potier-Ferry, and E.M. Daya, Review and assessment of various theories for modeling sandwich composites, Compos. Struct., vol. 84, pp. 282–292, 2008.

[20] D.K. Rao, Frequency and loss factor of sandwich beams under various boundary conditions, J. Mech. Eng. Sci., vol. 20, pp. 271–282, 1978.

[21] B.A Ma and J.F. He, Finite element analysis of viscoelastically damped sandwich plates, J. Sound Vib., vol. 52, pp. 107–123, 1992.

[22] E.M. Daya and M. Potier-Ferry, A shell finite element for viscoelastically damped sandwich structures, Rev. Euro. Elém. Finis, vol. 11, pp. 39–56, 2002.

[23] L. Duigou, E.M, Daya, and M. Potier-Ferry, Iterative algorithms for nonlinear eigenvalue problems: Application to vibrations of viscoelastic shells. Comput. Methods Appl. Mech. Eng., vol. 192, pp. 1323–1335, 2003.

[24] N. Büchter, E. Ramm, and D. Roehl, Three dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept, Int. J. Numer. Methods Eng., vol. 37, pp. 2551–2568, 1994.

[25] H. Mottaqui, B. Braikat, and N. Damil, Discussion about parametrization in the asymptotic numerical method: Application to nonlinear elastic shell, Comput. Methods Appl. Mech. Eng., vol. 199, no. 25–28, pp. 1701–1709, 2010.

[26] F. Boumediene, A. Miloudi, J.M. Cadou, L. Duigou, and E.H. Boutyour, Nonlinear forced vibration of damped plates by an asymptotic numerical method, Comput. Struct., vol. 87, pp. 1508–1515, 2009.

[27] B. Cochelin, N. Damil, and M. Potier-Ferry, Asymptotic numerical methods and Padé approximants for non-linear elastic structures, Int. J. Numer. Methods Eng., vol. 37, pp. 1187–213, 1994.

[28] F. Boumediene, J.M. Cadou, L. Duigou, and E.-M. Daya, A reduction model for eigensolutions of damped viscoelastic sandwich structures, Mech. Res. Commun., vol. 57, pp. 74–81, 2014.

[29] F. Boumediene, L. Duigou, E.H. Boutyour, A. Miloudi, and J.M. Cadou, Nonlinear forced vibration of damped plates coupling asymptotic numerical method and reduction models, Comput. Mech., vol. 47, no. 4, pp. 359–377, 2011.

[30] J.M. Cadou, N. Damil, M. Potier-Ferry, and B. Braikat, Projection technique to improve high order iterative correctors, Finite Elem. Anal. Des., vol. 41, no. 3, pp. 285–309, 2004.

[31] J. Heyman, G. Girault, Y. Guevel, C. Allery, A. Hamdouni, and J.M. Cadou, Computation of Hopf bifurcations coupling reduced order models and the asymptotic numerical method, Comput. Fluids, vol. 76, pp. 73–85, 2013.

## Appendix A: Constant Young's modulus

The governing equation, for a structure excited by a force, F, is given with respect to circular frequency and displacement as follows:

$$[\mathbf{K}_0 + E \cdot \mathbf{K}_v - \omega^2 \mathbf{M}] \mathbf{U} = \mathbf{F}. \quad (A.1)$$

In this case, the unknown is developed into power series as follows:

$$\begin{align*}
\mathbf{U} &= \mathbf{U}_0 + \sum_{p=1}^{n} a^p * \mathbf{U}_p, \\
\omega^2 &= C_0 + \sum_{p=1}^{n} a^p * C_p, \tag{A.2}
\end{align*}$$

where $\mathbf{U}_0$ and $C_0$ are initial known solutions; n is the tronce order; and $\mathbf{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 identified as the projection of the displacement increment ($\mathbf{U}-\mathbf{U}_0$), and the frequency increment ($C-C_0$), on the tangent vector ($\mathbf{U}_1, C_1$):

$$a = \langle \mathbf{U}-\mathbf{U}_0, \mathbf{U}_1 \rangle + (C-C_0)C_1, \quad (A.3)$$

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:

$$\begin{align}
\mathbf{U} &= \mathbf{U}^R + i\mathbf{U}^I, \\
E &= E^R + iE^I. \tag{A.4}
\end{align}$$

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

$$[\mathbf{B}] \begin{Bmatrix} \mathbf{U}_0^R \\ \mathbf{U}_0^I \end{Bmatrix} = \begin{Bmatrix} \mathbf{F} \\ 0 \end{Bmatrix}, \quad (A.5)$$

where

$$[\mathbf{B}] = \begin{bmatrix} \mathbf{K} & -E_0^I \mathbf{K}_v \\ -E_0^I \mathbf{K}_v & -\mathbf{K} \end{bmatrix}; \quad [\mathbf{K}] = [\mathbf{K}_0 + E_0^R \mathbf{K}_v - C_0 \mathbf{M}]; \quad (A.6)$$

### order 1

$$[\mathbf{B}] \begin{Bmatrix} \mathbf{U}_1^R \\ \mathbf{U}_1^I \end{Bmatrix} = C_1 [M] \begin{Bmatrix} \mathbf{U}_0^R \\ \mathbf{U}_0^I \end{Bmatrix}, \\
\mathbf{U}_1^T \mathbf{U}_1 + C_1^2 = 1; \quad (A.7)$$

### order p

$$[\mathbf{B}] \begin{Bmatrix} \mathbf{U}_p^R \\ \mathbf{U}_p^I \end{Bmatrix} = \sum_{s=0}^{p-1} C_{p-s} [\mathbf{M}] \begin{Bmatrix} \mathbf{U}_s^R \\ \mathbf{U}_s^I \end{Bmatrix}, \\
\mathbf{U}_1^T \mathbf{U}_p + C_1 C_p = 0. \quad (A.8)$$

## 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:

$$[\mathbf{K}_0 + E(\omega) \mathbf{K}_v - \omega^2 \mathbf{M}] \mathbf{U} = \mathbf{F}. \quad (B.1)$$

Young's modulus variation is modeled by the generalized Maxwell model:

$$E(\omega) = k_0 + \eta_0 \cdot i\omega + \sum_{j=1}^{N_{\max}} i\omega / \left( \frac{i\omega}{k_j} + \frac{1}{\eta_j} \right). \quad (B.2)$$

The Young's modulus should be developed into a truncated power series with respect to a path parameter 'a':

$$E = E_0 + \sum_{p=1}^{n} a^p * E_p, \quad (B.3)$$

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:

$$E_j = i\omega / \left( \frac{i\omega}{k_j} + \frac{1}{\eta_j} \right), \quad (B.4)$$
<table><caption>**Table B.1.** Real $E_{pj}^R$ and imaginary $E_{pj}^I$ parts of the Young's modulus at each order.</caption><thead><tr><th>Order</th><th>E<sup>R</sup><sub>pj</sub></th><th>E<sup>I</sup><sub>pj</sub></th></tr></thead><tbody><tr><td>0</td><td>C<sub>0</sub>k<sub>j</sub> / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td><td>&omega;<sub>0</sub>(k<sub>j</sub><sup>2</sup>/&eta;<sub>j</sub>) / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td></tr><tr><td>1</td><td>C<sub>1</sub> (k<sub>j</sub> - E<sup>R</sup><sub>0j</sub>) / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td><td>(&omega;<sub>1</sub> (k<sub>j</sub><sup>2</sup>/&eta;<sub>j</sub>) - C<sub>1</sub>E<sup>I</sup><sub>0j</sub>) / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td></tr><tr><td>p</td><td>[C<sub>p</sub> (k<sub>j</sub> - E<sup>R</sup><sub>0j</sub>) - &sum;<sub>s=1</sub><sup>p-1</sup> C<sub>s</sub>E<sup>R</sup><sub>(p-s)j</sub>] / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td><td>[&omega;<sub>p</sub>(k<sub>j</sub><sup>2</sup>/&eta;<sub>j</sub>) - &sum;<sub>s=1</sub><sup>p-1</sup> C<sub>s</sub>E<sup>I</sup><sub>(p-s)j</sub>] / ((<sup>k<sub>j</sub></sup>&frasl;<sub>&eta;<sub>j</sub></sub>)<sup>2</sup> + C<sub>0</sub>)</td></tr></tbody></table>

By replacing series (B.3) into (B.1) and identifying a similar power of 'a', the real $E_{pj}^R$ and imaginary $E_{pj}^I$ parts of Young's modulus are found at each order p, where:

$$E_{pj} = E_{pj}^{R} + iE_{pj}^{I}. \qquad (B.5)$$

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':

$$\omega = \omega_0 + \sum_{p=1}^{n} a^p * \omega_p, \qquad (B.6)$$

where $\omega_0$ is the initial known solution and $\omega_p$ is the new unknown at each order $p$ ($p = 1, n$) to be computed.

In this case, the path parameter used in the series (A.2), (B.3), and (B.6) can be identified as the projection of the displacement increment $(\mathbf{U}-\mathbf{U}_0)$, and the frequency increment $(\omega - \omega_0)$, on the tangent vector $(\mathbf{U}_1, \omega_1)$.

$$a = \langle \mathbf{U} - \mathbf{U}_0, \mathbf{U}_1 \rangle + (\omega - \omega_0)\omega_1, \qquad (B.7)$$

where <.,>, designates the Euclidian scalar product.

By introducing expressions (A.2), (B.3), and (B.6) into Eqs. (B.1), (B.2), and (B.7) and equating like powers of 'a', the following set of linear problems is obtained:

#### order 0

$$ \begin{aligned} [\mathbf{B}] \begin{Bmatrix} \mathbf{U}_0^R \\ \mathbf{U}_0^I \end{Bmatrix} &= \begin{Bmatrix} \mathbf{F} \\ \mathbf{0} \end{Bmatrix}, \\ C_0 &= \omega\omega_0^2, \\ E_0^R &= \sum_{j=1}^{\max} \left[ \omega_0^2 k_j / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right) \right], \\ E_0^I &= \eta_0\omega_0 + \sum_{j=1}^{\max} \left[ \omega_0 \left( \frac{k_j^2}{\eta_j} \right) / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right) \right], \end{aligned} \qquad (B.8) $$

where **B** is the same as given before (A.6).

Preparing a second member for the following order:

$$ \begin{Bmatrix} \mathbf{F}_1^R \\ \mathbf{F}_1^I \end{Bmatrix} = \begin{Bmatrix} 2\omega_0 \cdot \mathbf{M} \cdot \mathbf{U}_0^R - E_1^{R*} \mathbf{K}_v \cdot \mathbf{U}_0^R + E_1^{I*} \mathbf{K}_v \cdot \mathbf{U}_0^I \\ -2\omega_0 \cdot \mathbf{M} \cdot \mathbf{U}_0^I + E_1^{R*} \mathbf{K}_v \cdot \mathbf{U}_0^R + E_1^{I*} \mathbf{K}_v \cdot \mathbf{U}_0^I \end{Bmatrix}, $$

$$ \begin{aligned} E_1^{R*} &= \sum_{j=1}^{\max} \left[ 2\omega_0 (k_j - E_{0j}^R) / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right) \right], \\ E_1^{I*} &= \eta_0 + \sum_{j=1}^{\max} \left[ \left( (k_j^2/\eta_j) - 2\omega_0 E_{0j}^I \right) / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right) \right]; \end{aligned} \qquad (B.9) $$

#### order 1

$$ \begin{aligned} [\mathbf{B}] \begin{Bmatrix} \mathbf{U}_1^R \\ \mathbf{U}_1^I \end{Bmatrix} &= \omega_1 \cdot \begin{Bmatrix} \mathbf{F}_1^R \\ \mathbf{F}_1^I \end{Bmatrix}, \\ \mathbf{U}_1^T \mathbf{U}_1 + \omega_1^2 &= 1, \\ C_1 &= 2\omega_0\omega_1, \\ E_1^R &= \omega_1 \cdot E_1^{R*}, \\ E_1^I &= \omega_1 \cdot E_1^{I*}. \end{aligned} \qquad (B.10) $$

Preparing a second member for the following order:

$$ \begin{aligned} E_p^{R*} &= \sum_{j=1}^{\max} \left( C_p^* (k_j - E_{0j}^R) - \sum_{s=1}^{p-1} C_s E_{(p-s)j}^R \right) / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right), \\ E_p^{I*} &= - \sum_{j=1}^{\max} \left( E_{0j}^I C_p^* + \sum_{s=1}^{p-1} C_s E_{(p-s)j}^I \right) / \left( \left( \frac{k_j}{\eta_j} \right)^2 + C_0 \right), \\ C_p^* &= \sum_{s=1}^{p-1} \omega_s \omega_{p-s}. \end{aligned} \qquad (B.11) $$

#### order p

$$ \begin{aligned} [\mathbf{B}] \begin{Bmatrix} \mathbf{U}_p^R \\ \mathbf{U}_p^I \end{Bmatrix} &= \omega_p \cdot \begin{Bmatrix} \mathbf{F}_1^R \\ \mathbf{F}_1^I \end{Bmatrix} + \begin{Bmatrix} \mathbf{F}_p^R \\ \mathbf{F}_p^I \end{Bmatrix}, \\ \mathbf{U}_1^T \mathbf{U}_n + \omega_1 \omega_n &= 0, \\ C_p &= 2\omega_0\omega_p + C_p^*, \\ E_p^R &= \omega_p \cdot E_1^{R*} + E_p^{R*}, \\ E_p^I &= \omega_p E_1^{I*} + E_p^{I*}. \end{aligned} \qquad (B.12) $$