Dynamic Transient Response of Viscoelastic Structures

Abstract

This chapter focuses on investigating the dynamic transient response of viscoelastic structures. First, the influence of nonviscous modes on the vibrational response of exponentially damped systems has been studied on lumped parameter systems where the analytical solution has been derived by modal superposition and by means of Laplace transformation. Then, the analytical solution is obtained by modal superposition and compared to two numerical solutions derived for continuous systems by finite element formulations. These numerical solutions have been solved by modal superposition and by direct integration applying through a particularly built method together with the Newmark method. Finally, an experimental procedure for studying the influence that geometrical properties of viscoelastic joints have on the vibrational response of a metallic beam doubly supported on viscoelastic adhesive joints has been developed.

1. Introduction

This chapter focuses on investigating the dynamic transient response of viscoelastic structures. Mechanical behavior of viscoelastic materials is not only related to the instantaneous stress, but it is also a consequence of the past history of the stress. If dissipative forces in a structural system arise from viscoelastic materials, the classical governing equations of structural dynamics is not reached and direct methods or internal variables should be employed. Nevertheless, direct methods do not provide any information about the contribution of each vibration mode, information of main importance for engineering applications.

For the point of view of the practical application, viscoelastic adhesive joints are used in structural noise control due to its capability to introduce effective modal damping below 1 kHz. In structures under dynamic loads, the transmission of noise and vibration is governed by joint behavior.

Hence, an experimental procedure is presented to analyze the influence of geometry of viscoelastic joint over the dynamic response of the low-order flexural modes of adhesively bonded beams.

2. Analysis of exponentially damped systems

This chapter is aimed at investigating the influence of nonviscous modes on vibrational response of viscoelastic systems. Thus, exponential damping models are considered.

2.1. Lumped parameters systems

In short, the main objective of the present section is to study the influence of nonviscous modes on the vibrational response of exponentially damped systems. The analysis is carried out in time domain over a single degree of freedom (dof). The conditions for the equivalence between Zener and exponential damping models are stated. Next, the time response for free vibration is obtained. On the one hand, the need for solving internal variables is avoided through the Laplace transform properties. On the other hand, internal variables are used, and modal superposition is applied.

Figure 1.

Single dof system: (a) using a Zener model, (b) using relaxation functions.

2.1.1. Analysis via a single degree-of-freedom system

The equations of motion for the free vibration of a single dof system using a Zener and an exponential model are obtained in this section. The conditions for the full equivalence between both models are also deduced. The analytical response is reached by two ways: first, the internal variables are eliminated and the transient response is computed by means of Laplace transformation; and second, the internal variables are kept, allowing to apply traditional modal superposition.

2.1.1.1. Solution without solving internal variables

Next, the single dof systems represented in Figure 1 are analyzed, whose dissipative forces are modeled by means of Zener and exponential damping models, respectively.

The equation of motion for the lumped mass m results in

$m{\ddot{u}}_{\mathrm{z}}+c_1\left({\stackrel{.}{u}}_{\mathrm{z}}-\stackrel{.}{y}\right)+k_0{u}_{\mathrm{z}}=0$E1

where $c_1$ denotes the damping coefficient, $k_0$ the stiffness of the parallel spring, u_z represents the displacement of the mass, and y is an internal variable needed to solve the problem. This internal variable represents the displacement of the connection point between the dashpot and the in-series spring. In fact, it is necessary to write the force equilibrium equation for the internal variable y, yielding

$k_{1}y=c_1\left({\stackrel{.}{u}}_{\mathrm{z}}-\stackrel{\text{.}}{y}\right)$E2

where $k_1$ denotes the stiffness of the spring in series with the damper. Combining Eqs. (1) and (2) the linear differential equation for the displacement u_z is obtained as

$\frac{c_1}{k_1}m{\stackrel{⃛}{u}}_{\mathrm{z}}+m{\ddot{u}}_{\mathrm{z}}+c_1\left(1+\frac{k_0}{k_1}\right){\stackrel{.}{u}}_{\mathrm{z}}+k_0{u}_{\mathrm{z}}=0$E3

To solve this third-order differential equation, the initial displacement $u_{\mathrm{z}}\left(0\right)=u_0$, the initial velocity

, and the initial acceleration ${\ddot{u}}_{\mathrm{z}}\left(0\right)={\ddot{u}}_0$ are needed. (In structural dynamics just initial displacement and the initial velocity are needed.) However, this third condition (initial acceleration ${\ddot{u}}_{\mathrm{z}}\left(0\right)={\ddot{u}}_0$) can be obtained as a function of the initial displacement.

One way to solve Eq. (3) is making use of the Laplace transform properties, the response ū_z(s) being

${\overline{u}}_{\mathrm{z}}\left(s\right)=\frac{\frac{c_1}{k_1}m{u}_0{s}^2+m\left(\frac{c_1}{k_1}{\stackrel{.}{u}}_0+u_0\right)s+m{\stackrel{.}{u}}_0+c_1{u}_0+\frac{c_1}{k_1}\left(m{\ddot{u}}_0+k_0{\stackrel{.}{u}}_0\right)}{\frac{c_1}{k_1}m{s}^3+m{s}^2+c_1\left(1+\frac{k_0}{k_1}\right)s+k_0}$E4

where ū_z(s) represents the Laplace transform for the time response u_z(t). Consequently, this time response u_z(t) can be obtained by means of the inverse Laplace transform of Eq. (4).

However, the conditions for the equivalence between both considered damping models are deduced first. Regarding the exponential model, the motion equation [1] for the 1 dof system yields

$m{\ddot{u}}_{\mathrm{e}}\left(t\right)+{\displaystyle {\int }_0^{t}c\mu {\mathrm{e}}^{-\mu \left(t-\tau \right)}{\stackrel{.}{u}}_{\mathrm{e}}\left(\tau \right)\mathrm{d}\tau }+k_0{u}_{\mathrm{e}}\left(t\right)=0$E5

where u_e represents the displacement for the exponential formulation and c the damping coefficient. In contrast to the previous case, to solve Eq. (5) only two initial conditions are needed, the initial displacement $u_{\mathrm{e}}\left(0\right)=u_0$ and the initial velocity

.

Accordingly to the Zener formulation, by transforming Eq. (5) into the Laplace domain, the Laplace transform of the displacement ū_e(s) satisfies

${\overline{u}}_{\mathrm{e}}\left(s\right)=\frac{\frac{1}{\mu }m{u}_0{s}^2+m\left(\frac{1}{\mu }{\stackrel{.}{u}}_0+u_0\right)s+m{\stackrel{.}{u}}_0+c{u}_0}{\frac{1}{\mu }m{s}^3+m{s}^2+c\left(1+\frac{k_0}{c\mu }\right)s+k_0}$E6

At this point, by comparing Eqs. (4) and (6), a useful observation can be made. These two considered damping models are fully equivalents if three conditions are fulfilled. The first condition to be hold is that the damping coefficients of both models must be the same, $c_1=c$. The second condition is that the relaxation parameter μ of the exponential model and the parameters $k_1$ and $c_1$ of the Zener one are related according to $\mu =k_1/c_1$.

And finally, the third condition for the fully equivalence between both formulations results from the relationship between initial displacement and acceleration

${\ddot{u}}_0=-\frac{k_0}{m}{u}_0$E7

implying that in reality, only two independent initial conditions must be taken into account. Having verified that both formulations are equivalent, from now on, a unique response u(t) = u_z(t) = u_e(t) is considered. Thus, the transient response u(t) can be obtained by means of the inverse Laplace transform of Eq. (4) or (6), either, resulting in

$u\left(t\right)={\displaystyle \sum _{i=1}^3\left\{\frac{u_0\left({s_i}^{2}m+c\mu \right)+m\left({\stackrel{.}{u}}_0\mu +u_0{s}_i\mu +{\stackrel{.}{u}}_0{s}_i\right)}{3m{s_i}^2+2sm\mu +c\mu +k_0}{\mathrm{e}}^{s_{i}t}\right\}}$E8

where s_i denotes the roots of the third-order characteristic equation

$s^3+\mu s^2+\frac{c\mu +k_0}{m}s+\frac{k_0\mu }{m}=0$E9

These three roots may be solved by means of the Cardan method [2], yielding

$s_{1,2}^{\ast }=-\frac{c}{3}-\frac{1}{2}\left(u+v\right)\pm \mathrm{i}\frac{\sqrt[3]{3}}{2}\left(u-v\right)$E10

and

