# Vibroacoustical Analysis of Multiple-Layered Structures with Viscoelastic Damping Cores

Authors: Fei Lin; Mohan D. Rao
Journal: ISRN Mechanical Engineering
 (2013)

---

## Abstract

This paper presents a modeling technique to study the vibroacoustics of multiple-layered viscoelastic laminated beams using the Biot damping model. In this work, a complete simulation procedure for studying the structural acoustics of the system using a hybrid numerical model is presented. The boundary element method (BEM) was used to model the acoustical cavity, whereas the finite element method (FEM) was the basis for vibration analysis of the multiple-layered beam structure. Through the proposed procedure, the analysis can easily be extended to another complex geometry with arbitrary boundary conditions. The nonlinear behavior of viscoelastic damping materials was represented by the Biot damping model taking into account the effects of frequency, temperature, and different damping materials for individual layers. The curve-fitting procedure used to obtain the Biot constants for different damping materials for each temperature is explained.  The results from structural vibration analysis for selected beams agree with published closed-form results, and results for the radiated noise for a sample beam structure obtained using a commercial BEM software are compared with the acoustical results of the same beam by using the Biot damping model.


---

## Body

## 1. Introduction

The traditional designs of free-layer, constrained-layer or sandwich-layer, damping treatment using viscoelastic materials have been around for over forty years. Recent improvements in the understanding and application of the damping principles, together with advances in materials science and manufacturing, have led to many successful applications and the development of patch damping and multiple-layered damping structures. The key point in any design is to recognize that the damping material must be applied in such a way that it is significantly strained whenever the structure is deformed in the vibration mode under investigation. Numerous researchers have successfully implemented the passive constrained layer (PCL) and active constrained layer (ACL) systems. In 1959, Kerwin [1] and Ross et al. [2] presented a general analysis of viscoelastic material structure. The damping was attributed to the extension and shear deformations of the viscoelastic layers. Ditaranto [3] developed sixth-order equations of motion in terms of axial displacements and developed a closed-form solution. Mead and Markus [4] extended the sixth-order equations of motion for transverse displacement to include various boundary conditions. A paper by Rao [5] presented the equations of motion of viscoelastic sandwich beams with various boundary conditions using the energy method. The equations were solved numerically, and a practical design guideline was presented. Similar to Rao’s theory, Cottle [6] used Hamilton’s principle to derive equations of motion. The damping could also be increased by adding passive stand-off layer (PSOL) and slotted stand-off layer (SSOL) to the layered systems. Falugi [7] and Parin et al. [8] conducted theoretical and experimental work on a four-layered panel and a five-layered beam with PSOL treatment. Rogers and Parin [9] and Yellin et al. [10] have performed experimental investigations and demonstrated that PSOL treatment increased the damping significantly in aeronautical structures and beams. Yellin and colleagues [11, 12] also developed normalized equations of motion for beam, fully treated with PSOL using nonideal stand-off layer assumption. The equations were solved using the method of distributed transfer functions [13].In addition to the closed-form analytical approach, many researchers have used the finite element method (FEM), the most popular numerical modeling method in building the numerical model of the multiple layers system. In 2000, Chen and Chan [14] studied four different types of integral FEM models with the viscoelastic cores. The numerical stability and accuracy as well as for convergence issue of these four different FEM models were demonstrated by comparing the numerical results with those from experiments. Lesieutre and Lee [15] proposed a 3-node, 10 DOF FEM model for the three-layer ACL damping beam. This FEM model is advantageous in active control application due to its features of nonshear locking and adaptability to segmented constraining layers.Other than the Hamilton and FEM numerical methods in building the models, other researchers have proposed many irregular modeling techniques for the numerical representation of continuous/discontinuous systems. Kung and Singh [16] calculated the natural frequencies and loss factor using the Rayleigh-Ritz energy method and modal-strain energy technique in modeling a 3-layer patch damping structure. These approximate modeling methods were also extended to rectangular damping patch of plates and shells with viscoelastic cores. Zhang and Sainsbury [17] combined the Gerlerkin orthogonal function with the traditional finite element method and successfully applied to the vibration analysis of the damped sandwich plates.While the FEM is used widely in the modeling of the structure, many researchers sought for proper mathematical models to represent the damping behavior of the viscoelastic material, as well as incorporating the damping model in commercial FEM software packages. Currently, many FEM commercial software incorporate damping models based mostly on viscous/hysteretic damping. Some allow incorporation of damping energy dissipation in the time domain using the Prony series. None of these damping models, however, is suitable to capture the damping behavior in the frequency domain, which is the most important issue in predicting the vibro-acoustical response of complex structural systems. The drawback of these damping models raised considerable interest and motivation in the development of damping models of viscoelastic material in the frequency domain compatible with FEM software. These damping models can be classified as derivative type and integral type.The “Fractional Derivative” is essentially the representative damping model in the derivative form family proposed by Bagley and Torvik [18] in 1983. This damping model not only described the material properties of viscoelastic damping but so established the closed-form equation compatible with the FEM technique. Compared with the other integral-form models, the fractional derivative is only able to capture the relatively weak frequency-dependent information; however, it was an important milestone in the area of damping research.Lesieutre et al. [20] mathematically modeled the relaxation behavior of viscoelastic material in terms of augmenting thermodynamic field (ATF) in 1989. Initially, introducing a single augment field, this damping model provided the ability to represent the light-damping behavior, with the application of a 1D viscoelastic structure. In the subsequent research, using a series of augment fields, the ATF model is able to model the damping material of higher loss factor with the weak frequency dependence. Remedying the limitation of 1D application, Lesieutre and Lee [15] proposed an anelastic displacement field (ADF) technique in 1996 and successfully extended its application from the 1-D problem to the 3-layer sandwich beam and 3-D problems.As far as the mini-oscillator damping models are concerned, the complex shear modulus which is a function of both frequency and temperature can be expressed by a series of mini-oscillation perturbations. Biot [21] first proposed the first-order relaxation function with the introduction of the “dissipative variables” into the dynamic equations using the theory of irreversible thermodynamics. In 2007, Zhang and Zheng [22] utilized the Biot model to describe the dynamic behavior of a viscoelastic structure. The dimension reduction technique and nonlinear curve-fitting procedure were discussed in the paper. McTavish [23] developed another mini-oscillator damping model called “GHM” by the usage of Second-order relaxation function. Compared with the Biot model, the GHM model has a more complicated expression and also requires better performance of the computational tool.The popularity of these integral-form damping models in recent years brought two research interests: nonlinear curve fitting and dimension reduction. The advanced curve-fitting techniques in the damping models guarantee the accuracy of the numerical representation of the actual shear modulus data from the experiment. The dimension reduction technique increases the computational efficiency due to the additional orders of equation in order to gain the frequency independence of the frequency-form damping model.Zhang et al. [24] converted the nonlinear curve-fitting problem in frequency domain with respect to the GHM parameters into the constrained nonlinear optimization problem. The efficiency and correctness were demonstrated for a commercial viscoelastic material.Park et al. [25] examined the GHM damping model with the application to the FEM method associated with the Guyan reduction technique. The numerical example in this research leads to an FEM model applied to the GHM dynamic equation quantitatively without increasing the number of order.Hao and Rao [26] carried out the optimum design of a three-layer sandwich beam for the vibration analysis in 2005. In this research, the numerical model is a comprehensive formulation for a three-layer unsymmetrical sandwich beam with two different damping materials adjacent to each other. The criterion of the optimization is to minimize the mass of the structure while maximizing the system damping. In 2008, Lee [27] published the semicoupled vibroacoustical analysis and optimization of a simply supported three-layer sandwich beam. The modal superposition method was used to investigate the vibration problem with the fractional derivative damping model. The interior acoustical problem was studied by BEM numerical technique, and the optimization problem was established through the appropriate sizing parameters of the sandwich beam.The objective of this paper is to extend the previous work by the authors [28] on the vibration analysis of a multiple-layered beam structure incorporating the Biot damping model to solve the acoustic problem to predict the radiated noise. In this paper, we present a complete numerical procedure for the vibroacoustical analysis and design for a multiple-layer laminated damping beam. Results obtained from the proposed vibration analysis are compared with the previous closed-form results to show the validity of this approach. The radiated noise spectrum at selected field point shows good agreement between the 2-D BEM acoustical analysis and the result without system damping calculated by commercial software for a sample viscoelastic damping structure. The acoustical solution is demonstrated, and the correlation between sound pressure level (SPL) and the loss factor is also highlighted.


## 2. FEM Modeling and the Biot Dynamic Equation

The FEM modeling procedure and the establishment of the Biot dynamic equation will be discussed in this section. The structure chosen for illustration is a seven-layer viscoelastic sandwich beam. The elastic beam and the constrained damping layer are the two fundamental components in this FEM-modeling technique. The concept of transfer matrix is used to convert the local coordinates to the global coordinates in order to assemble and construct the complete model of the sandwich damping structure with arbitrary number of layers. The Biot viscoelastic damping model will be used to describe the damping behavior. Through the use of the FEM, the structure is discretized which will enable the use of the Biot damping model for different damping layers in the structure. The reader is referred to the nomenclature for the definition of different variables used in the derivation.


### 2.1. FEM Modeling of Component I: The Elastic Layer

Figure1 shows the elastic layer in the FEM model, containing 2 nodes and 6 degrees of freedom (DOF). The element displacements of each node can be expressed as follows:
(1){δ}elastice=(wiθiui∣wjθjuj)T︸6×1.Figure 1
Configuration of the elastic layer showing the DOF.The stiffness matrix can be derived based on the following energy method:(2)[Ke]elastice=∫01EAl[∂Ne∂ξ]︸6×1T[∂Ne∂ξ]︸1×6dξ+∫01EIl[∂2Nf∂2ξ]T[∂2Nf∂ξ2]dξ︸,6×6
as the shape functions are the following:
(3)[Nf]=[1-3ξ2+2ξ3(ξ-2ξ2+ξ3)l03ξ2-2ξ3(-ξ2+ξ3)l0],[Ne]=[001-ξ00ξ],
in which ξ: the local coordinate, ξ=x/l, ξ∈[0,1],  l: longitudinal length of elastic layer,  A: cross-sectional area of the elastic layer,  E: Young’s modulus of the elastic layer, and  I: moment of inertia of elastic layer.Similarly, the element mass matrix can be expressed as:(4)[M]elastice︸6×6=∫01ml([Nf]︸6×1T[Nf]︸1×6+[Ne]T[Ne])dξ.


### 2.2. FEM Modeling of Fundamental Component II: The Constrained Damping Layer

The FEM model of the constrained layout containing the damping layer sandwiched between two outer layers is shown in Figure2. This Figure illustrates each element consisting of 2 nodes and 8 DOF, where the nodal displacement vector is as follows:
(5){δ}conse=(wiθiu1iu3i∣wjθju1ju3j)T︸8×1.Figure 2
Configuration of the constrained damping layer showing the DOF.Through the introduction of transfer matrix,(6)[TI]=[e1e2e3e5e6e7]T;[TIII]=[e1e2e4e5e6e8]T,
in which each ei means the following vector: ei=(00⋯⋯1︸ithplace⋯⋯0)T; the element elastic stiffness and the element viscoelastic stiffness matrix for this 3-layer component, respectively, are the following:
(7)[Ke]conse=[TI]T·[Ke]elastic·[TI]+[TIII]T·[Ke]elastic·[TIII],[Kv]conse=∫01G2A2lkh2[Ne1-Ne3h2+h0h2·1l·∂Nf1∂ξ]T×[Ne1-Ne3h2+h0h2·1l·∂Nf1∂ξ]dξ,[Ne1]︸1×8=[Ne]︸1×6[TI]︸6×8;[Ne3]=[Ne][TIII],
where A2: cross-sectional area of the damping layer, G2: long-term shear modulus of the damping layer, and k: correction factor of the shear strain energy, for the rectangular cross-section, k=1.2.Also, the element mass matrix for this 3-layer component is(8)[M]conse=[TI]T·[M]elastic·[TI]+[M]cons,2e+[TIII]T·[M]elastic·[TIII],
where
(9)[M]cons,2e=∫01m2l[Nf1]T[Nf1]dξ,[Nf1]=[Nf][TI].


### 2.3. FEM Modeling of a Seven-Layer Constrained Damping Beam

The seven-layer sandwich beam consists of seven alternating layers—four elastic layers and three damping layers. Figure3 shows the FEM model of a seven-layer sandwich beam containing 2 nodes and 10 DOF, and the node displacement vector is as follows:
(10){δ}7layere=(wiθiu1iu3iu5iu7i∣wjθju1ju3ju5ju7j)T︸12×1.Figure 3
Configuration of a seven-layer damping structure showing the DOF.The transfer matrix to obtain the element stiffness and the mass matrix when the 1st, 3rd, 5th, and 7th layers are elastic are follows:(11)[T1]=(e1e2e3e7e8e9)T︸6×12,[T3]=(e1e2e4e7e8e10)T︸6×12,[T5]=(e1e2e5e7e8e11)T︸6×12,[T7]=(e1e2e6e7e8e12)T︸6×12.Similarly, the element stiffness and the mass matrix for the 2nd, 4th, and 6th layers of the constrained damping layer can be derived through the transfer matrix:(12)[T2]=(e1e2e3e4e7e8e9e10)T︸8×12,[T4]=(e1e2e4e5e7e8e10e11)T︸8×12,[T6]=(e1e2e5e6e7e8e11e12)T︸8×12,
where the notation ei means
(13)ei=(00⋯⋯1︸ithplace⋯⋯0)T︸12×1.Based on the above equations and design parameters of each layer, the element mass/stiffness/damping matrix of the seven-layer sandwich damping beam can be expressed as follows:(14)[Ke]e=T1T[Ke1]T1+T3T[Ke3]T3+T5T[Ke5]T5+T7T[Ke7]T7,[Kv]e=T2T[Kv2]T2+T4T[Kv4]T4+T6T[Kv6]T6,[Me]e=∑i=17TiT[Mei]Ti.Thus, the element matrices can be assembled to obtain the global mass/stiffness/damping matrix and can be applied to the boundary condition through the conventional FEM technique. Taking into the consideration of the viscoelastic damping properties, the global matrices need to be manipulated as a portion of the Biot dynamic equation.