$s_3=-\frac{c}{3}+\left(u+v\right)$E11

where

, the parameters p and q being $p=\frac{3B-C^2}{9}$ and $q=\frac{9CB-27A-2C^3}{54}$, respectively, with A = μ, $B=\frac{1}{m}\left(c\mu +k_0\right)$, and $C=\frac{k_0\mu }{m}$.

It should be pointed out that, in contrast to a viscous 1 dof system, three roots have been found instead of two, the third extra root being always real, involving an overdamped vibration mode. This fact implies that the system could not oscillate even if the roots $s_{1,2}^{\ast }$ are complex (this question was extensively studied by Muller [3] and Adhikari [4]).

2.1.1.2. Solution solving internal variables

Next, the time response for the system modeled in Figure 1(a) is obtained by means of classical modal superposition, using the internal variable y(t). Rewriting Eqs. (1) and (2) in matrix form, a classical second-order differential equation for free vibration

$\mathbf{M}\ddot{\mathbf{z}}\left(t\right)+\mathbf{C}\stackrel{.}{\mathbf{z}}\left(t\right)+\mathbf{K}\mathbf{z}\left(t\right)=\mathbf{0}$E12

is reached, where M, C, and K are the mass, damping, and stiffness matrices, given by

$\mathbf{M}=\left[\begin{array}{cc}\hfill m\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right];\mathbf{C}=\left[\begin{array}{cc}\hfill c_1\hfill & \hfill -c_1\hfill \\ \hfill -c_1\hfill & \hfill c_1\hfill \end{array}\right];\mathbf{K}=\left[\begin{array}{cc}\hfill k_0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill k_1\hfill \end{array}\right]$E13

where z(t) vector satisfies $\mathbf{z}\left(t\right)={\left\{\begin{array}{cc}\hfill u\left(t\right)\hfill & \hfill y\left(t\right)\hfill \end{array}\right\}}^{\mathrm{T}}$.

Solving the eigenproblem for Eq. (12), the eigenvalues are obtained from the characteristic equation

$s^3+\frac{k_1}{c_1}{s}^2+\frac{1}{m}\left(k_1+k_0\right)s+\frac{k_1{k}_0}{c_{1}m}=0$E14

As it was expected, Eqs. (9) and (14) are equivalent, and the same three eigenvalues $s_1^{*}$, $s_2^{*}$, and $s_3$ given by Eqs. (9) and (10) are obtained. Associated to the ith eigenvalue s_i, the ith eigenvector Z_i satisfies

${\mathbf{Z}}_i=\left\{\begin{array}{c}\hfill U_i\hfill \\ \hfill Y_i\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill \frac{1}{1+\frac{k_1}{s_i{c}_1}}\hfill \end{array}\right\}$E15

Thus, applying modal superposition, the time response z(t) can be written as

$\mathbf{z}\left(t\right)=q_1{\mathbf{Z}}_1{\mathrm{e}}^{s_1^{\ast }t}+q_2{\mathbf{Z}}_2{\mathrm{e}}^{s_2^{\ast }t}+q_3{\mathbf{Z}}_3{\mathrm{e}}^{s_{3}t}$E16

where $q_1$, $q_2$, and $q_3$ denote the modal participation factors. To solve them, initial conditions $\mathbf{z}\left(0\right)$ and $\stackrel{.}{\mathbf{z}}\left(0\right)$ have to be employed. Nevertheless, to establish the initial value for the internal variable $y\left(0\right)$ and its derivative

, some physical assumptions can be made. In order to respect the force equilibrium at the initial instant, the value of the internal variable $y\left(0\right)$ must be zero and its time derivative

must be the same as the initial velocity

. Hence, the internal variable does not introduce extra energy into the system: the spring of stiffness $k_1$ is not deformed and the dashpot presents a rigid movement dissipating no energy. Therefore, the initial conditions are given by $z(0)=z_0={\left\{\begin{array}{cc}u(0)& y(0)\end{array}\right\}}^{\text{T}}={\left\{\begin{array}{cc}{u}_0& 0\end{array}\right\}}^{\text{T}}$ and $\dot{z}(0)={\dot{z}}_0={\left\{\begin{array}{cc}\dot{u}(0)& \dot{y}(0)\end{array}\right\}}^{\text{T}}={\left\{\begin{array}{cc}{\dot{u}}_0& {\dot{u}}_0\end{array}\right\}}^{\text{T}}$.

However, applying these initial conditions, a four equation system with three unknowns is reached. Nevertheless, the system has a unique solution because the rank of the resulting system matrix is 3. This fact is due to the linear combination among the internal variable y(t), its time derivative

, and the velocity

according to Eq. (2). Thus, the three equations needed to solve the modal participation factors yield

$\left\{\begin{array}{c}\hfill u_0\hfill \\ \hfill {\stackrel{.}{u}}_0\hfill \\ \hfill 0\hfill \end{array}\right\}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill s_1^{\ast }\hfill & \hfill s_2^{\ast }\hfill & \hfill s_3\hfill \\ \hfill Y_1\hfill & \hfill Y_2\hfill & \hfill Y_3\hfill \end{array}\right]\left\{\begin{array}{c}\hfill q_1\hfill \\ \hfill q_2\hfill \\ \hfill q_3\hfill \end{array}\right\}$E17

Therefore, once the modal participation factors are solved, the transient response can be achieved using Eq. (16) by modal superposition.

In short, the conditions for the equivalence between the Zener and the exponential models have been established. Using these nonviscous damping models, a 1 dof system has been analyzed and three roots have been obtained. The analytical response has been reached without solving and solving the internal variable. For the former, Laplace transform properties have been employed. For the latter, using the initial value of the internal variable, modal superposition has been applied.

2.2. Continuous systems

Next, different procedures for solving the dynamics of an exponentially damped rod are described. First, the analytical solution is derived by means of modal superposition, solving and without solving internal variables. Then, two different finite element formulations are proposed; one is derived for direct integration methods, and the other to apply modal superposition procedures. The latter allows analyzing the impact of nonviscous modes into the global response.

Figure 2.

Representation of a continuous and homogeneous rod.

2.2.1. Analytical solution

Regarding the analytical solution, two different procedures are presented. The first developed procedure solves internal variables and the solution is obtained by modal superposition. In the second one, internal variables are avoided and a third-order time derivative field equation is found. The solution conditions are obtained and applied reducing the system equation to that reached using internal variables and enabling the application of modal superposition (Figure 2).

2.2.1.1. Solution using internal variables

To solve the displacement field u(x, t) of the rod, the material behavior law

$\sigma \left(t\right)+\frac{{\overline{c}}_1}{E_1}\frac{\mathrm{d}\sigma \left(t\right)}{\mathrm{d}t}=E_0\epsilon \left(t\right)+{\overline{c}}_1\left(1+\frac{E_0}{E_1}\right)\frac{\mathrm{d}\epsilon \left(t\right)}{\mathrm{d}t}$E18

has to be taken into account, where $\epsilon \left(x,t\right)=\frac{\partial u\left(x,t\right)}{\partial x}$ represents the strain field and where σ denotes the stress, ε the strain, ${\overline{c}}_1$ represents the damping coefficient, t is the current time, and $E_1$ and $E_0$ are stiffness parameters. The governing equations are obtained using an internal variable field y(x, t). Its gradient ${\epsilon }_y\left(x,t\right)=\frac{\partial y\left(x,t\right)}{\partial x}$ is related with its time derivative

and the strain rate $\stackrel{.}{\epsilon }\left(x,t\right)$ as

$E_1{\epsilon }_y\left(x,t\right)=c_1\left[\stackrel{.}{\epsilon }\left(x,t\right)-{\stackrel{.}{\epsilon }}_y\left(x,t\right)\right]$E19

Besides, the force equilibrium equation satisfies

$\sigma \left(x,t\right)-c_1\left[\stackrel{.}{\epsilon }\left(x,t\right)-{\stackrel{.}{\epsilon }}_y\left(x,t\right)\right]-E_0\epsilon \left(x,t\right)=0$E20

Thus, the field equation can be written in matrix form yielding