### 2.4. Introduction of the Biot Dynamic Equation

To consider the vibration problem numerically, the dynamic equation discretized by FEM technique needs to be expressed by the following second-order ordinary differential equation (ODE) form:(15)Mx¨+Cx˙+Kx=f(t).The Biot viscoelastic damping model numerically represents the complex shear modulus with a series of mini-oscillator perturbing terms:(16)sG~(s)=G∞[1+∑k=1makss+bk],
in which G∞ is the long-term shear moduli; ak and bk are the Biot constants. These parameters are positive and can be determined by nonlinear curve fitting from the experimental data. The curve-fitting procedure will be discussed in Section 3.Substituting the Biot damping model into (15), the dynamic equation with m terms of the Biot parameters for the first viscoelastic material and n terms for the second viscoelastic material can be developed as follows:
(17)G1∞Λv1=Λ1,G2∞Λv2=Λ2,R1=G1∞Rv1Λv1,R2=G2∞Rv2Λv2,
where Rv and Λv are the eigenvector and diagonal eigenvalue matrices, respectively, from the damping matrix C. Additionally, a11⋯a1m, b11⋯b1m, and z11⋯z1m denote m terms of the Biot parameters and the dissipative coordinates, respectively, for first viscoelastic material.Similarly,a21⋯a2n, b21⋯b2n, and z21⋯z2n denote n terms of the Biot parameters and the dissipation coordinates, respectively, for second viscoelastic material. A detailed derivation can be found in the previous publication [28].


## 2.1. FEM Modeling of Component I: The Elastic Layer

Figure1 shows the elastic layer in the FEM model, containing 2 nodes and 6 degrees of freedom (DOF). The element displacements of each node can be expressed as follows:
(1){δ}elastice=(wiθiui∣wjθjuj)T︸6×1.Figure 1
Configuration of the elastic layer showing the DOF.The stiffness matrix can be derived based on the following energy method:(2)[Ke]elastice=∫01EAl[∂Ne∂ξ]︸6×1T[∂Ne∂ξ]︸1×6dξ+∫01EIl[∂2Nf∂2ξ]T[∂2Nf∂ξ2]dξ︸,6×6
as the shape functions are the following:
(3)[Nf]=[1-3ξ2+2ξ3(ξ-2ξ2+ξ3)l03ξ2-2ξ3(-ξ2+ξ3)l0],[Ne]=[001-ξ00ξ],
in which ξ: the local coordinate, ξ=x/l, ξ∈[0,1],  l: longitudinal length of elastic layer,  A: cross-sectional area of the elastic layer,  E: Young’s modulus of the elastic layer, and  I: moment of inertia of elastic layer.Similarly, the element mass matrix can be expressed as:(4)[M]elastice︸6×6=∫01ml([Nf]︸6×1T[Nf]︸1×6+[Ne]T[Ne])dξ.


## 2.2. FEM Modeling of Fundamental Component II: The Constrained Damping Layer

The FEM model of the constrained layout containing the damping layer sandwiched between two outer layers is shown in Figure2. This Figure illustrates each element consisting of 2 nodes and 8 DOF, where the nodal displacement vector is as follows:
(5){δ}conse=(wiθiu1iu3i∣wjθju1ju3j)T︸8×1.Figure 2
Configuration of the constrained damping layer showing the DOF.Through the introduction of transfer matrix,(6)[TI]=[e1e2e3e5e6e7]T;[TIII]=[e1e2e4e5e6e8]T,
in which each ei means the following vector: ei=(00⋯⋯1︸ithplace⋯⋯0)T; the element elastic stiffness and the element viscoelastic stiffness matrix for this 3-layer component, respectively, are the following:
(7)[Ke]conse=[TI]T·[Ke]elastic·[TI]+[TIII]T·[Ke]elastic·[TIII],[Kv]conse=∫01G2A2lkh2[Ne1-Ne3h2+h0h2·1l·∂Nf1∂ξ]T×[Ne1-Ne3h2+h0h2·1l·∂Nf1∂ξ]dξ,[Ne1]︸1×8=[Ne]︸1×6[TI]︸6×8;[Ne3]=[Ne][TIII],
where A2: cross-sectional area of the damping layer, G2: long-term shear modulus of the damping layer, and k: correction factor of the shear strain energy, for the rectangular cross-section, k=1.2.Also, the element mass matrix for this 3-layer component is(8)[M]conse=[TI]T·[M]elastic·[TI]+[M]cons,2e+[TIII]T·[M]elastic·[TIII],
where
(9)[M]cons,2e=∫01m2l[Nf1]T[Nf1]dξ,[Nf1]=[Nf][TI].


## 2.3. FEM Modeling of a Seven-Layer Constrained Damping Beam

The seven-layer sandwich beam consists of seven alternating layers—four elastic layers and three damping layers. Figure3 shows the FEM model of a seven-layer sandwich beam containing 2 nodes and 10 DOF, and the node displacement vector is as follows:
(10){δ}7layere=(wiθiu1iu3iu5iu7i∣wjθju1ju3ju5ju7j)T︸12×1.Figure 3
Configuration of a seven-layer damping structure showing the DOF.The transfer matrix to obtain the element stiffness and the mass matrix when the 1st, 3rd, 5th, and 7th layers are elastic are follows:(11)[T1]=(e1e2e3e7e8e9)T︸6×12,[T3]=(e1e2e4e7e8e10)T︸6×12,[T5]=(e1e2e5e7e8e11)T︸6×12,[T7]=(e1e2e6e7e8e12)T︸6×12.Similarly, the element stiffness and the mass matrix for the 2nd, 4th, and 6th layers of the constrained damping layer can be derived through the transfer matrix:(12)[T2]=(e1e2e3e4e7e8e9e10)T︸8×12,[T4]=(e1e2e4e5e7e8e10e11)T︸8×12,[T6]=(e1e2e5e6e7e8e11e12)T︸8×12,
where the notation ei means
(13)ei=(00⋯⋯1︸ithplace⋯⋯0)T︸12×1.Based on the above equations and design parameters of each layer, the element mass/stiffness/damping matrix of the seven-layer sandwich damping beam can be expressed as follows:(14)[Ke]e=T1T[Ke1]T1+T3T[Ke3]T3+T5T[Ke5]T5+T7T[Ke7]T7,[Kv]e=T2T[Kv2]T2+T4T[Kv4]T4+T6T[Kv6]T6,[Me]e=∑i=17TiT[Mei]Ti.Thus, the element matrices can be assembled to obtain the global mass/stiffness/damping matrix and can be applied to the boundary condition through the conventional FEM technique. Taking into the consideration of the viscoelastic damping properties, the global matrices need to be manipulated as a portion of the Biot dynamic equation.


## 2.4. Introduction of the Biot Dynamic Equation