$\left[\begin{array}{cc}\hfill \rho S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \ddot{u}\left(x,t\right)\hfill \\ \hfill \ddot{y}\left(x,t\right)\hfill \end{array}\right\}-\left[\begin{array}{cc}\hfill c_{1}S\hfill & \hfill -c_{1}S\hfill \\ \hfill -c_{1}S\hfill & \hfill c_{1}S\hfill \end{array}\right]\frac{{\partial }^2}{\partial x^2}\left\{\begin{array}{c}\hfill \stackrel{.}{u}\left(x,t\right)\hfill \\ \hfill \stackrel{.}{y}\left(x,t\right)\hfill \end{array}\right\}-\left[\begin{array}{cc}\hfill E_{0}S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E_{1}S\hfill \end{array}\right]\frac{{\partial }^2}{\partial x^2}\left\{\begin{array}{c}\hfill u\left(x,t\right)\hfill \\ \hfill y\left(x,t\right)\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right\}$E21

where ρ is material density and S denotes cross-sectional area. It should be remarked the correspondence between Eqs. (19) and (21), to those relating the 1 dof case (Eqs. (1)–(3)). Thus, applying variable separation for u(x, t) and y(x, t) as

$\mathbf{z}\left(x,t\right)=\left\{\begin{array}{c}\hfill u\left(x,t\right)\hfill \\ \hfill y\left(x,t\right)\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill U\left(x\right)\hfill \\ \hfill Y\left(x\right)\hfill \end{array}\right\}{\mathrm{e}}^{st}$E22

U(x) being the spatial component, Y(x) the component for the internal variable y, and s a complex variable, the time dependence can be eliminated yielding

$s^2\left[\begin{array}{cc}\hfill \rho \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\left\{\begin{array}{c}\hfill U\left(x\right)\hfill \\ \hfill Y\left(x\right)\hfill \end{array}\right\}+\left[\begin{array}{cc}\hfill -s{c}_1-E_0\hfill & \hfill s{c}_1\hfill \\ \hfill -s{c}_1\hfill & \hfill s{c}_1+E_1\hfill \end{array}\right]\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\left\{\begin{array}{c}\hfill U\left(x\right)\hfill \\ \hfill Y\left(x\right)\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right\}$E23

The solution of the eigenproblem (23) provides the eigenfunctions ${\left\{\begin{array}{cc}\hfill U_j\left(x\right)\hfill & \hfill Y_j\left(x\right)\hfill \end{array}\right\}}^{\mathrm{T}}$ and the eigenvalues λ_j given by

$\left\{\begin{array}{c}\hfill U_j\left(x\right)\hfill \\ \hfill Y_j\left(x\right)\hfill \end{array}\right\}=q_j\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill A_j\hfill \end{array}\right\}sin\left({\lambda }_{j}x\right)$E24

and

${\lambda }_j=\frac{\mathrm{\pi }}{2\ell }\left(2j-1\right)$E25

respectively, j being the mode number, q_j the modal participation factor, A_j the amplitude relationship, and ℓ the rod length. As it can be verified through Eq. (25), all eigenvalues λ_j are real. Actually, proportional damping can be considered based on the homogeneity of the material. Therefore, the system presents normal modes, the eigenfunctions U_j(x) being those of the undamped one (see, e.g., Muller [3], for the solution of an undamped rod).

Using eigenfunctions ${\left\{\begin{array}{cc}\hfill U_j\left(x\right)\hfill & \hfill Y_j\left(x\right)\hfill \end{array}\right\}}^{\mathrm{T}}$ and the eigenvalues λ_j, the characteristic equation for Eq. (23) yields

$s^3+s^2\frac{E_1}{c_1}+s\frac{E_1+E_0}{\rho }{\lambda }_j^2+\frac{E_0{E}_1}{c_1\rho }{\lambda }_j^2=0$E26

which is a third-order equation, analogous to that of 1 dof case. Consequently, the solution is configured by infinite groups of three roots $s_{j,1}^{*}$, $s_{j,2}^{*}$, and $s_{j,3}$, associated to each eigenfunction ${\left\{\begin{array}{cc}\hfill U_j\left(x\right)\hfill & \hfill Y_j\left(x\right)\hfill \end{array}\right\}}^{\mathrm{T}}$ and eigenvalue λ_j. Thus, in each vibration mode, two elastic and one nonviscous components are involved.

Thus, the time response by means of modal superposition can be written as

$\mathbf{z}\left(x,t\right)={\displaystyle \sum _{j=1}^{\infty }\left(q_{j,1}{\mathbf{Z}}_{j,1}{\mathrm{e}}^{s_{j,1}^{\ast }t}+q_{j,2}{\mathbf{Z}}_{j,2}{\mathrm{e}}^{s_{j,2}^{\ast }t}+q_{j,3}{\mathbf{Z}}_{j,3}{\mathrm{e}}^{s_{j,3}t}\right)sin\left({\lambda }_{j}x\right)}$E27

where q_j,k ($k=1,2,3$) are the modal participation factors and Z_j,k are eigenvectors given by

${\mathbf{Z}}_{j,k}=\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill A_{j,k}\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill 1\hfill \\ \hfill \frac{1}{1+\frac{E_1}{s_{j,k}{c}_1}}\hfill \end{array}\right\}$E28

It should be noted that the resulting expression for the eigenvectors of the continuous system hold in Eq. (28) is analogous to the one corresponding to the 1 dof case (see Eq. (15)). Therefore, the time response can be obtained by modal superposition by solving, for each jth mode the modal participation factors q_j,k. However, from the equations provided by the initial conditions ${\mathbf{z}}_0\left(x\right)$ and ${\stackrel{.}{\mathbf{z}}}_0\left(x\right)$,

${\mathbf{z}}_0\left(x\right)={\displaystyle \sum _{j=1}^{\infty }\left(q_{j,1}{\mathbf{Z}}_{j,1}+q_{j,2}{\mathbf{Z}}_{j,2}+q_{j,3}{\mathbf{Z}}_{j,3}\right)sin\left({\lambda }_{j}x\right)}$E29

and

${\stackrel{.}{\mathbf{z}}}_0\left(x\right)={\displaystyle \sum _{j=1}^{\infty }\left(s_{j,1}^{*}{q}_{j,1}{\mathbf{Z}}_{j,1}+s_{j,2}^{*}{q}_{j,2}{\mathbf{Z}}_{j,2}+s_{j,3}{q}_{j,3}{\mathbf{Z}}_{j,3}\right)sin\left({\lambda }_{j}x\right)}$E30

these q_j,k cannot be directly solved because infinite unknowns are present. To overcome this situation, the consideration of proportional damping can be recalled, concluding that the eigenfunctions U_j(x) have to be orthogonal with respect to the mass and stiffness operators. Hence, the modes can be decoupled according to (see, e.g., Ref. [5] for details)

${\displaystyle {\int }_{\ell }{U}_i\left(x\right)\rho S{U}_j\left(x\right)\mathrm{d}x=0},\mathrm{f}\mathrm{o}\mathrm{r}i\ne j$E31

Accordingly, applying Eqs. (29)–(31), it is obtained

${\displaystyle {\int }_{\ell }{\mathbf{z}}_0\left(x\right)sin\left({\lambda }_{j}x\right)\mathrm{d}x=\left(q_{j,1}{\mathbf{Z}}_{j,1}+q_{j,2}{\mathbf{Z}}_{j,2}+q_{j,3}{\mathbf{Z}}_{j,3}\right){\displaystyle {\int }_{\ell }{sin}^2\left({\lambda }_{j}x\right)}\mathrm{d}x}$E32

and

${\displaystyle {\int }_{\ell }{\stackrel{.}{\mathbf{z}}}_0\left(x\right)sin\left({\lambda }_{j}x\right)\mathrm{d}x=\left(s_{j,1}^{\ast }{q}_{j,1}{\mathbf{Z}}_{j,1}+s_{j,2}^{\ast }{q}_{j,2}{\mathbf{Z}}_{j,2}+s_{j,3}{q}_{j,3}{\mathbf{Z}}_{j,3}\right){\displaystyle {\int }_{\ell }{sin}^2\left({\lambda }_{j}x\right)}\mathrm{d}x}$E33

which constitute a four-equation and three unknown system. Here, it should be pointed out (see Eq. (19)) that the internal variable y(x, t) and its time derivative $\dot{y}(x,t)$ are linearly combined together with the displacement u(x, t). Consequently, the rank of the equation system is 3. Therefore, only the initial conditions are needed, namely $u_0\left(x\right)$ and ${\dot{u}}_0(x)$, together with the initial value for the internal variable field $y_0\left(x\right)$ or together with its time derivative ${\dot{y}}_0(x)$. The initial condition for $y_0\left(x\right)$ and its time derivative ${\dot{y}}_0(x)$ have to satisfy $y_0\left(x\right)=0$ and ${\dot{y}}_0(x)={\dot{u}}_0(x)$, respectively. Hence, the equation to resolve each group $q_{j,1}$, $q_{j,2}$, and $q_{j,3}$ of modal participation factors is given by

$\frac{2}{\ell }{\displaystyle {\int }_{\ell }\left\{\begin{array}{c}\hfill u_0\left(x\right)\hfill \\ \hfill {\stackrel{.}{u}}_0\left(x\right)\hfill \\ \hfill 0\hfill \end{array}\right\}sin\left({\lambda }_{j}x\right)\mathrm{d}x}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill s_{j,1}^{\ast }\hfill & \hfill s_{j,2}^{\ast }\hfill & \hfill s_{j,3}\hfill \\ \hfill A_{j,1}\hfill & \hfill A_{j,2}\hfill & \hfill A_{j,3}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill q_{j,1}\hfill \\ \hfill q_{j,2}\hfill \\ \hfill q_{j,3}\hfill \end{array}\right\}$E34

Thus, once the modal participation factors are solved for each mode, the transient response by modal superposition can be achieved, giving

$u\left(x,t\right)={\displaystyle \sum _{j=1}^{\infty }\left(q_{j,1}{\mathrm{e}}^{s_{j,1}^{\ast }t}+q_{j,2}{\mathrm{e}}^{s_{j,2}^{\ast }t}+q_{j,3}{\mathrm{e}}^{s_{j,3}t}\right)}sin\left({\lambda }_{j}x\right)$E35

2.2.1.2. Solution without using internal variables

Next, the analytical solution for the displacement field u(x, t) of the rod is reached without introducing internal variables. Thus, recalling the material behavior law (Eq. (18)), the field equation can be written as

$\frac{c_1}{E_1}\rho S\stackrel{⃛}{u}\left(x,t\right)+\rho S\ddot{u}\left(x,t\right)-c_{1}S\left(1+\frac{E_0}{E_1}\right)\frac{{\partial }^2\stackrel{.}{u}\left(x,t\right)}{\partial x^2}-E_{0}S\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}=0$E36

Eq. (36) is a third-order equation in time, analogous to the one relating the 1 dof case (Eq. (3)). Applying variable separation for u(x, t) as u(x, t) = U(x) e^st it yields

$-\left[s^3\rho \frac{c_1}{E_1}+s^2\rho \right]U\left(x\right)+\left[s{c}_1\left(1+\frac{E_0}{E_1}\right)+E_0\right]\frac{{\partial }^{2}U\left(x\right)}{\partial x^2}=0$E37

The solution of the eigenproblem (37) provides the eigenfunctions U_j(x), given by U_j(x) = B_j sin(λ_jx) where the eigenvalues λ_j are those of Eq. (25), j being the mode number and B_j the amplitude. Therefore, the characteristic Eq. (26) is reached also.

Hence, if the solution is going to be derived without using internal variables, it is formed by infinite groups of three roots $s_{j,1}^{*}$, $s_{j,2}^{*}$, and $s_{j,3}$, associated to each eigenfunction U_j(x) and eigenvalue λ_j. Although, three initial conditions are needed to solve Eq. (36). These are the initial displacement $u\left(x,0\right)=u_0\left(x\right)$, the initial velocity

, and the initial acceleration $\ddot{u}\left(x,0\right)={\ddot{u}}_0\left(x\right)$, but, it can be proved that this initial acceleration ${\ddot{u}}_0\left(x\right)$ can be written as a function of the initial displacement $u_0\left(x\right)$.

In fact, using the Boltzmann superposition principle [6], the memory of a viscoelastic material can be properly modeled using hereditary models. The stress can be evaluated using relaxation functions R(t) through convolution integrals given by

$\sigma \left(t\right)={\displaystyle {\int }_0^{t}E\left(t-\tau \right)\stackrel{.}{\epsilon }\left(\tau \right)\mathrm{d}\tau }=E_0\epsilon \left(t\right)+{\displaystyle {\int }_0^{t}R\left(t-\tau \right)\stackrel{.}{\epsilon }\left(\tau \right)\mathrm{d}\tau }$E38

where E(t) is the relaxation modulus of the material, τ corresponds to the retardation time, and $\left(\stackrel{}{•}\right)$ represents time derivative. A relaxation function R(t) widely used in the literature [7–9] is the exponential model R(t) = cμ e^− μ t, where μ is the relaxation parameter and c represents the damping coefficient. Hence, using Eq. (38) and writing the force equilibrium for the continuous rod, it is obtained

$\rho \ddot{u}\left(x,t\right)=E_0\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+{\displaystyle {\int }_0^{t}R\left(t-\tau \right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\mathrm{d}\tau }$E39

Assuming causality, the material presents no memory prior to initial conditions and therefore, the initial acceleration ${\ddot{u}}_0\left(x\right)$ can be written as

${\ddot{u}}_0\left(x\right)=\frac{E_0}{\rho }\frac{{\partial }^{2}u\left(x,0\right)}{\partial x^2}$E40

Hence, recalling Eq. (37), from Eq. (40) the initial acceleration yields

${\ddot{u}}_0\left(x\right)=-\frac{E_0}{\rho }{\lambda }^2{u}_0\left(x\right)$E41

Thus, the time response can be written also as

$u\left(x,t\right)={\displaystyle \sum _{j=1}^{\infty }\left(B_{j,1}{\mathrm{e}}^{s_{j,1}^{\ast }t}+B_{j,2}{\mathrm{e}}^{s_{j,2}^{\ast }t}+B_{j,3}{\mathrm{e}}^{s_{j,3}t}\right)}sin\left({\lambda }_{j}x\right)$E42

where B_j,k ($k=1,2,3$) are the modal participation factors. Thus, these modal participation factors B_j,k must be solved for each jth mode in order to solve the time response by means of modal superposition. However, from the equations derived from the initial conditions $u_0\left(x\right)$, ${\dot{u}}_0(x)$, and ${\ddot{u}}_0\left(x\right)$

$u_0\left(x\right)={\displaystyle \sum _{j=1}^{\infty }\left(B_{j,1}+B_{j,2}+B_{j,3}\right)}sin\left({\lambda }_{j}x\right)$E43

${\stackrel{.}{u}}_0\left(x\right)={\displaystyle \sum _{j=1}^{\infty }\left(s_{j,1}^{*}{B}_{j,1}+s_{j,2}^{*}{B}_{j,2}+s_{j,3}{B}_{j,3}\right)}sin\left({\lambda }_{j}x\right)$E44

and

${\ddot{u}}_0\left(x\right)={\displaystyle \sum _{j=1}^{\infty }\left[{\left(s_{j,1}^{*}\right)}^2{B}_{j,1}+{\left(s_{j,2}^{*}\right)}^2{B}_{j,2}+{\left(s_{j,3}\right)}^2{B}_{j,3}\right]}sin\left({\lambda }_{j}x\right)$E45

these B_j,k cannot be directly solved because infinite unknowns are present. Hence, applying the orthogonality conditions, the modes can be decoupled according to Eq. (31) where variable separation together with Eq. (25) is needed. Consequently, it is obtained

${\displaystyle {\int }_{\ell }{u}_0\left(x\right)sin\left({\lambda }_{j}x\right)\mathrm{d}x=\left(B_{j,1}+B_{j,2}+B_{j,3}\right){\displaystyle {\int }_{\ell }{sin}^2\left({\lambda }_{j}x\right)}\mathrm{d}x}$E46

${\displaystyle {\int }_{\ell }{\stackrel{.}{u}}_0\left(x\right)sin\left({\lambda }_{j}x\right)\mathrm{d}x=\left(s_{j,1}^{*}{B}_{j,1}+s_{j,2}^{*}{B}_{j,2}+s_{j,3}{B}_{j,3}\right){\displaystyle {\int }_{\ell }{sin}^2\left({\lambda }_{j}x\right)}\mathrm{d}x}$E47

and

${\displaystyle {\int }_{\ell }{\ddot{u}}_0\left(x\right)sin\left({\lambda }_{j}x\right)\mathrm{d}x=\left[{\left(s_{j,1}^{*}\right)}^2{B}_{j,1}+{\left(s_{j,2}^{*}\right)}^2{B}_{j,2}+{\left(s_{j,3}\right)}^2{B}_{j,3}\right]{\displaystyle {\int }_{\ell }{sin}^2\left({\lambda }_{j}x\right)}\mathrm{d}x}$E48

which constitute a 3-equation and 3-unknown system. Hence, the equation to solve each group $B_{j,1}$, $B_{j,2}$, and $B_{j,3}$ of modal participation factors is given by

$\frac{2}{\ell }{\displaystyle {\int }_{\ell }\left\{\begin{array}{c}\hfill u_0\left(x\right)\hfill \\ \hfill {\stackrel{.}{u}}_0\left(x\right)\hfill \\ \hfill {\ddot{u}}_0\left(x\right)\hfill \end{array}\right\}sin\left({\lambda }_{j}x\right)\mathrm{d}x}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill s_{j,1}^{\ast }\hfill & \hfill s_{j,2}^{\ast }\hfill & \hfill s_{j,3}\hfill \\ \hfill {\left(s_{j,1}^{\ast }\right)}^2\hfill & \hfill {\left(s_{j,2}^{\ast }\right)}^2\hfill & \hfill {\left(s_{j,3}\right)}^2\hfill \end{array}\right]\left\{\begin{array}{c}\hfill B_{j,1}\hfill \\ \hfill B_{j,2}\hfill \\ \hfill B_{j,3}\hfill \end{array}\right\}$E49