To consider the vibration problem numerically, the dynamic equation discretized by FEM technique needs to be expressed by the following second-order ordinary differential equation (ODE) form:(15)Mx¨+Cx˙+Kx=f(t).The Biot viscoelastic damping model numerically represents the complex shear modulus with a series of mini-oscillator perturbing terms:(16)sG~(s)=G∞[1+∑k=1makss+bk],
in which G∞ is the long-term shear moduli; ak and bk are the Biot constants. These parameters are positive and can be determined by nonlinear curve fitting from the experimental data. The curve-fitting procedure will be discussed in Section 3.Substituting the Biot damping model into (15), the dynamic equation with m terms of the Biot parameters for the first viscoelastic material and n terms for the second viscoelastic material can be developed as follows:
(17)G1∞Λv1=Λ1,G2∞Λv2=Λ2,R1=G1∞Rv1Λv1,R2=G2∞Rv2Λv2,
where Rv and Λv are the eigenvector and diagonal eigenvalue matrices, respectively, from the damping matrix C. Additionally, a11⋯a1m, b11⋯b1m, and z11⋯z1m denote m terms of the Biot parameters and the dissipative coordinates, respectively, for first viscoelastic material.Similarly,a21⋯a2n, b21⋯b2n, and z21⋯z2n denote n terms of the Biot parameters and the dissipation coordinates, respectively, for second viscoelastic material. A detailed derivation can be found in the previous publication [28].


## 3. Parametric Determination of the Biot Damping Model

A curve-fitting technique is used to provide the accurate Biot constants to the dynamic equation and to establish the dynamic characteristics of the viscoelastic materials. In this section, the nonlinear curve-fitting procedure for the complex shear modulus in the frequency domain is converted into a nonlinear constrained optimization problem.The complex shear modulus with the Biot damping model form can be broken into real and imaginary parts separately:(18)sG~(jω)=G∞[1+∑i=1Naiϖ2bi2+ϖ2]+jG∞[1+∑i=1Naibiϖbi2+ϖ2].The Biot parameters—G∞, ai, and bi —are estimated from experimental data with the certain fitting frequency range, on real part and imaginary parts separately. Generally speaking, one set of the Biot parameters needs to be determined for each ambient temperature independently. In (18), N is the number of the Biot perturbing items, defining the capability of this numerical approximation. As the Biot terms (N) are increased, the relative error between the experimental data and the curve-fitting result reduces.Assumingx1=G∞; x2=a1; x3=b1; x4=a2; x5=b2; … with the constraint condition: xk≥0; k=1,2,…,num, the target equation of the optimization problem is the following:
(19)minxF(x)=∑j=1P|Gj*(x)-G0j|2.In the target equation (19), G0j stands for the complex shear modulus from the experimental data with P interested points (larger than the number of unknowns). The 3 M ISD-110/112 viscoelastic polymer is selected in this simulation. The experimental data is obtained by the Arrhenius empirical equation from [19]. With a specific fitting range at a particular temperature, the complex shear modulus can be synthesized from one set of the Arrhenius coefficients. The number of terms (N) in (18) needs to be determined to ensure the precision of this approximation. The curve fitting of the experimental data is accomplished using the commercial software package Auto2fit on the real and imaginary parts simultaneously. Using the Biot terms equal to six and four with respect to two commercial damping materials 3 M ISD-110 and 112, respectively, the results are shown in Table 1 for ambient temperature (T) equal to 45°C and frequency range of 500 Hz.Table 1
The Biot constants of 3M ISD-110/112 45°C.
ISD110
ISD112
G
∞
55000 (Pa)
172000 (Pa)
a1
1.809517
5.699386303
a2
14.53095
0.596843249
a3
3.221535
1.000560485
a4
52.01026
0.577694736
a5
19768.22
a6
6.561162
b1
5.410993
4268.18097
b2
1093.778
70.26089968
b3
60.36544
501.5607814
b4
4319.613
1.969150769
b5
2840958
b6
298.0672Figures4(a) and 4(b) show the comparison between the Arrhenius data and curve-fitting data for the real and imaginary parts, respectively. Figure 5 shows the relative error in the fitting range.(a) Comparison between the Arrhenius and curve-fitting data for the real part of the shear modulus (3 M-ISD-110, 45°C). (b) Comparison between the Arrhenius and curve-fitting data for the imaginary part of the shear modulus (3 M-ISD-110, 45°C).
(a)
(b)Figure 5
Relative error between Arrhenius and curve-fitting data.As shown in Figures4(a) and 4(b), the Biot parametric determination technique estimates the dynamic properties of 3 M ISD-110/112 at 45°C with almost zero error. The constants determined using the above procedure along with the FEM model of sandwich beam will now be incorporated to solve the complete Biot dynamic equation using the decoupling transformation technique.


## 4. Decoupling Transformation and Dynamic Solution

In this section, the algorithm used to obtain the frequency response function (FRF) will be discussed with respect to the vibroacoustical problem for a multiple-layer viscoelastic damping structure. In this research, the damping matrixD- in (15) does not have a proportional relationship with the mass and stiffness matrix. Thus, a decoupling transformation is needed to construct the first-order state equation by introducing the auxiliary equation Mq˙-Mq˙=0 as follows:
(20)Ay˙+By=f^,
where
(21)A=[D-M-M-0],B=[K-00-M-],y={qq˙},f^={f-0}.Here,N is the number of DOF in the M-, D-, and K- matrices, the DOF of A and B matrices is 2N.Firstly, the free vibration of (20) will be considered. Assuming f^=0, the following form of solution is obtained: (22a)(Aλ+By)Φ=0,
or
(22b)(Aλ+By){ΨΨλ}=0,where λ matrix stands for 2N complex conjugate eigenvalues including the natural frequencies and loss factors information:
(23)[⋱λ⋱]=[λ1λ1*0⋱0λNλN*].It must be noted that zero items will appear in the eigenvalue matrix if the damping matrixD- is not fully ranked. The mode shape vector Ψ for the vector q can be extracted from the eigenvector matrix Φ with respect to the vector y:
(24)[Φ]={[Ψ]1[Ψ]1*⋯[Ψ]N[Ψ]N*λ1[Ψ]1λ1*[Ψ]1*⋯λN[Ψ]NλN*[Ψ]N*}.In addition, (22b) can be numerically solved by λ{Ψ}=-[A]-1[B]{Ψ} using mathematical software package such as MATLAB or Mathematica.Secondly, the forced vibration solution of (20) in the time domain will be discussed. Assuming f^={F0}, the variable substitution can be made by assuming y=Φp, converting the state-space equation from the time space to the modal space. By left multiplying of ΦTwith the substitution of y, we get:
(25)ΦTAΦp˙+ΦTBΦp=ΦTf^.The diagonal modal mass and stiffness matrix are:(26)ΦTAΦ=[⋱Mp⋱];ΦTBΦ=[⋱Kp⋱].Then rewrite the equation with the diagonal mass and stiffness matrices(27)[⋱Mp⋱]p˙+[⋱Kp⋱]p=ΦTf^.The FRF in the frequency domain can be easily determined through the complex conjugate eigenvalue matrixλ, eigenvector matrix Φ, and the modal mass matrix Mp. The modal scaling factor matrix can be calculated through the following:
(28)[⋱Q⋱]=[⋱Mp⋱]-1.Thus, FRF can be established through the modal parameters, being expressed in partial fraction form in terms of the residue vector and system poles as follows:(29)[H(jω)]=X(jω)F(jω)=∑i=1N[QiΨiΨiT(jω-λi)+QiΨi*Ψi*T(jω-λi*)].The system velocity can now be obtained from the above equation by a simple Fourier transformation. By doing so, the vibration problem can be extended to an FRF-based acoustical problem and the combination of these two analyses is the particle velocities information calculated by the following:(30)X˙(jω)=jω[H(jω)]F(jω)=jω[F(jω)]·∑i=1N[QiΨiΨiT(jω-λi)+QiΨi*Ψi*T(jω-λi*)].