or taking into account Eq. (41), it can be also written as

$\frac{2}{\ell }{\displaystyle {\int }_{\ell }\left\{\begin{array}{c}\hfill u_0\left(x\right)\hfill \\ \hfill {\stackrel{.}{u}}_0\left(x\right)\hfill \\ \hfill -\frac{E_0}{\rho }{\lambda }_j^2{u}_0\left(x\right)\hfill \end{array}\right\}sin\left({\lambda }_{j}x\right)\mathrm{d}x}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill s_{j,1}^{\ast }\hfill & \hfill s_{j,2}^{\ast }\hfill & \hfill s_{j,3}\hfill \\ \hfill {\left(s_{j,1}^{\ast }\right)}^2\hfill & \hfill {\left(s_{j,2}^{\ast }\right)}^2\hfill & \hfill {\left(s_{j,3}\right)}^2\hfill \end{array}\right]\left\{\begin{array}{c}\hfill B_{j,1}\hfill \\ \hfill B_{j,2}\hfill \\ \hfill B_{j,3}\hfill \end{array}\right\}$E50

It should be remarked that the system equations provided in Eq. (34) to solve the modal participation factors $q_{j,1}$, $q_{j,2}$, and $q_{j,3}$ are equivalent to that provided in Eq. (50) to solve those of the system without internal variables $B_{j,1}$, $B_{j,2}$, and $B_{j,3}$. In fact, it can be verified that the third row of Eq. (50) is a linear combination of the third and first rows of Eq. (34) as

${\mathrm{row}}_{3,39}={\alpha }_j\times {\mathrm{row}}_{1,22}+{\beta }_j\times {\mathrm{row}}_{3,22}$E51

where

and

.

As a conclusion it should be highlighted that if internal variables are avoided, the system equation reached (Eq. (49)) is equivalent to the one obtained using internal variables (Eq. (34)). Accordingly, the response of Eq. (34) can be derived by modal superposition through Eq. (35).

As a conclusion, the governing equations for a nonviscously damped rod have been obtained, using and without using internal variables. First, considering internal variables, the analytical response has been reached applying modal superposition. Then, without considering internal variables a third-order equation system in time is found where three initial conditions needed (see solve Eq. (36)). However, it has been proved that the initial acceleration ${\ddot{u}}_0\left(x\right)$ is a function of the initial displacement $u_0\left(x\right)$. Also, the equality between the equation systems (49) and (34) has been validated. Consequently, the response has been solved by modal superposition.

2.2.2. Finite element formulation

In this section, two finite element formulations are developed. These formulations allow to compute the motion (see Eq. (21) or (36)) of continuous rods with exponential damping. One is conceived to apply modal superposition and the other for direct integration methods.

The one for modal superposition allows to compute the time response (Eq. (21)) by solving two internal variables per finite element. The one created for direct integration, solves Eq. (36) by direct methods. This solution is accurate and efficient (no internal variables are solved).

2.2.2.1. Finite element formulation for modal superposition

Next, a finite element formulation to solve Eq. (21) is presented. It is a two-node formulation, leading to a four dof nodal displacement vector ${\mathbf{z}}_{\mathrm{e}}\left(t\right)={\left\{\begin{array}{cccc}\hfill u_1\left(t\right)\hfill & \hfill y_1\left(t\right)\hfill & \hfill u_2\left(t\right)\hfill & \hfill y_2\left(t\right)\hfill \end{array}\right\}}^{\mathrm{T}}$ considering the nodal displacements $u_1\left(t\right)$ and $u_2\left(t\right)$ and the corresponding internal variables $y_1\left(t\right)$ and $y_2\left(t\right)$. The approximate displacement z_e(x, t) can be written as

${\mathbf{z}}_{\mathrm{e}}\left(x,t\right)=\mathbf{N}\left(x\right){\mathbf{z}}_{\mathrm{e}}\left(t\right)=\left[\begin{array}{cccc}\hfill N_1\left(x\right)\hfill & \hfill 0\hfill & \hfill N_2\left(x\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill N_1\left(x\right)\hfill & \hfill 0\hfill & \hfill N_2\left(x\right)\hfill \end{array}\right]\left\{\begin{array}{c}\hfill u_1\left(t\right)\hfill \\ \hfill y_1\left(t\right)\hfill \\ \hfill u_2\left(t\right)\hfill \\ \hfill y_2\left(t\right)\hfill \end{array}\right\}$E52

where N(x) is the matrix of interpolation functions, in which linear interpolation ones are employed (see, e.g., Refs. [5, 10] for details about the finite element method). Thus, the weighted residual technique from the Galerkin point of view gives

$\begin{array}{l}{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\left[\begin{array}{cc}\hfill \rho S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]{\ddot{\mathbf{z}}}_{\mathrm{e}}\left(x,t\right)}\mathrm{d}x-{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\left[\begin{array}{cc}\hfill c_{1}S\hfill & \hfill -c_{1}S\hfill \\ \hfill -c_{1}S\hfill & \hfill c_{1}S\hfill \end{array}\right]\frac{{\partial }^2{\stackrel{.}{\mathbf{z}}}_{\mathrm{e}}\left(x,t\right)}{\partial x^2}}\mathrm{d}x-\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\left[\begin{array}{cc}\hfill E_{0}S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E_{1}S\hfill \end{array}\right]\frac{{\partial }^2{\mathbf{z}}_{\mathrm{e}}\left(x,t\right)}{\partial x^2}}\mathrm{d}x=\mathbf{0}\text{.}\end{array}$E53

Solving by parts the second and the third integrals, Eq. (53) yields

$\begin{array}{l}{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\left[\begin{array}{cc}\hfill \rho S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\mathbf{N}\left(x\right)\mathrm{d}x}{\ddot{\mathbf{z}}}_{\mathrm{e}}\left(t\right)+{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}\frac{\mathrm{d}{\mathbf{N}}^{\mathrm{T}}\left(x\right)}{\mathrm{d}x}\left[\begin{array}{cc}\hfill c_{1}S\hfill & \hfill -c_{1}S\hfill \\ \hfill -c_{1}S\hfill & \hfill c_{1}S\hfill \end{array}\right]\frac{\mathrm{d}\mathbf{N}\left(x\right)}{\mathrm{d}x}\mathrm{d}x}{\stackrel{.}{\mathbf{z}}}_{\mathrm{e}}\left(t\right)+\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}\frac{\mathrm{d}{\mathbf{N}}^{\mathrm{T}}\left(x\right)}{\mathrm{d}x}\left[\begin{array}{cc}\hfill E_{0}S\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E_{1}S\hfill \end{array}\right]\frac{\mathrm{d}\mathbf{N}\left(x\right)}{\mathrm{d}x}\mathrm{d}x}{\mathbf{z}}_{\mathrm{e}}\left(t\right)={\mathbf{F}}_{\mathrm{e}}\left(t\right)\text{.}\end{array}$E54

Therefore, a classical second-order equation with constant coefficients

${\mathbf{M}}_{\mathrm{e}}{\ddot{\mathbf{z}}}_{\mathrm{e}}\left(t\right)+{\mathbf{C}}_{\mathrm{e}}{\stackrel{.}{\mathbf{z}}}_{\mathrm{e}}\left(t\right)+{\mathbf{K}}_{\mathrm{e}}{\stackrel{.}{\mathbf{z}}}_{\mathrm{e}}\left(t\right)={\mathbf{F}}_{\mathrm{e}}\left(t\right)$E55

is reached. Consequently, the mass, damping and stiffness elementary matrices yield

${\mathbf{M}}_{\mathrm{e}}=\frac{\rho S{\ell }_{\mathrm{e}}}{6}\left[\begin{array}{cccc}\hfill 2\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$E56

${\mathbf{C}}_{\mathrm{e}}=\frac{c_{1}S}{{\ell }_{\mathrm{e}}}\left[\begin{array}{cccc}\hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \end{array}\right]$E57

and

${\mathbf{K}}_{\mathrm{e}}=\frac{S}{{\ell }_{\mathrm{e}}}\left[\begin{array}{cccc}\hfill E_0\hfill & \hfill 0\hfill & \hfill -E_0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E_1\hfill & \hfill 0\hfill & \hfill -E_1\hfill \\ \hfill -E_0\hfill & \hfill 0\hfill & \hfill E_0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -E_1\hfill & \hfill 0\hfill & \hfill E_1\hfill \end{array}\right]$E58

Assembling the elementary matrices, the state-space equation

$\left[\begin{array}{cc}\hfill \mathbf{C}\hfill & \hfill \mathbf{M}\hfill \\ \hfill \mathbf{M}\hfill & \hfill \mathbf{0}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \stackrel{.}{\mathbf{z}}\left(t\right)\hfill \\ \hfill \ddot{\mathbf{z}}\left(t\right)\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill \mathbf{F}\left(t\right)\hfill \\ \hfill \mathbf{0}\hfill \end{array}\right\}+\left[\begin{array}{cc}\hfill -\mathbf{K}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{M}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \mathbf{z}\left(t\right)\hfill \\ \hfill \stackrel{.}{\mathbf{z}}\left(t\right)\hfill \end{array}\right\}$E59

is reached. It should be remarked that the size of the eigenproblem obtained from Eq. (59) becomes $4N$. However, reminding that the extra degrees of freedom have no associated mass, the rank of the matrix system is $3N$. Hence, the response yields

$\mathbf{z}\left(t\right)={\displaystyle \sum _{j=1}^{3N}{q}_j{\mathbf{Z}}_j{\mathrm{e}}^{s_{j}t}}$E60

where q_j represents the modal participation factors, Z_j the system eigenvectors, and s_j the system roots. To compute the $3N$ modal participation factors, $2N$ equations can be written from the initial conditions ${\mathbf{u}}_0=\mathbf{u}\left(0\right)$ and

, where the N remaining equations can be derived from the initial values for the internal variables, satisfying

. Therefore, once the modal participation factors q_j are solved, the response can be computed by modal superposition through Eq. (60).

2.2.2.2. Finite element formulation for direct integration

Next, a finite element formulation to solve Eq. (36) is presented. For that, the rod is discretized in two-node finite elements of length ℓ_e with linear interpolation functions. The weighted residual technique from the Galerkin point of view gives

$\begin{array}{l}{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\frac{c_1}{E_1}\rho S\mathbf{N}\left(x\right)\mathrm{d}x}{\stackrel{⃛}{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\rho S\mathbf{N}\left(x\right)\mathrm{d}x}{\ddot{\mathbf{u}}}_{\mathrm{e}}\left(t\right)-\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)c_{1}S\left(1+\frac{E_0}{E_1}\right)\frac{{\partial }^2{\stackrel{.}{\mathbf{u}}}_{\mathrm{e}}\left(x,t\right)}{\partial x^2}\mathrm{d}x}-\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)E_{0}S\frac{{\partial }^2{\mathbf{u}}_{\mathrm{e}}\left(x,t\right)}{\partial x^2}\mathrm{d}x}=\mathbf{0}\text{,}\end{array}$E61

where

denotes the transposition operator, N(x) is the matrix of interpolation functions, and ${\mathbf{u}}_{\mathrm{e}}\left(t\right)={\left\{\begin{array}{cc}\hfill u_1\left(t\right)\hfill & \hfill u_2\left(t\right)\hfill \end{array}\right\}}^{\mathrm{T}}$ is the nodal displacement vector. Solving by parts the third and the fourth integrals of Eq. (61), it gives