## 5. Acoustical Boundary Element Method (BEM) Analysis

### 5.1. Introduction of Acoustical BEM Theory

In Section4, the vibration problem of the multiple-layer sandwich beam is solved through the time-domain dynamic ordinary differential equation of the Biot damping model with numerical analysis by the FEM technique. The vibration problem can be extended to the acoustical problem by the semicoupled method: the vibration will induce a change in sound pressure, yet the sound pressure will not cause the vibration. In this section, the acoustical interior problem will be numerically solved by 2D boundary element method (BEM) technique [29] in a bounded fluid domain V as shown in Figure 6.Figure 6
Notations of 2-D BEM interior problem in a fluid domainV.The sound pressure distribution (p) of the time-harmonic wave in the domain V satisfies the governing differential equation, well known as the Helmholtz equation, associated with the boundary conditions on boundary Γ(=Γp∪Γvn∪ΓZ) as follows:
(31)∇2p(x)+k2p(x)=0,whenp(x)=p0,x∈Γpvn(x)≡-1jwρ0∂p∂n=vn0,x∈Γvnz(x)≡pvn=z0,x∈ΓZ.Here,k is equal to ω/c, which means that the wave number is equal to the radiant frequency over the speed of sound; vn, ρ0, z stand for the normal velocity, density of the fluid V (normally the air), and acoustical impedance of the fluid V, respectively.In this work, the link between the vibration and the acoustics analysis is the normal velocity at the acoustical boundaries. Recalling the dynamic solution of the decoupling transformation, the particle velocity in the time domain at each node can be calculated through (30) if the multiple-layer sandwich beam is discretized by the FEM; alternatively the FRF, the complex ratio between the output and input response in the frequency domain, can be determined through (29). Once the input signal is given, the particle velocity of the system displacement versus frequency relationship can be conveniently obtained through the FRF.To solve the governing differential equation (31) in the bounded fluid domain V, the Helmholtz Equation can be transformed into the integral equation, converting the 2-D area integration to the 1-D curve integration around the area:
(32)c(ξ)p(ξ)+∫Γ∂Ψ(ξ,x)∂np(ξ)dΓ=-∫Γiρ0ωΨ(ξ,x)vn(x)dΓ.in whichc(ξ): geometry-dependent coefficient, normally c(ξ)=0 when ξ is in the domain V and c(ξ)=0.5 when ξ is on the smooth boundary Γ, p(ξ): sound pressure at source point ξ, Ψ(ξ,x): x is the field point and Ψ=-(i/4)H0(2)(kr) for the 2D BEM problem r the Euclidian distance between x and ξ, H0(2) and the Second-type Henkel function, n: normal vector pointing away to the fluid domain V.By discretizing the boundary into a series of curve-linear elements through the introduction of the shape functions, the integral equation can be calculated numerically by solving the following linear matrix:(33)HP=GVN,
where H comes from the terms of c(ξ) and ∫Γp(ξ)(∂Ψ(ξ,x)/∂n)dΓ, G is derived from -∫Γiρ0ωvn(x)Ψ(ξ,x)dΓ, and the vector P and VN include sound pressure and particle velocity values, both unknowns and known from the boundary condition.Thus, each set of node velocities due to the force input results in one set of solutions on the sound pressure by BEM discussed in this section. In sum, through the proposed acoustical BEM, it is possible to compute the time-harmonic sound pressure distribution corresponding to each single frequency point in the frequency spectrum.


### 5.2. Calculation Details in This BEM Analysis

For this particular acoustical BEM interior problem, the boundary of acoustical cavity is discretized as 18 quadratic equally spaced boundary elements. The quadratic curvilinear element has three nodes, and the interpolation between each node represents the geometry of each element. The shape functions are as folows:(34)N1=12ξ(ξ-1),N2=(ξ+1)(ξ-1),N3=12ξ(ξ+1),
with respect to the following element coordinates:
(35)x=∑i=13xiNi(ξ),y=∑i=13yiNi(ξ),
where xi and yi are the coordinates at each nodal point, and ξ stands for the local coordinate between −1 and 1 on a master element.When the seven-layered sandwich beam (L=1 m) is simply supported at the bottom of the acoustical cavity, the sound pressure level at the field point (x=0.5 m, y=0.4 m) is calculated through this proposed method, and the calculation results are presented in Section 6. Figure 7 demonstrates the detailed layout of this 2D acoustical cavity problem. The anechoic boundary condition is applied on the inside of the acoustical cavity, and the thickness of the multiple-layered beam is neglected.Figure 7
Layout of BEM acoustical cavity problem.


## 5.1. Introduction of Acoustical BEM Theory

In Section4, the vibration problem of the multiple-layer sandwich beam is solved through the time-domain dynamic ordinary differential equation of the Biot damping model with numerical analysis by the FEM technique. The vibration problem can be extended to the acoustical problem by the semicoupled method: the vibration will induce a change in sound pressure, yet the sound pressure will not cause the vibration. In this section, the acoustical interior problem will be numerically solved by 2D boundary element method (BEM) technique [29] in a bounded fluid domain V as shown in Figure 6.Figure 6
Notations of 2-D BEM interior problem in a fluid domainV.The sound pressure distribution (p) of the time-harmonic wave in the domain V satisfies the governing differential equation, well known as the Helmholtz equation, associated with the boundary conditions on boundary Γ(=Γp∪Γvn∪ΓZ) as follows:
(31)∇2p(x)+k2p(x)=0,whenp(x)=p0,x∈Γpvn(x)≡-1jwρ0∂p∂n=vn0,x∈Γvnz(x)≡pvn=z0,x∈ΓZ.Here,k is equal to ω/c, which means that the wave number is equal to the radiant frequency over the speed of sound; vn, ρ0, z stand for the normal velocity, density of the fluid V (normally the air), and acoustical impedance of the fluid V, respectively.In this work, the link between the vibration and the acoustics analysis is the normal velocity at the acoustical boundaries. Recalling the dynamic solution of the decoupling transformation, the particle velocity in the time domain at each node can be calculated through (30) if the multiple-layer sandwich beam is discretized by the FEM; alternatively the FRF, the complex ratio between the output and input response in the frequency domain, can be determined through (29). Once the input signal is given, the particle velocity of the system displacement versus frequency relationship can be conveniently obtained through the FRF.To solve the governing differential equation (31) in the bounded fluid domain V, the Helmholtz Equation can be transformed into the integral equation, converting the 2-D area integration to the 1-D curve integration around the area:
(32)c(ξ)p(ξ)+∫Γ∂Ψ(ξ,x)∂np(ξ)dΓ=-∫Γiρ0ωΨ(ξ,x)vn(x)dΓ.in whichc(ξ): geometry-dependent coefficient, normally c(ξ)=0 when ξ is in the domain V and c(ξ)=0.5 when ξ is on the smooth boundary Γ, p(ξ): sound pressure at source point ξ, Ψ(ξ,x): x is the field point and Ψ=-(i/4)H0(2)(kr) for the 2D BEM problem r the Euclidian distance between x and ξ, H0(2) and the Second-type Henkel function, n: normal vector pointing away to the fluid domain V.By discretizing the boundary into a series of curve-linear elements through the introduction of the shape functions, the integral equation can be calculated numerically by solving the following linear matrix:(33)HP=GVN,
where H comes from the terms of c(ξ) and ∫Γp(ξ)(∂Ψ(ξ,x)/∂n)dΓ, G is derived from -∫Γiρ0ωvn(x)Ψ(ξ,x)dΓ, and the vector P and VN include sound pressure and particle velocity values, both unknowns and known from the boundary condition.Thus, each set of node velocities due to the force input results in one set of solutions on the sound pressure by BEM discussed in this section. In sum, through the proposed acoustical BEM, it is possible to compute the time-harmonic sound pressure distribution corresponding to each single frequency point in the frequency spectrum.


## 5.2. Calculation Details in This BEM Analysis

For this particular acoustical BEM interior problem, the boundary of acoustical cavity is discretized as 18 quadratic equally spaced boundary elements. The quadratic curvilinear element has three nodes, and the interpolation between each node represents the geometry of each element. The shape functions are as folows:(34)N1=12ξ(ξ-1),N2=(ξ+1)(ξ-1),N3=12ξ(ξ+1),
with respect to the following element coordinates:
(35)x=∑i=13xiNi(ξ),y=∑i=13yiNi(ξ),
where xi and yi are the coordinates at each nodal point, and ξ stands for the local coordinate between −1 and 1 on a master element.When the seven-layered sandwich beam (L=1 m) is simply supported at the bottom of the acoustical cavity, the sound pressure level at the field point (x=0.5 m, y=0.4 m) is calculated through this proposed method, and the calculation results are presented in Section 6. Figure 7 demonstrates the detailed layout of this 2D acoustical cavity problem. The anechoic boundary condition is applied on the inside of the acoustical cavity, and the thickness of the multiple-layered beam is neglected.Figure 7
Layout of BEM acoustical cavity problem.


## 6. Numerical Results and Discussion

### 6.1. Design Parameter of Sandwich Beam and Vibration Analysis Result

In Figure8, a seven-layer sandwich beam with viscoelastic cores is shown, with the design parameters listed in Table 2.Table 2
Design parameters of seven-layer structure.
Length: 1 m
Thickness: 0.1 m
Number of element: 12
Number of nodes: 13
Number of layer
Height of layer
Elastic/viscoelastic properties
Material density
1st
h
1
=
1 mm
E
1
=
210 GPa
ρ
1
=
7800 kg/m3
2nd
h
2
=
0.8 mm
G
2: Biot
ρ
2
=
970 kg/m3
3rd
h
3
=
1 mm
E
3
=
210 GPa
ρ
3
=
7800 kg/m3
4th
h
4
=
0.8 mm
G
4: Biot
ρ
4
=
970 kg/m3
5th
h
5
=
1 mm
E
5
=
210 GPa
ρ
5
=
7800 kg/m3
6th
h
6
=
0.8 mm
G
6: Biot
ρ
6
=
970 kg/m3
7th
h
7
=
1 mm
E
7
=
210 GPa
ρ
7
=
7800 kg/m3Figure 8
Seven-layer sandwich structure with viscoelastic cores.The data presented in Table2 are used to predict the vibration performance of the system using the numerical simulation method presented in this paper, and the results are compared with the closed-form solution of Hao [19]. The curve-fitting results for the damping material 3M ISD-110 at 45°C discussed earlier are selected for the shear modulus of the viscoelastic layers in this example. The results are shown in Table 3. It shows that the simulation presented in this paper conforms to the closed-form solution. This validates the analysis methodology proposed in the paper.Table 3
Comparison of results for simply supported boundary condition.
Number of  mode
Damping model
Hao [19]
FEM model of this paper
ISD110-45°C Arrhenius
ISD110-45°C6-term Biot
1st
Frequency
4.7443 Hz
4.5834 Hz
Loss factor
0.6248
0.7916
2nd
Frequency
13.902 Hz
13.9489 Hz
Loss factor
0.6008
0.6824
3rd
Frequency
27.661 Hz
27.7721 Hz
Loss factor
0.5317
0.5632
4th
Frequency
46.1548 Hz
47.1053 Hz
Loss factor
0.4715
0.4681
5th
Frequency
69.3118 Hz
68.82259 Hz
Loss factor
0.421
0.4226


### 6.2. Frequency-Spectrum Analysis under the Arbitrary Input