$\begin{array}{l}{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\frac{c_1}{E_1}\rho S\mathbf{N}\left(x\right)\mathrm{d}x}{\stackrel{⃛}{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+{\displaystyle {\int }_{{\ell }_{\mathrm{e}}}{\mathbf{N}}^{\mathrm{T}}\left(x\right)\rho S\mathbf{N}\left(x\right)\mathrm{d}x}{\ddot{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}\frac{\mathrm{d}{\mathbf{N}}^{\mathrm{T}}\left(x\right)}{\mathrm{d}x}{c}_{1}S\left(1+\frac{E_0}{E_1}\right)\frac{\mathrm{d}\mathbf{N}\left(x\right)}{\mathrm{d}x}\mathrm{d}x}{\stackrel{.}{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+\\ {\displaystyle {\int }_{{\ell }_{\mathrm{e}}}\frac{\mathrm{d}{\mathbf{N}}^{\mathrm{T}}\left(x\right)}{\mathrm{d}x}{E}_{0}S\frac{\mathrm{d}\mathbf{N}\left(x\right)}{\mathrm{d}x}\mathrm{d}x}{\mathbf{u}}_{\mathrm{e}}\left(t\right)={\mathbf{F}}_{\mathrm{e}}\left(t\right)\text{,}\end{array}$E62

where F_e(t) represents the external forces nodal vector. Thus, the third order in time matrix system

${\mathbf{J}}_{\mathrm{e}}{\stackrel{⃛}{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+{\mathbf{M}}_{\mathrm{e}}{\ddot{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+{\mathbf{C}}_{\mathrm{e}}{\stackrel{.}{\mathbf{u}}}_{\mathrm{e}}\left(t\right)+{\mathbf{K}}_{\mathrm{e}}{\mathbf{u}}_{\mathrm{e}}\left(t\right)={\mathbf{F}}_{\mathrm{e}}\left(t\right)$E63

is reached, in which the elementary matrices to be assembled are

${\mathbf{J}}_{\mathrm{e}}=\frac{\rho S{\ell }_{\mathrm{e}}}{6}\frac{c_1}{E_1}\left[\begin{array}{cc}\hfill 2\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 2\hfill \end{array}\right]$E64

${\mathbf{M}}_{\mathrm{e}}=\frac{\rho S{\ell }_{\mathrm{e}}}{6}\left[\begin{array}{cc}\hfill 2\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 2\hfill \end{array}\right]$E65

${\mathbf{C}}_{\mathrm{e}}=\frac{c_{1}S}{{\ell }_{\mathrm{e}}}\left(1+\frac{E_0}{E_1}\right)\left[\begin{array}{cc}\hfill 1\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 1\hfill \end{array}\right]$E66

and

${\mathbf{K}}_{\mathrm{e}}=\frac{E_{0}S}{{\ell }_{\mathrm{e}}}\left[\begin{array}{cc}\hfill 1\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 1\hfill \end{array}\right]$E67

Assembling all the finite element matrices, Eq. (63) yields

$\mathbf{J}\stackrel{⃛}{\mathbf{u}}\left(t\right)+\mathbf{M}\ddot{\mathbf{u}}\left(t\right)+\mathbf{C}\stackrel{.}{\mathbf{u}}\left(t\right)+\mathbf{K}\mathbf{u}\left(t\right)=\mathbf{F}\left(t\right)$E68

To solve Eq. (68) by means of traditional direct methods of structural dynamics, the procedure employed in the Cortes Mateos Elejabarrieta (CME) method [11] is used next. For that, the backward definition of the first derivative $\stackrel{.}{\mathbf{u}}\left(t\right)$ of the displacement u(t) at the time $t_{n+1}$ can be approximated as

, where $\Delta t=t_{n+1}-t_n$ is a finite time step. Hence, Eq. (68) yields

${\mathbf{M}}_{\mathrm{eq}}{\ddot{\mathbf{u}}}_{n+1}+\mathbf{C}{\stackrel{.}{\mathbf{u}}}_{n+1}+\mathbf{K}{\mathbf{u}}_{n+1}={{\mathbf{F}}_{\mathrm{eq},}}_{n+1}$E69

with ${\mathbf{M}}_{\mathrm{eq}}=\mathbf{M}+\Delta t^{-1}\mathbf{J}$ and F_eq,n + 1 = F_n + 1 + Δt^− 1Jü_n representing an equivalent second-order forced system. The system response can be computed through direct integration methods as, for example, the Newmark method.

However, assuming causality, there is no memory prior the initial conditions are applied and therefore, the initial acceleration ü₀ = ü(0) is computed from

$\mathbf{M}{\ddot{\mathbf{u}}}_0+\mathbf{K}{\mathbf{u}}_0={\mathbf{F}}_0$E70

3. Experimental analysis of the vibrational response of an adhesively bonded beam

In this section, an experimental procedure for studying the influence that geometrical properties of adhesive joints have on the vibrational response of a metallic beam doubly supported on viscoelastic adhesive joints is presented. A test bench has been specifically constructed for the experimental program. At a first step, the modal shapes have been experimentally identified. Regarding the experiments, the beams are seismically excited and the influence of joint thickness and overlapping length on the beam motion is studied by computing the root-mean-square (rms) value of 21 transmissibility functions obtained along the length of the beam. The analysis is carried out on resonance frequencies, peak amplitudes, and modal loss factors.

3.1. Experimental program

The objective of the experiments is to obtain sets of frequency f dependent transmissibility functions T(f) to study the influence that overlapping length ${\ell }_0$ and joint thickness h have on the beam vibrational response.

Item	Use
PC 1	Controls the dSPACE card
dSPACE (+charge amplifier)	Activates the piezoelectric actuator
Piezoelectric actuator PPA40M	Imposes the motion to the adhesive base
Tri-axial accelerometer ICP PCB 356A16	Measures the adhesive base acceleration
PULSE B&K acquisition system	Samples and processes the transducer time signal
PC 2	Controls the PULSE B&K acquisition system
Laser interferometer (+signal amplifier)	Measures the beam velocity response

Table 1.

Measurement equipment.

3.1.1. Equipment

Figure 1 shows a schematic diagram of the experimental setup, whereas Table 1 gives the details about measuring equipment. Figure 1 shows three groups of components: excitation of the system, data acquisition and processing, and the adhesively bonded beam itself (Figure 3).

The base motion is imposed by a PPA40M piezoelectric actuator [12] controlled by a DS1104 dSPACE real-time control card [13]. This card governs the piezoelectric actuator by a control program developed in Simulink® that generates a digital control signal. This is converted into an analogical one by the dSPACE card. Finally, it is amplified and sent to the piezoelectric actuator.

Figure 3.

Experimental setup: (a) computer 1, (b) real-time dSPACE card, (c) signal conditioner, (d) piezoelectric device, (e) triaxial accelerometer (1), (f) beam specimen, (g) triaxial accelerometer (2), (h) PULSE acquisition system, (i) computer 2, (j) laser signal amplifier, (k) laser interferometer OFV-505.

Figure 4.

Test bench: general view.

The data acquisition is done by a B&K PULSE acquisition system [14]. Two response sensors are used: a triaxial accelerometer ICP 356A16 of PCB electronics [15] for measuring the adhesive base motion and a laser interferometer OFV-505 of POLYTEC [16] for measuring the beam response. Hence, to validate the excitation, the acceleration of the adhesive base is measured in three directions, whereas the beam motion is measured just transversally. This velocity signal

is derived to obtain the corresponding acceleration

.

The designed test bench ensures a repetitive test procedure. Mainly, it is composed by a moving and a steady support assembled on a rigid base. Figure 4 illustrates a general view of the test bench while Figure 5 shows the moving support.

Figure 5.

Test bench: detailed view of the excitation system.

Figure 4 shows the test bench and Figure 5 shows the excitation system. Particularly, Figure 4 shows (i) the test bench base, (ii) where the excitation support, (iii) the steady support, and (iv) an adhesively bonded beam can be appreciated. Figure 5 shows the components of the excitation support; these are: (i) the actuator base, (ii) the piezoelectric device used as excitatory (PPA40M of CEDRAT technologies), (iii) the base of the adhesive, and (iv) the car-guide set. It should be noted that the seismic motion is directly imposed to the adhesive base by the piezoelectric actuator.

Two groups of metallic pieces are used to ensure the overlapping length ${\ell }_0$. Figure 6 shows these pieces. It should be noted that the evenly spaced slots are used to ensure five different overlapping lengths. The procedure is analogous for both sides and therefore left and right overlapping lengths are always nominally identical. Experimental results are presented for three different overlapping lengths, ${\ell }_0=10$, ${\ell }_0=30$, and ${\ell }_0=50\mathrm{mm}$, respectively.

Figure 6.

Components for controlling the overlapping length ${\ell }_0$.

Figure 7.

Components for controlling the thickness h.

Figure 8.

Detailed view of the hollow for the adhesive joint.

The adhesive joint thickness h is determined by a set of pieces that are mounted onto the ones used for determining the overlapping length. Analogously, it is repeated exactly and simultaneously in both supports ensuring that both adhesive joints are nominally identical. These devices support the beam during the adhesive curing. Figure 7 shows one of them.

Figure 7 shows placement strategy where five sets have been manufactured. Each set has a particular height related to a particular joint thickness.

The components determining overlapping length and joint thickness are assembled together onto the adhesive bases before the adhesive material is applied. Once they are mounted, the hollow for the adhesive joint is formed among them and the beam itself. This can be appreciated in Figure 8.

These components are removed prior to the measurements.

3.1.2. Experiments and data processing

The experiment consists on exciting the adhesive base and measuring the beam vibration response. The response is sequentially measured at 21 beam locations using the laser vibrometer.

Figure 9.

Imposed displacement to the adhesive base: time domain response $w_1\left(t\right)$.

Figure 10.

Imposed displacement to the adhesive base: signal-to-noise ratio for autospectrum $W_{11}\left(f\right)$.

The piezoelectric actuator receives a 0.8 V step signal that induces the displacement of the adhesive base shown in Figure 9. The signal-to-noise ratio for the autospectrum of $W_{11}\left(f\right)$ represented in Figure 10. This signal-to-noise ratio is calculated as the ratio for the signal at the adhesive base while the beam is vibrating, to the signal while the beam is at rest. The adhesive base displacement was computed by numerical integration from the acceleration measured signal. The measure was carried out for a sample period $T_s=2\mathrm{s}$ , the time resolution being $\Delta t=0.0025\mathrm{s}$, involving 800 data lines.

From the experimental data, transmissibility functions T_j(f) between adhesive base and each jth point beam are obtained by means of the $H_1\left(f\right)$ definition [14]

$H_1\left(f\right)=\frac{W_{12}\left(f\right)}{W_{11}\left(f\right)}$E71

where $W_{12}\left(f\right)$ represents the cross-spectrum between input $w_1\left(t\right)$ and output $w_2\left(t\right)$ and $W_{11}\left(f\right)$ the input autospectrum. The motion $w_2\left(t\right)$ is computed by numerical integration from the velocity

measured by means of the laser vibrometer [16] for each of the 21 points of the beam. The study is performed in the 0–200 Hz frequency range, the frequency resolution being $\Delta f=0.5\mathrm{Hz}$.

Prior to the experiments and aimed at understanding the flexural behavior of the vibrating beam, the first three flexural modal shapes are obtained from these transmissibility functions [17]. This is done just for one of the specimens.

For the experimental analysis, the rms value (for the set of 21 points) is computed as

$\mathrm{r}\mathrm{m}\mathrm{s}\left(T\left(f\right)\right)=\sqrt{\frac{1}{N}{\displaystyle \sum _{j=1}^N\left({\left|T_j\left(f\right)\right|}^2\right)}}$E72

with $N=21$. From the rms(T(f)), the ith resonance frequency f_i, peak amplitude A_i, and modal loss factor η_i of each specimen are studied.

3.1.3. Materials and specimens

The considered adhesive is a flexible one, concretely a Bostik™ modified silane commercially named ISR 70-03. Its mechanical behavior was analyzed in earlier works [18, 19], where the relaxation modulus and the complex modulus were determined. It was done under tensile strain where the strain level imposed was of 0.5%. The test specimens were obtained from plates of the cured adhesive that were manufactured using casts of 50 mm × 70 mm × h, where h represents the nominal thickness. Three casts were manufactured with different thickness h of 0.5, 1, and 1.5 mm using Teflon™, which ensures that after the curing of the adhesive, a plate of solid material can be easily demolded without degradation. For the present case, the curing time was 48 h for all plates; the curing process was carried out at room temperature and no specific equipment was employed [20].

Figure 11 represents the complex modulus E*(f) for the adhesive in the form of storage modulus E′(f) and loss factor η(f).

Figure 11.

Tested adhesive material complex modulus E*(f) in the form of storage modulus E′(f) and loss factor η(f).

The metallic beam specimens were manufactured from the same stainless steel sheet by water-jet cutting in order to ensure homogeneous properties and undeformed specimens. The nominal beam dimensions are 400 mm × 10 mm × 1 mm , representing length, width, and thickness, respectively; the experimental values being 397 mm × 10 mm × 1 mm.

To prepare the adhesively bonded beam specimens an adhesive curing time of 72 h was established. This curing process was carried out at normal conditions (room temperature and atmospheric pressure). The test sample preparation procedure is outlined as follows:

• At a first step, the devices determining overlapping length and joint thickness h (see Figures 6 and 7) are mounted onto the supports (see Figure 8).

• Then, the adhesive is applied. It should be emphasized that this is done in a single motion ensuring that the nozzle does not get in touch with the adhesion surface. Otherwise, void creation is promoted.

• Immediately, the beam is placed forcing the adhesive to fill the gap among the components and the beam itself. The spare adhesive material is removed once it is cured.

Five different configurations were tested with three specimens for each configuration. Table 2 shows, for each configuration, the nominal joint overlapping length and thickness together with the particular dimensions for each specimen. Also, the mean value is reported (in parenthesis). Three different overlapping lengths (BS1, BS2, and BS3) and three joint thicknesses (BS1, BS4, and BS5) were tested using the ISR 70–03 adhesive. It should be highlighted that constant length beam specimens are employed. This means that the higher the overlapping length, the lower the vibrating length.

Specimen group	Overlapping length ℓo(mm)			Joint thickness h (mm)
	Nominal	Excitation side individual (mean)	Steady side individual (mean)	Nominal	Excitation side individual (mean)	Steady side individual (mean)
BS1	55	55.2/55.0/55.0 (55.1)	55.1/55.1/55.2 (55.1)	0.5	0.6/0.6/0.5 (0.6)	0.5/0.5/0.6 (0.5)
BS2	35	35.2/35.1/35.2 (35.2)	35.1/35.2/35.0 (35.1)	0.5	0.6/0.5/0.6 (0.6)	0.6/0.5/0.5 (0.5)
BS3	15	15.2/15.2/15.1 (15.2)	15.1/15.1/15.2 (15.1)	0.5	0.5/0.5/0.6 (0.5)	0.6/0.6/0.6 (0.6)
BS4	55	55.1/55.1/55.2 (55.1)	55.0/55.3/55.1 (55.2)	1	1.2/1.0/1.2 (1.1)	1.1/1.1/1.0 (1.1)
BS5	55	55.1/55.1/55.0 (55.1)	55.1/55.2/55.1 (55.1)	2.5	2.6/2.5/2.6 (2.6)	2.5/2.5/2.6 (2.6)

Table 2.

Joint dimensions for the beam specimens: individual and mean values.

Figure 12.

Modal shapes below 200 Hz: rigid body mode.

Figure 13.

Modal shapes below 200 Hz: first bending mode.

Figure 14.

Modal shapes below 200 Hz: second bending mode.

Figure 15.

Modal shapes below 200 Hz: third bending mode.

3.2. Modal shapes and analysis

To understand the dynamic behavior of the adhesively bonded beams, the modal shapes are experimentally obtained from the transmissibility functions [17]. Figure 12 represents the rigid body mode and Figures 13–15 show the first, second, and third bending modes.

The analyzed system can be modeled as a pinned-pinned beam with torsion spring of stiffness C applied at each pinned joint (see Figure 16). These torsion springs represent the resistance to turn introduced by the joint to the beam.

Figure 16.

Model with pinned-pinned beam with torsion spring at the pinned joints.

For the analysis, only the bending modes are taken into account.

3.3. Results and discussion

Next, the influence of overlapping length and joint thickness h are presented and discussed. It should be reminded that the results shown represent the rms value of a set of 21 transmissibility functions where three different specimens are used for each configuration. The particular and average result is reported together with the standard deviation.

Table 3 collects the resonance frequencies, peak amplitudes, and modal loss factor of the each specimen tested. The loss factor η_i is computed by the Nyquist-circle method.

For the subsequently study, a distinction must be made between system stiffness and joint stiffness. System stiffness refers to the stiffness property that globally affects the mechanical behavior of the system, which is normally traduced in natural frequencies. The joint stiffness (corresponding to the previously mentioned torsion stiffness C represented in Figure 16) refers to a local stiffness property of the joint. Actually, it has not a significant influence on the global response (slightly affects natural frequencies), but affects to the vibration transmission mechanism.

Table 3.

Frequency domain results for the analyzed beam specimens.

Individual and mean values: (a) BS1; (b) BS2; (c) BS3; (d) BS4; (e) BS5; (f) BS6.

3.3.1. Overlapping length influence

Three overlapping lengths are analyzed: 55, 35, and 15 mm, the specimen groups being BS1, BS2, and BS3 (see Table 2), respectively, the mean thickness being around 0.5 mm.

Figure 17.

Transmissibility function amplitude |T(f)| for three different overlapping lengths ℓ_o.

Figure 17 shows the rms(T(f)) function for the set of the 21 transmissibility functions obtained.

From Figure 17, it should be pointed out that there are two resonances in the analyzed range except for, in which three can be appreciated.

• About resonance frequency. As expected, the bigger the overlapping length, the bigger the resonance frequencies. This is due to as the overlapping length increases, so does the system stiffness. This is based on the decrement of the vibrating length.

• About resonance peak amplitude. Analogously, the higher is the overlapping length, the higher are the resonance peak amplitudes. Besides, the amplitudes for the first mode are higher than those of the second mode.

• About modal loss factor. Involving modal loss factor (see Table 3), it should be remarked that the second modes present lower values than those of the first ones. It can be stated that the higher the overlapping length, the lower the modal loss factor. However, this difference is minimal.

3.3.2. Thickness influence

Three adhesive joint thickness values are analyzed. These are 0.5, 1.0, and 2.5 mm, the specimen groups being BS1, BS4, and BS5 (see Table 2), respectively. Small adhesive joint thickness is considered in order not to dramatically decrease the joint strength. Figure 18 illustrates the function for the corresponding sets.

Figure 18.

Transmissibility function amplitude |T(f)| for three joint thicknesses h.

• About resonance frequency. It can be concluded that the thinner the thickness, the higher the joint stiffness, and therefore, the resonance frequencies increase. But, according to the previously mentioned distinction between system and joint stiffness, by decrementing five times the joint thickness, natural frequencies increase only about 10%.

• About resonance peak amplitude. The higher amplitude is found for the first resonance of each case (see Table 3). For each mode, the amplitude presents an inverse evolution with thickness: the minimum amplitude is given for the highest thickness h₁ = 2.5 mm. It should be noted that decreasing thickness results in increasing joint stiffness (leading to amplitude decrement) and decrementing damping material amount (leading to amplitude augmentation). Nevertheless, the effect on damping results more significant concluding that the thinner the thickness the higher the amplitude.

• About modal loss factor. Within each mode, the lower loss factor corresponds to the thinnest case ($h_3=0.5\mathrm{mm}$) due to the joint stiffness being higher, the damping capacity decreases by material deformation. However, the modal loss factor presents a minimum for the intermediate thickness.

From this analysis, the following discussion can be made:

About stiffness

• The overlapping length modifies the system stiffness, whereas the joint thickness modifies the joint stiffness.

• To decrease the system stiffness implies a decrement in the resonance frequencies and peak amplitudes A_i.

• To decrease the joint stiffness decreases the amplitudes A_i, because the vibration is worse transmitted. However, the resonance frequencies remain practically unaffected.

About damping

• The overlapping length and joint thickness h modify damping capacity because they are related to the amount of viscoelastic material used.

• To decrease overlapping length, on the one hand, decreases joint damping capacity because of the amount of material is decreased. However, on the other hand, as the joint stiffness is decreased, the viscoelastic material deformation is augmented and so is the damping capacity.

• To decrease joint thickness decreases joint damping capacity because of the joint stiffness is increased, and therefore the deformations are dismissed. Besides, damping material amount is decreased.

As a general conclusion, it can be drawn that increasing overlapping length or joint thickness can lead to opposite effects. Hence, to know and optimize the behavior of these kinds of adhesive joints, exhaustive experiments would be needed, or failing that, numerical simulation should be required in order to separate and investigate the influence of each parameter.