Figure9 shows the transverse velocity of the middle node (node number 7) with a 10 N step input in the frequency domain vertically applied at the middle (node number 7) of the simply-supported seven-layer sandwich beam with the same design parameters as the previous example. The same curve-fitting results of 3 M ISD110 at the ambient temperature of 45°C for the shear modulus are used in this example. This pivotal result is the demonstration of extending the vibration to the acoustical problem in the frequency domain when an arbitrary force is applied on the structure.Figure 9
Transverse velocity under the impulse excitation (on middle node #7).


### 6.3. Acoustical BEM Results

Figure10 illustrates the contour plot (f=10 Hz) of SPL when the seven-layer sandwich beam (using the same design parameters as before) is subjected to a 10 N step input in the frequency domain at the middle node.Figure 10
Contour Plot of Sound Pressure Level (in dB) when the impulse force applied isN.The interpolation of each elements result in Figures10 and 11 shows the continuous sound pressure distribution in the acoustical cavity with an anechoic boundary condition. Figure 12 extracts the frequency spectrum of SPL at the filed point (0.5, 0.4 m) indicated by red dot in Figure 10. From the results of Figure 12, it can be found that the dominant contribution is due to the peak value of the first flexible vibration mode, which is in agreement with the frequency-spectrum analysis of the vibration problem.Figure 11
Element result of Sound Pressure Level in dB (10 Hz).Figure 12
Nodal frequency spectrum at field point (0.5, 0.4 m).


### 6.4. Validation Using a BEM Commercial Software

In this section, a hybrid FEM-BEM model of a beam without the viscoelastic damping was developed using the commercial software packages ANSYS ADPL and LMS Virtual. Lab Acoustics. The harmonic vibration analysis is conducted in ANSYS APDL module, and the frequency spectrum of field point SPL was calculated in Virtual.Lab Acoustics module for comparison with the SPL frequency spectrum presented in Section5. The analysis sequence consists of the following steps.(a)
Build the FEM model and apply appropriate boundary conditions in ANSYS ADPL. The 8-node element SOLID45 (element size  =10 mm for each layer) was used to build the 3D seven-layer model. The design parameters are identical with the parameters in Tables1 and 2 for the comparison and the geometry boundary conditions are simply supported. A 10 N force at each frequency is applied at the middle nodes.(b)
Conduct the harmonic vibration analysis in ANSYS ADPL. The harmonic analysis is used to calculate the nodal displacements for a forced vibration problem in the frequency domain. The frequency range is 0–200 Hz with a 2 Hz for step size, and the full method is being utilized in this analysis. The comparison of system frequencies between ANSYS modal results and calculation results by the Biot dynamic equation is shown in Table4. The results show that the 3D model built in ANSYS APDL has good correlation with the FEM model.(c)
Prepare the BEM mesh in LMS Virtual. Lab Pre-Acoustics module. It converts from a solid FEM model to a skin mesh that the BEM analysis requires. The BEM mesh, can be seen as a wrap around the structural mesh and usually the BEM mesh is coarser.(d)
Calculate the sound pressure in LMS Virtual.lab Acoustics module. Both acoustical and structural meshes are imported to VL Acoustics. The nodal displacement at each vibration mode calculated in ANSYS APDL is also imported and mesh-mapped to the acoustical skin mesh as the vibration boundary condition. The location of field plane and field point is consistent with the 2D BEM analysis in this research. The acoustical pressure is solved over the frequency range from 2 to 200 Hz.Table 4
Comparison of modal results with ANSYS simulation.
Number of  mode
Damping model
ANSYS 3D
FEM model of this research
No damping
3M-ISD110-45°C6-term Biot
1st
Frequency
6.4803 Hz
4.5834 Hz
Loss factor
0.7916
2nd
Frequency
14.338 Hz
13.9489 Hz
Loss factor
0.6824
3rd
Frequency
27.989 Hz
27.7721 Hz
Loss factor
0.5632As shown in Figure13, the peak frequency from the 2-D BEM calculation matches with the first dominant SPL peak obtained from the VL Acoustic result without the damping. Comparing the two results, it is clear that the introduction of viscoelastic damping not only causes almost a 20 dB reduction in the first peak SPL but also attenuates the sound at other peaks as well. This proves that the use of viscoelastic damping material will greatly attenuate the vibroacoustical response of the structure.Figure 13
Comparison of SPL between VL and 2D BEM calculation at field point (0.5, 0.4), 10 N input.


### 6.5. Acoustical Performance for a Combination of Several Viscoelastic Materials at Different Temperatures

The temperature is a significant external factor affecting the performance of viscoelastic damping material in a mechanical system. With an increase in temperature, the loss factor approaches its best performance towards the transition region and then decreases afterwards. In this example, the objective is to study the effects of both 3 M ISD110 material (that has a better damping performance) and the 3 M ISD112 over the chosen temperature between 40 and 60 degree Celsius. It is of interest to study the effect of the combination of these two materials on the damping of the structure.To introduce the different viscoelastic materials, the seven-layer sandwich beam (with the same parameters as in the previous example) is redesigned incorporating both damping materials (3 M ISD110 and ISD112). This system is compared to an identical structure with only one damping material (either 3 M ISD110 or ISD112). In the system including two viscoelastic materials, the outer damping layers (2nd and 6th) are 3 M ISD112 and the inner damping layer (5th) is the 3 M ISD110. The simply supported boundary condition is examined in this numerical example, and the temperature range is from 40 to 60 degree Celsius. The acoustical response is also calculated with the step input in the frequency domain (equivalent to impulse input in the time domain). Table5 shows the first order natural frequency, the system loss factor, and the corresponding peak value (dB) of the sound pressure level over the temperature range with the simply supported boundary condition applied to the FEM model.Table 5
First damped frequency, loss factor and SPL (simple supported BC).
Type
Temp-°C
Freq-rad/s
Loss factor
SPL-dB
3M ISD-110
40
33.17
1.085
50.35
45
28.82
0.7899
51.91
50
26.25
0.5847
54.15
55
24.38
0.464
56.24
60
23.5
0.387
57.45
3M ISD-112
40
32.65
0.1946
59.54
45
31.62
0.1652
61.25
50
30.8
0.1413
62.30
55
30.15
0.1218
62.90
60
29.61
0.1057
63.21
3M ISD 110 & 112
40
32.88
0.4427
54.31
45
30.76
0.3506
56.45
50
29.3
0.2776
57.68
55
28.22
0.2252
58.05
60
27.54
0.185
59.67It can be seen that for the same damping material, as the ambient temperature is increased, the value of SPL increases while the loss factor decreases.


## 6.1. Design Parameter of Sandwich Beam and Vibration Analysis Result

In Figure8, a seven-layer sandwich beam with viscoelastic cores is shown, with the design parameters listed in Table 2.Table 2
Design parameters of seven-layer structure.
Length: 1 m
Thickness: 0.1 m
Number of element: 12
Number of nodes: 13
Number of layer
Height of layer
Elastic/viscoelastic properties
Material density
1st
h
1
=
1 mm
E
1
=
210 GPa
ρ
1
=
7800 kg/m3
2nd
h
2
=
0.8 mm
G
2: Biot
ρ
2
=
970 kg/m3
3rd
h
3
=
1 mm
E
3
=
210 GPa
ρ
3
=
7800 kg/m3
4th
h
4
=
0.8 mm
G
4: Biot
ρ
4
=
970 kg/m3
5th
h
5
=
1 mm
E
5
=
210 GPa
ρ
5
=
7800 kg/m3
6th
h
6
=
0.8 mm
G
6: Biot
ρ
6
=
970 kg/m3
7th
h
7
=
1 mm
E
7
=
210 GPa
ρ
7
=
7800 kg/m3Figure 8
Seven-layer sandwich structure with viscoelastic cores.The data presented in Table2 are used to predict the vibration performance of the system using the numerical simulation method presented in this paper, and the results are compared with the closed-form solution of Hao [19]. The curve-fitting results for the damping material 3M ISD-110 at 45°C discussed earlier are selected for the shear modulus of the viscoelastic layers in this example. The results are shown in Table 3. It shows that the simulation presented in this paper conforms to the closed-form solution. This validates the analysis methodology proposed in the paper.Table 3
Comparison of results for simply supported boundary condition.
Number of  mode
Damping model
Hao [19]
FEM model of this paper
ISD110-45°C Arrhenius
ISD110-45°C6-term Biot
1st
Frequency
4.7443 Hz
4.5834 Hz
Loss factor
0.6248
0.7916
2nd
Frequency
13.902 Hz
13.9489 Hz
Loss factor
0.6008
0.6824
3rd
Frequency
27.661 Hz
27.7721 Hz
Loss factor
0.5317
0.5632
4th
Frequency
46.1548 Hz
47.1053 Hz
Loss factor
0.4715
0.4681
5th
Frequency
69.3118 Hz
68.82259 Hz
Loss factor
0.421
0.4226


## 6.2. Frequency-Spectrum Analysis under the Arbitrary Input

Figure9 shows the transverse velocity of the middle node (node number 7) with a 10 N step input in the frequency domain vertically applied at the middle (node number 7) of the simply-supported seven-layer sandwich beam with the same design parameters as the previous example. The same curve-fitting results of 3 M ISD110 at the ambient temperature of 45°C for the shear modulus are used in this example. This pivotal result is the demonstration of extending the vibration to the acoustical problem in the frequency domain when an arbitrary force is applied on the structure.Figure 9
Transverse velocity under the impulse excitation (on middle node #7).


## 6.3. Acoustical BEM Results

Figure10 illustrates the contour plot (f=10 Hz) of SPL when the seven-layer sandwich beam (using the same design parameters as before) is subjected to a 10 N step input in the frequency domain at the middle node.Figure 10
Contour Plot of Sound Pressure Level (in dB) when the impulse force applied isN.The interpolation of each elements result in Figures10 and 11 shows the continuous sound pressure distribution in the acoustical cavity with an anechoic boundary condition. Figure 12 extracts the frequency spectrum of SPL at the filed point (0.5, 0.4 m) indicated by red dot in Figure 10. From the results of Figure 12, it can be found that the dominant contribution is due to the peak value of the first flexible vibration mode, which is in agreement with the frequency-spectrum analysis of the vibration problem.Figure 11
Element result of Sound Pressure Level in dB (10 Hz).Figure 12
Nodal frequency spectrum at field point (0.5, 0.4 m).


## 6.4. Validation Using a BEM Commercial Software

In this section, a hybrid FEM-BEM model of a beam without the viscoelastic damping was developed using the commercial software packages ANSYS ADPL and LMS Virtual. Lab Acoustics. The harmonic vibration analysis is conducted in ANSYS APDL module, and the frequency spectrum of field point SPL was calculated in Virtual.Lab Acoustics module for comparison with the SPL frequency spectrum presented in Section5. The analysis sequence consists of the following steps.(a)
Build the FEM model and apply appropriate boundary conditions in ANSYS ADPL. The 8-node element SOLID45 (element size  =10 mm for each layer) was used to build the 3D seven-layer model. The design parameters are identical with the parameters in Tables1 and 2 for the comparison and the geometry boundary conditions are simply supported. A 10 N force at each frequency is applied at the middle nodes.(b)
Conduct the harmonic vibration analysis in ANSYS ADPL. The harmonic analysis is used to calculate the nodal displacements for a forced vibration problem in the frequency domain. The frequency range is 0–200 Hz with a 2 Hz for step size, and the full method is being utilized in this analysis. The comparison of system frequencies between ANSYS modal results and calculation results by the Biot dynamic equation is shown in Table4. The results show that the 3D model built in ANSYS APDL has good correlation with the FEM model.(c)
Prepare the BEM mesh in LMS Virtual. Lab Pre-Acoustics module. It converts from a solid FEM model to a skin mesh that the BEM analysis requires. The BEM mesh, can be seen as a wrap around the structural mesh and usually the BEM mesh is coarser.(d)
Calculate the sound pressure in LMS Virtual.lab Acoustics module. Both acoustical and structural meshes are imported to VL Acoustics. The nodal displacement at each vibration mode calculated in ANSYS APDL is also imported and mesh-mapped to the acoustical skin mesh as the vibration boundary condition. The location of field plane and field point is consistent with the 2D BEM analysis in this research. The acoustical pressure is solved over the frequency range from 2 to 200 Hz.Table 4
Comparison of modal results with ANSYS simulation.
Number of  mode
Damping model
ANSYS 3D
FEM model of this research
No damping
3M-ISD110-45°C6-term Biot
1st
Frequency
6.4803 Hz
4.5834 Hz
Loss factor
0.7916
2nd
Frequency
14.338 Hz
13.9489 Hz
Loss factor
0.6824
3rd
Frequency
27.989 Hz
27.7721 Hz
Loss factor
0.5632As shown in Figure13, the peak frequency from the 2-D BEM calculation matches with the first dominant SPL peak obtained from the VL Acoustic result without the damping. Comparing the two results, it is clear that the introduction of viscoelastic damping not only causes almost a 20 dB reduction in the first peak SPL but also attenuates the sound at other peaks as well. This proves that the use of viscoelastic damping material will greatly attenuate the vibroacoustical response of the structure.Figure 13
Comparison of SPL between VL and 2D BEM calculation at field point (0.5, 0.4), 10 N input.


## 6.5. Acoustical Performance for a Combination of Several Viscoelastic Materials at Different Temperatures

The temperature is a significant external factor affecting the performance of viscoelastic damping material in a mechanical system. With an increase in temperature, the loss factor approaches its best performance towards the transition region and then decreases afterwards. In this example, the objective is to study the effects of both 3 M ISD110 material (that has a better damping performance) and the 3 M ISD112 over the chosen temperature between 40 and 60 degree Celsius. It is of interest to study the effect of the combination of these two materials on the damping of the structure.To introduce the different viscoelastic materials, the seven-layer sandwich beam (with the same parameters as in the previous example) is redesigned incorporating both damping materials (3 M ISD110 and ISD112). This system is compared to an identical structure with only one damping material (either 3 M ISD110 or ISD112). In the system including two viscoelastic materials, the outer damping layers (2nd and 6th) are 3 M ISD112 and the inner damping layer (5th) is the 3 M ISD110. The simply supported boundary condition is examined in this numerical example, and the temperature range is from 40 to 60 degree Celsius. The acoustical response is also calculated with the step input in the frequency domain (equivalent to impulse input in the time domain). Table5 shows the first order natural frequency, the system loss factor, and the corresponding peak value (dB) of the sound pressure level over the temperature range with the simply supported boundary condition applied to the FEM model.Table 5
First damped frequency, loss factor and SPL (simple supported BC).
Type
Temp-°C
Freq-rad/s
Loss factor
SPL-dB
3M ISD-110
40
33.17
1.085
50.35
45
28.82
0.7899
51.91
50
26.25
0.5847
54.15
55
24.38
0.464
56.24
60
23.5
0.387
57.45
3M ISD-112
40
32.65
0.1946
59.54
45
31.62
0.1652
61.25
50
30.8
0.1413
62.30
55
30.15
0.1218
62.90
60
29.61
0.1057
63.21
3M ISD 110 & 112
40
32.88
0.4427
54.31
45
30.76
0.3506
56.45
50
29.3
0.2776
57.68
55
28.22
0.2252
58.05
60
27.54
0.185
59.67It can be seen that for the same damping material, as the ambient temperature is increased, the value of SPL increases while the loss factor decreases.


## 7. Conclusions

A framework for conducting vibro-acoustical analysis for multiple-layer beam structures containing different types of viscoelastic materials is presented in this paper. Several observations and conclusions can be drawn from the results of this research.(1)
The vibration section of the proposed analysis consists of FEM model of multiple-layered damping beam incorporating the Biot damping model. The FEM model of the beam structure can be extended to more complicated damping structures using the same procedure. The nonlinear curve-fitting technique accurately estimates the Biot constants. The Biot damping model can then be solved using the decoupling transformation to yield the frequency-spectrum analysis.(2)
The Biot damping model is also capable of improving a structure’s damping performance by adding new features such as different viscoelastic materials and the variation of operating temperature. The result obtained through the procedure of vibration analysis discussed in this paper compares well to the closed-form solution from a previous work. The first peak from the frequency spectrum is the predominant cause of the vibration issue in this damping structure.(3)
The direct boundary element method of analysis for acoustical cavity applied under anechoic boundaries was chosen as the basis for predicting the particle velocity from the frequency-spectrum analysis. The acoustical result validates the frequency-spectrum result from vibration analysis and has good agreement with the predicted SPL spectrum of the identical sandwich beam without damping calculated by commercial software.



---

*Source: 645232-2013-03-30.xml*