# Dynamic responses of flexible-link mechanisms with passive/active damping treatment

Jean-François Deü, Ana Cristina Galucio, Roger Ohayon

► To cite this version:

Jean-François Deü, Ana Cristina Galucio, Roger Ohayon. Dynamic responses of flexible-link mechanisms with passive/active damping treatment. Computers & Structures, 2008, 86 (3-5), pp.258-265. 10.1016/j.compstruc.2007.01.028. hal-01519381

HAL Id: hal-01519381

https://hal.science/hal-01519381

Submitted on 7 May 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Dynamic responses of flexible-link mechanisms with
passive/active damping treatment

J.-F. Deü, A.C. Galucio, R. Ohayon

Conservatoire National des Arts et Métiers (CNAM), Structural Mechanics and Coupled Systems Laboratory,
Case 353, 292 rue Saint Martin, 75141 Paris Cedex 03, France

This work presents a finite element formulation for non-linear transient dynamic analysis of adaptive beams. The main contribution of this work concerns the development of an original co-rotational sandwich beam element, which allows large displacements and rotations, and takes active/passive damping into account. This element is composed of a viscoelastic core and elastic/piezoelectric laminated faces. The latter are modeled using classical laminate theory, where the electromechanical coupling is considered by modifying the stiffness of the piezoelectric layers. For the core, a four-parameter fractional derivative model is used to characterize its viscoelastic dissipative behavior. Equations of motion are solved using an incremental-iterative method based on the Newmark direct time integration scheme in conjunction with the Gru¨nwald approximation of fractional derivatives, and the Newton-Raphson algorithm.

Keywords: Piezoelectric material; Viscoelasticity; Fractional derivatives; Sandwich beam; Co-rotational formulation; Non-linear dynamics

# 1. Introduction

In the context of linear dynamics, the active constrained layer damping treatment is widely used due to its beneficial performance in attenuating structural vibrations. The investigation of such a treatment in the geometric non-linear case, more specifically for flexible multibody systems, is very sparsely analyzed. This paper presents a non-linear finite element formulation for dynamic analysis of adaptive sandwich beams composed of a viscoelastic core constrained by elastic/piezoelectric laminated faces. It is important to mention that flexible beams undergoing finite rotations have been investigated by several approaches namely the floating frame approach, the inertial frame approach or the co-rotational approach. Most of these methods are described in standard textbooks and the inter-

ested reader can be referred to the review article [10] which compares these methods and summarizes the recent developments in computational modeling of flexible multibody systems. In this paper, following the works of Crisfield [2], a co-rotational approach is used to describe the large displacement analysis of the sandwich beam. One of the main difficulties in modeling such structures is the characterization of the damping properties of the viscoelastic layer. In this study, a four-parameter fractional derivative model [1] is used to describe the frequency-dependence of the viscoelastic material. The advantage of employing fractional operators is related to the small number of parameters required to represent the material damping over a broad range of frequencies. Concerning the active component, the electromechanical coupling is taken into account by modifying the stiffness matrix of the piezoelectric layers (sensor) and by applying a mechanical load written in terms of the applied tension (actuator). Consequently, in the final finite element formulation, no electrical degree-of-freedom explicitly appears [5]. Elastic/piezoelectric faces are modeled using classical laminate theory whereas a standard three-layer sandwich theory is used for the face/core/
face system [9]. In order to implement the fractional derivative constitutive equations, the Grünwald formalism is used in conjunction with a time discretization scheme based on Newmark method. One of the particularities of the proposed algorithm lies on the storage of displacement history only. This considerably reduces the numerical efforts related to the non-locality of fractional operators. Finally, numerical examples, including a flexible robot arm and a slider-crank mechanism, are considered in order to show the effectiveness of the present co-rotational formulation.

## 2. Active/passive damping treatment

Consider a sandwich beam composed of a viscoelastic core and elastic/piezoelectric laminated faces. Hypotheses proposed in this section are valid in a local coordinate system, which continuously rotates and translates with the beam element. Furthermore, a linear strain definition is used. The first part is addressed to the kinematical (for the mechanical displacement) and electrostatic (for the electric potential) assumptions. In the second part, constitutive equations for the piezoelectric and viscoelastic materials are outlined. Finally, the internal virtual energy of the sandwich beam is derived.

### 2.1. Planar sandwich beam

The sandwich beam is modeled using Euler-Bernoulli assumptions for the faces and Timoshenko ones for the core. The elastic/piezoelectric faces are modeled using classical laminated theory with linear electric potential through the thickness of each piezoelectric sublayer. The mechanical displacement field within the ith layer can be written as

$$u_{xi}(x, z, t) = u_i(x, t) - (z - z_i)\theta_i(x, t) \quad (1a)$$

$$u_{zi}(x, z, t) = w(x, t) \quad (1b)$$

where the subscript i = a, b, c stands for upper (composed of $N_a$ sublayers), lower (composed of $N_b$ sublayers) and middle layers, respectively; $u_{xi}$ and $u_{zi}$ are the axial and transverse displacements of each layer; $u_i$ and $\theta_i$ are the axial displacement of the center line and the fiber rotation of each layer; and w is the transverse displacement (Fig. 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Sandwich beam kinematics.</figcaption></figure>

Let us introduce the mean and relative axial displacements given by $\bar{u} = (u_a + u_b)/2$ and $\tilde{u} = u_a - u_b$. Euler-Bernoulli hypotheses for the faces lead to $\theta_k = w'$ for $k = a, b$ with $(\cdot)' = \partial(\cdot)/\partial x$. As all layers are supposed to be perfectly bonded, the displacement continuity conditions at interface layers can be written as $u_{xa} = u_{xc}$ at $z = h_c/2$ and $u_{xb} = u_{xc}$ at $z = -h_c/2$. Therefore, axial displacements of the centerlines and rotations of each layer can be written in terms of $w'$ and the above defined variables $\bar{u}$ and $\tilde{u}$ as

$$\begin{align}
u_a &= \bar{u} + \frac{\tilde{u}}{2}, & u_b &= \bar{u} - \frac{\tilde{u}}{2}, & u_c &= \bar{u} + \frac{\tilde{h}}{4}w', \nonumber \\
\theta_a &= \theta_b = w', & \theta_c &= -\frac{\tilde{u} + \tilde{h}w'}{h_c} \tag{2}
\end{align}$$

$$\text{where } \tilde{h} = (h_a + h_b)/2 \text{ and } \tilde{h} = h_a - h_b.$$

From (1) and (2), and taking the hypothesis of plane stress state into account, we can write the axial strain of the ith layer $\epsilon_{1i}$ and the shear strain of the core $\epsilon_{5c}$ as follows:

$$\epsilon_{1i} = \epsilon_i + (z - z_i)\kappa_i, \quad \epsilon_{5c} = \gamma_c \quad (3)$$

where $\epsilon_i$ and $\kappa_i$ are the membrane strain and curvature of the ith layer, and $\gamma_c$ the shear strain of the core.

Concerning the electrostatic aspects, two main assumptions are taken under consideration. The first one concerns the electric potential, which is supposed to be linear within the thickness of each piezoelectric sublayer

$$\phi_{kj}(x, z, t) = \bar{\phi}_{kj}(x, t) + (z - z_{kj}) \frac{V_{kj}(x, t)}{h_{kj}} \quad (4)$$

where $\bar{\phi}_{kj} = (\phi_{kj}^{+} + \phi_{kj}^{-})/2$ and $V_{kj} = \phi_{kj}^{+} - \phi_{kj}^{-}$. The quantities $\phi_{kj}^{+}$ and $\phi_{kj}^{-}$ are the electrical boundary conditions at top ($z = z_{kj} + h_{kj}/2$) and at bottom ($z = z_{kj} - h_{kj}/2$) surfaces. The local z-axis of the $k_j$th face sublayer is situated at (for $k = a(+), b(-)$)

$$z_{kj} = \pm \frac{h_{kj} + h_c}{2} \pm \sum_{r=1}^{j-1} h_{kr}$$

with $h_{kj}$ being the thickness of the $k_j$th layer. Secondly, the axial component of the electrical field can be neglected because its contribution to the electromechanical energy is small if compared to that of the transverse one [5]. The two above assumptions imply a constant transverse electrical field within the $k_j$th piezoelectric sublayer

$$E_{3kj} = - \frac{\partial \phi_{kj}}{\partial z} = - \frac{V_{kj}}{h_{kj}} \quad (5)$$

### 2.2. Fractional derivative viscoelastic model

The one-dimensional constitutive equation introduced by Bagley and Torvik [1] is adopted in this work to describe the viscoelastic behavior of the core

$$\sigma_i + \tau^\alpha D^\alpha \sigma_i = \bar{c}_{ii} \left( \varepsilon_i + \tau^\alpha \frac{E_\infty}{E_0} D^\alpha \varepsilon_i \right) \quad (6)$$
where $\sigma_i$ and $\epsilon_i$ are axial ($i=1$) and shear ($i=5$) stress and strain. Since the core is supposed to be isotropic, the elastic constants are defined by: $\bar{c}_{11} = E_0/(1-v^2)$ and $\bar{c}_{55} = E_0/(2+2v)$. Furthermore, $E_0$ and $E_\infty$ are the relaxed and non-relaxed elastic moduli, $\tau$ is the relaxation time, $\alpha$ is the fractional order of the time derivative, and $\mathcal{D}^\alpha$ denotes the operator of fractional derivation of $\alpha$th order. The statements $0 < \alpha < 1$, $\tau > 0$ and $E_\infty > E_0$ fulfill the second law of thermodynamics.

This four-parameter fractional derivative model has been shown to be an effective tool to describe the weak frequency-dependence of most viscoelastic materials [1,7]. Its behavior in frequency domain is described between two asymptotic values: the static modulus of elasticity $E_0$ and the high-frequency limit value of the dynamic modulus $E_\infty$. The Fourier transform of Eq. (6) leads to the calculation of the complex modulus of elasticity, which is given by

$$E^*(\omega) = \frac{E_0 + E_\infty(i\omega\tau)^\alpha}{1 + (i\omega\tau)^\alpha} \quad (7)$$

The fractional operator $\mathcal{D}^\alpha$, appearing in the constitutive Eq. (6), can be approximated by several methods. One of them is the Grünwald definition, which is often adopted in literature since it is valid for all values of $\alpha$ and easy to implement numerically. The finite difference approximation of the Grünwald definition is given by

$$(D^\alpha f)_n \approx \frac{1}{\Delta t^\alpha} \sum_{j=0}^{N_1} A_{j+1} f_{n-j} \quad (8)$$

where $\Delta t$ is the time step increment of the numerical scheme (the function $f_n$ is approximated by $f(t_n)$, with $t_n = n\Delta t$), $N_1$ is the truncation number of the series, and $A_{j+1}$ represents the Grünwald coefficients given either in terms of the gamma function or by a recurrence formula

$$A_{j+1} = \frac{\Gamma(j-\alpha)}{\Gamma(-\alpha)\Gamma(j+1)} = \frac{j-\alpha-1}{j} A_j$$

## 2.3. Piezoelectric constitutive equations

The piezoelectric sublayers of the laminated faces are poled in the thickness direction with an electrical field applied parallel to this polarization. Such a configuration is characterized by the electromechanical coupling between the axial strain $\epsilon_1$ and the transverse electrical field $E_3$. The three-dimensional constitutive equations can be reduced to

$$\sigma_1 = \bar{c}_{11}\epsilon_1 - \bar{e}_{31}E_3 \quad (9a)$$

$$D_3 = \bar{e}_{31}\epsilon_1 + \bar{d}_{33}E_3 \quad (9b)$$

where $\sigma_1$ and $D_3$ are the axial stress and the transverse electrical displacement. Modified elastic, piezoelectric and dielectric constants are respectively given by

$$\bar{c}_{11} = c_{11} - \frac{c_{13}^2}{c_{33}}, \quad \bar{e}_{31} = e_{31} - \frac{c_{13}e_{33}}{c_{33}}, \quad \bar{d}_{33} = d_{33} + \frac{e_{33}^2}{c_{33}} \quad (10)$$

For elastic faces, the piezoelectric constants vanish. Moreover, if the material is isotropic, $\bar{c}_{11} = E/(1-v^2)$, where $E$ and $v$ are the elastic modulus and the Poisson's ratio.

## 2.4. Virtual strain energy

The virtual strain energy of the three-layer sandwich beam is the sum of the variations of the internal energy of each layer: $\delta U = \delta U_a + \delta U_b + \delta U_c$. This will be used in the finite element formulation in order to evaluate the internal force vector and stiffness matrix in the local coordinate system.

### 2.4.1. Viscoelastic core

For the core, the variation of the strain energy is given by

$$\delta U_c = \int_{\Omega_c} (\sigma_{1c} \delta \epsilon_{1c} + \sigma_{5c} \delta \epsilon_{5c}) dV \quad (11)$$

In order to take into account the fractional derivative viscoelastic behavior, let us introduce the internal variable $\bar{e}_{ic}$ such that

$$\sigma_{ic} = \bar{c}_{ii}^c \frac{E_\infty}{E_0} (\bar{e}_{ic} - \bar{e}_{ic}) \quad (12)$$

The constitutive equation (6) can then be rewritten as

$$\bar{e}_{ic} + \tau^\alpha D^\alpha \bar{e}_{ic} = \frac{E_\infty - E_0}{E_\infty} \epsilon_{ic} \quad (13)$$

This variable change implies that Eq. (13) contains only one fractional derivative term instead of two as in (6).

Using the Grünwald approximation (8) and noting that, by definition, $A_1 = 1$, relation (13) takes the following discretized form:

$$\bar{e}_{ic}^{n+1} = (1 - c_\alpha) \frac{E_\infty - E_0}{E_\infty} \epsilon_{ic}^{n+1} - c_\alpha \sum_{j=1}^{N_1} A_{j+1} \bar{e}_{ic}^{n+1-j} \quad (14)$$

where $c_\alpha$ is a dimensionless constant given by $c_\alpha = \tau^\alpha/(\tau^\alpha + \Delta t^\alpha)$.

Replacing this expression in Eq. (12), the axial and shear stresses in the core at time $t_{n+1}$ are given by

$$\sigma_{ic}^{n+1} = \bar{c}_{ii}^c \left\{ \left( 1 + c_\alpha \frac{E_\infty - E_0}{E_0} \right) \epsilon_{ic}^{n+1} + c_\alpha \frac{E_\infty}{E_0} \sum_{j=1}^{N_1} A_{j+1} \bar{e}_{ic}^{n+1-j} \right\} \quad (15)$$

The variation of the strain energy of the core is thus composed of two terms: the first one can be associated with a modified internal force vector or stiffness matrix, and the second one, which depends on the anelastic deformation history, can be shifted to the right-hand side of the governing equation modifying in this way the external force vector.

It should be stated that the Grünwald coefficients in Eq. (14), which are strictly decreasing when $j$ increases, describe the fading memory phenomena. In other words, the behavior of the viscoelastic material at a given time step depends
more strongly on the recent time history than on the distant one.

For an elastic material, since the constant $c_\alpha$ is zero, using strain relations (3) and integrating over the cross section, Eq. (11) gives

$$ \delta U_c = \int_0^L [\bar{e}_{11}^c (A_c \epsilon_c \delta \epsilon_c + I_c \kappa_c \delta \kappa_c) + \bar{e}_{55}^c A_c \gamma_c \delta \gamma_c] \, dx \quad (16) $$

where $A_c = bh_c$ and $I_c = bh_c^3/12$ are the cross section area and moment of inertia of the core. Moreover, $L$ and $b$ are the total length and width of the beam.

### 2.4.2. Piezoelectric laminated faces

For the piezoelectric/elastic laminated faces ($k = a, b$), the virtual internal energy is

$$ \begin{aligned} \delta U_k &= \sum_{j=1}^{N_k} \int_{\Omega_{k,j}} (\sigma_{1k_j} \delta \epsilon_{1k} - D_{3k_j} \delta E_{3k_j}) \, d\Omega \\ &= \delta U_k^M + \delta U_k^{ME} + \delta U_k^{EM} + \delta U_k^E \end{aligned} \quad (17) $$

where $\delta U_k^M$ and $\delta U_k^E$ are the virtual mechanical and dielectrical internal energies, and $\delta U_k^{ME}$ and $\delta U_k^{EM}$ comprise the virtual piezoelectric internal energy, i.e., the electromechanical coupling.

Using strain relations (3), piezoelectric constitutive equations (9), and electrical field expression (5), each term in the right-hand side of Eq. (17) can be evaluated as described below.

The virtual mechanical internal energy can be written as

$$ \delta U_k^M = b \int_0^L [A^k \epsilon_k \delta \epsilon_k + B^k (\kappa_k \delta \epsilon_k + \epsilon_k \delta \kappa_k) + D^k \kappa_k \delta \kappa_k] \, dx \quad (18) $$

where $A^k$ is the extensional stiffness, $D^k$ is the bending stiffness, and $B^k$ is the bending-extensional coupling stiffness defined as

$$ [A^k, B^k, D^k] = \sum_{j=1}^{N_k} \bar{e}_{11}^{kj} \int_{z_{k_j} - h_{k_j} / 2}^{z_{k_j} + h_{k_j} / 2} [1, (z - z_k), (z - z_k)^2] \, dz \quad (19) $$

The three last terms in the Eq. (17) are related to the electromechanical coupling and the dielectrical quantities in the virtual internal energy of the faces. These terms are defined as

$$ \delta U_k^{ME} = b \sum_{j=1}^{N_k} \int_0^L \bar{e}_{31}^{kj} V_{k_j} [\delta \epsilon_k + (z_{k_j} - z_k) \delta \kappa_k] \, dx \quad (20a) $$

$$ \delta U_k^{EM} = b \sum_{j=1}^{N_k} \int_0^L \bar{e}_{31}^{kj} [\epsilon_k + (z_{k_j} - z_k) \kappa_k] \delta V_{k_j} \, dx \quad (20b) $$

$$ \delta U_k^E = -b \sum_{j=1}^{N_k} \int_0^L \bar{d}_{33}^{kj} \frac{V_{k_j}}{h_{k_j}} \delta V_{k_j} \, dx \quad (20c) $$

It is easy to see in Eqs. (20a) and (20b) that the electromechanical coupling is taken into account by means of the piezoelectric constant $\bar{e}_{31}^{kj}$. This electromechanical coupling

describes the so-called inverse and direct piezoelectric effects.

An *actuator configuration* consists of applying an electrical field to a piezoelectric ceramic in order to induce a mechanical deformation. This electrical field is imposed here by means of an external force, which is written as a function of the voltage applied to the piezoelectric patch. Consequently, the variation of the difference of electrical potential of the $k_j$th lamina vanishes and the above mentioned external force is extracted from Eq. (20a) (see [5] for more details).

The case where an electrical field is induced by a mechanical deformation in the piezoelectric material corresponds to a *sensor configuration*. The unknown is the difference of electric potential of the $k_j$th piezoelectric layer. Using Eqs. (20b) and (20c), it can be shown that the virtual electromechanical internal energy of the $k_j$th lamina corresponds to an increase of the internal force vector or stiffness matrix of the beam due to the direct piezoelectric effect [5].

## 3. Co-rotational finite element formulation

The co-rotational finite element formulation provides a simple method for large displacement analysis. The key idea is to separate the motion of each finite element in a rigid body part and a deformational part, which is small relative to a local coordinate system [3,6]. The equations of motion are defined in the global frame while strains are measured in the co-rotational coordinate system of the element. In this study, the linear sandwich beam theory is used in a coordinate system, which rotates and translates with the element but do not deform with it (Fig. 2).

Stiffness and mass matrices of the individual elements are constructed in the element coordinate system, and then transformed to the global coordinate system using standard procedures [2]. The resulting equations of motion are the same as those typically arising in non-linear

<figure><img src="image_3.png" /><figcaption>Fig. 2. Co-rotational sandwich beam coordinate system.</figcaption></figure>
structural dynamics. In this section, we first derive relations between the local and global virtual displacement vectors. Then we define the internal force vector as well as the tangent stiffness matrix in the global coordinate system. Finally, the resolution method, based on the Newmark direct integration scheme and the Newton-Raphson algorithm, is briefly described.

### 3.1. Virtual displacement vectors

The linear finite element formulation of the sandwich beam composed of elastic/piezoelectric laminated faces and a viscoelastic core is thoroughly described in [5]. We just recall that the displacements are discretized with linear (axial displacement) and cubic (deflection) shape functions. We propose in this section to determine relations between the local and global virtual displacement vectors. These relations will be used to compute the internal force vector and tangent stiffness matrix.

The vectors of global and local displacements are respectively defined by

$$ \mathbf{q} = [\bar{u}_1 \quad w_1 \quad \beta_1 \quad \tilde{u}_1 \quad | \quad \bar{u}_2 \quad w_2 \quad \beta_2 \quad \tilde{u}_2]^T \quad (21) $$

and

$$ \hat{\mathbf{q}} = [\hat{u}_1 \quad \hat{w}'_1 \quad \hat{u}_1 \quad | \quad \hat{u}_2 \quad \hat{w}'_2 \quad \hat{u}_2]^T \quad (22) $$

Using the notations defined in Fig. 2, the displacements and rotations for nodes 1 and 2 in the local coordinate system $(\hat{\mathbf{x}}, \hat{\mathbf{z}})$ are given, in terms of global quantities, by

$$ \begin{aligned} \hat{u}_{c1} &= 0, & \hat{\theta}_{c1} &= \theta_{c1} - \phi, & \hat{w}'_1 &= \beta_1 - \phi \\ \hat{u}_{c2} &= \ell - \ell_0, & \hat{\theta}_{c2} &= \theta_{c2} - \phi, & \hat{w}'_2 &= \beta_2 - \phi \end{aligned} \quad (23) $$

where $\ell_0$ and $\ell$ denote the initial and current lengths of the element.

In these equations, local components of axial displacements and rotations of the core can be written in terms of the local degrees of freedom (for node $i = 1, 2$)

$$ \hat{u}_{c_i} = \hat{u}_i + h_1 \hat{w}'_i \quad (24) $$

$$ \hat{\theta}_{c_i} = -\frac{1}{h_c} \hat{u}_i - \left(\frac{h_2}{h_c} - 1\right) \hat{w}'_i \quad (25) $$

where $h_1 = \tilde{h}/4$ and $h_2 = h_c + \bar{h}$.

Using previous equations and geometrical considerations, it can be shown that the variations of $\ell$ and $\phi$ are given by

$$ \delta\ell = \mathbf{r}^T \delta\mathbf{q} \quad (26) $$

$$ \delta\phi = \frac{1}{\ell} \mathbf{z}^T \delta\mathbf{q} \quad (27) $$

where the following notations are introduced:

$$ \begin{aligned} \mathbf{r} &= [-c \quad -s \quad -h_1 c \quad 0 \quad c \quad s \quad h_1 c \quad 0]^T \\ \mathbf{z} &= [s \quad -c \quad h_1 s \quad 0 \quad -s \quad c \quad -h_1 s \quad 0]^T \end{aligned} $$

with $c = \cos(\phi + \phi_0)$ and $s = \sin(\phi + \phi_0)$.

The virtual local displacements are obtained through differentiation of the components of the local degrees of

freedom vector. Hence, the transformation matrix $\mathbf{T}$, which connects local and global virtual displacement vectors

$$ \delta\hat{\mathbf{q}} = \mathbf{T}\delta\mathbf{q} \quad (28) $$

is given by

$$ \mathbf{T} = \frac{1}{\ell} \begin{bmatrix} h_1 s & -h_1 c & -\ell h_1 + h_1^2 s & 0 & -h_1 s & h_1 c & -h_1^2 s & 0 \\ -s & c & \ell - h_1 s & 0 & s & -\ell c & h_1 s & 0 \\ h_2 s & -h_2 c & h_1 h_2 s & \ell & -h_2 s & h_2 c & -h_1 h_2 s & 0 \\ -\ell c + h_1 s & -\ell s - h_1 c & -\ell h_1 c + h_1^2 s & 0 & \ell c - h_1 s & \ell s + h_1 c & \ell h_1 (c-1) - h_1^2 s & 0 \\ -s & c & -h_1 s & 0 & s & -c & \ell + h_1 s & 0 \\ h_2 s & -h_2 c & h_1 h_2 s & 0 & -h_2 s & h_2 c & -h_1 h_2 s & 1 \end{bmatrix} $$

### 3.2. Internal force vector

The local and global internal force vectors are calculated by using the internal virtual work in both local and global systems

$$ \delta U = \delta \mathbf{q}^T \mathbf{F} = \delta \hat{\mathbf{q}}^T \hat{\mathbf{F}} = \delta \mathbf{q}^T \mathbf{T}^T \hat{\mathbf{F}} \quad (29) $$

From the internal virtual energy in the local coordinate system $\delta U = \delta U_a + \delta U_b + \delta U_c$ and using linear and cubic shape functions, the local internal force vector

$$ \hat{\mathbf{F}} = \hat{\mathbf{F}}_a + \hat{\mathbf{F}}_b + \hat{\mathbf{F}}_c = [\bar{N}_1 \quad M_1 \quad \tilde{N}_1 \quad \bar{N}_2 \quad M_2 \quad \tilde{N}_2]^T \quad (30) $$

can be evaluated analytically. It can be noted that the resultant forces $\bar{N}$ and $\tilde{N}$ are associated with mean and relative axial displacements.

Since Eq. (29) must hold for any arbitrary $\delta\mathbf{q}$, the global internal force vector is simply given by

$$ \mathbf{F} = \mathbf{T}^T \hat{\mathbf{F}} \quad (31) $$

### 3.3. Tangent stiffness matrix

The global tangent stiffness matrix $\mathbf{K}_T$ is obtained by differentiation of Eq. (31)

$$ \delta \mathbf{F} = \mathbf{T}^T \delta \hat{\mathbf{F}} + \delta \mathbf{T}^T \hat{\mathbf{F}} = \mathbf{K}_M \delta \mathbf{q} + \mathbf{K}_G \delta \mathbf{q} = \mathbf{K}_T \delta \mathbf{q} \quad (32) $$

where $\mathbf{K}_M$ and $\mathbf{K}_G$ are material and geometrical stiffness matrices.

The global material stiffness matrix $\mathbf{K}_M$ is computed by

$$ \mathbf{K}_M = \mathbf{T}^T \hat{\mathbf{K}} \mathbf{T} \quad (33) $$

where the $6 \times 6$ local stiffness matrix is such that

$$ \delta \hat{\mathbf{F}} = \hat{\mathbf{K}} \delta \hat{\mathbf{q}} \quad (34) $$

Using previous definitions for $\hat{\mathbf{q}}$ and $\hat{\mathbf{F}}$, this matrix can be evaluated analytically. For example, if the three layers are elastic homogeneous, we have

$$ \hat{\mathbf{K}} = \begin{bmatrix} \hat{K}_{11} & \hat{K}_{12} & \hat{K}_{13} & -\hat{K}_{11} & -\hat{K}_{12} & -\hat{K}_{13} \\ \hat{K}_{22} & \hat{K}_{23} & -\hat{K}_{12} & \hat{K}_{25} & -\hat{K}_{23} & -\hat{K}_{13} \\ \hat{K}_{33} & -\hat{K}_{13} & -\hat{K}_{23} & \hat{K}_{36} & \hat{K}_{12} & \hat{K}_{13} \\ & & & \hat{K}_{11} & \hat{K}_{12} & \hat{K}_{22} & \hat{K}_{23} \\ & & & & \hat{K}_{33} & & \hat{K}_{36} \end{bmatrix} \quad (35) $$
with the following components:

$$ \hat{K}_{11} = \frac{1}{\ell_0} (\bar{c}_{11}^a A_a + \bar{c}_{11}^b A_b + \bar{c}_{11}^c A_c) $$

$$ \hat{K}_{12} = \frac{h_1}{\ell_0} \bar{c}_{11}^c A_c $$

$$ \hat{K}_{13} = \frac{1}{2\ell_0} (\bar{c}_{11}^a A_a - \bar{c}_{11}^b A_b) $$

$$ \hat{K}_{22} = \frac{4}{\ell_0} (\bar{c}_{11}^a I_a + \bar{c}_{11}^b I_b) + \frac{4(h_2 - h_c)^2}{\ell_0 h_c^2} \bar{c}_{11}^c I_c + \frac{4h_1^2}{\ell_0} \bar{c}_{11}^c A_c + \frac{2h_2^2 \ell_0}{15h_c^2} \bar{c}_{55}^c A_c $$

$$ \hat{K}_{23} = \frac{h_2 - h_c}{\ell_0 h_c^2} \bar{c}_{11}^c I_c + \frac{h_2 \ell_0}{12h_c^2} \bar{c}_{55}^c A_c $$

$$ \hat{K}_{25} = \frac{2}{\ell_0} (\bar{c}_{11}^a I_a + \bar{c}_{11}^b I_b) + \frac{2(h_2 - h_c)^2}{\ell_0 h_c^2} \bar{c}_{11}^c I_c + \frac{2h_1^2}{\ell_0} \bar{c}_{11}^c A_c - \frac{h_2^2 \ell_0}{30h_c^2} \bar{c}_{55}^c A_c $$

$$ \hat{K}_{33} = \frac{1}{4\ell_0} (\bar{c}_{11}^a A_a + \bar{c}_{11}^b A_b) + \frac{1}{\ell_0 h_c^2} \bar{c}_{11}^c I_c + \frac{\ell_0}{3h_c^2} \bar{c}_{55}^c A_c $$

$$ \hat{K}_{36} = -\frac{1}{4\ell_0} (\bar{c}_{11}^a A_a - \bar{c}_{11}^b A_b) - \frac{1}{\ell_0 h_c^2} \bar{c}_{11}^c I_c + \frac{\ell_0}{6h_c^2} \bar{c}_{55}^c A_c $$

Moreover, the evaluation of δB gives the expression of the global geometrical tangent stiffness matrix

$$ \mathbf{K}_G = \frac{1}{\ell} \mathbf{z} \mathbf{z}^T \overline{\mathbf{N}_2} - \frac{1}{\ell^2} (\mathbf{r} \mathbf{z}^T + \mathbf{z} \mathbf{r}^T) \times [h_1(\overline{\mathbf{N}_1} + \overline{\mathbf{N}_2}) - (M_1 + M_2) + h_2(\overline{\mathbf{N}_1} + \overline{\mathbf{N}_2})] \quad (36) $$

### 3.4. Resolution algorithm

An incremental-iterative method based on the Newmark direct integration scheme (with $\beta = 1/4$ and $\gamma = 1/2$) and Newton-Raphson algorithm is employed here. Assuming the equilibrium configuration at time $t_n$ is known, one searches the solution $q^{n+1}$ at time $t_{n+1}$ such that the following non-linear equation is satisfied:

$$ \mathbf{R}^{n+1} = \mathbf{M}\mathbf{q}^{n+1} + \mathbf{F}(\mathbf{q}^{n+1}) - \mathbf{G}^{n+1} - \tilde{\mathbf{G}}^{n+1} - \overline{\mathbf{G}}^{n+1} = \mathbf{0} \quad (37) $$

where $\mathbf{R}$ is the residual vector, $\mathbf{M}$ is the mass matrix, $\mathbf{F}$ the internal force vector, $\mathbf{G}$ the mechanical external load, $\tilde{\mathbf{G}}$ the electrical force associated with the actuator configuration for the faces and $\overline{\mathbf{G}}$ the dissipative force associated with the viscoelastic behavior of the core. All these quantities are defined in the global reference system using transformation matrices and their local definitions given in [5] and thoroughly explained in [4].

It is well known that in a Newmark-Newton algorithm, an iteration loop on equilibrium is performed within the time loop. The non-linear equation (37) is then solved at each time step. Moreover, at each iteration, the displacement correction $\Delta\mathbf{q}$ is computed using the linearized equation $\mathbf{S}(\mathbf{q}^{n+1})\Delta\mathbf{q} = -\mathbf{R}^{n+1}$, where the iteration matrix is defined by $\mathbf{S} = \mathbf{K}_T + \mathbf{M}/(\beta\Delta t^2)$.

## 4. Numerical examples

### 4.1. Flexible robot arm

This simulation is concerned with the repositioning of a sandwich flexible beam rotating horizontally about a vertical

axis passing through one end. This classical example was studied for example in [6,8] for a single layer beam. Here, the beam is composed of a viscoelastic core (ISD112 at 27 °C) constrained by symmetrical elastic faces (Aluminum). Geometry data and material properties of the structure are shown in Fig. 3 and Table 1, respectively. Model parameters for the ISD112 at 27 °C are identified in [4], assuming that the Poisson's ratio is frequency-independent.

The robot arm is first repositioned to an angle of $\pi/2$ rad from its initial position. This is achieved by prescribing the rotation angle $\psi(t)$ as a linear function of time, as shown in Fig. 3. The beam is discretized by a regular mesh with 10 elements. The study is performed up to 200 ms with a fixed time step $\Delta t = 0.05$ ms. The sequence of motion during the repositioning stage is depicted in Fig. 4a. Once $\psi(t) = \pi/2$ rad for all time $t \ge t_1 = 20$ ms, the robot arm undergoes finite vibrations as shown in Fig. 4b.

Fig. 5a illustrates the attenuation of the tip elastic displacements in x and z directions due to the constrained layer damping treatment. The damped oscillations for x component can also be observed in the phase-space diagram (Fig. 5b).

### 4.2. Slider-crank mechanism

This second example corresponds to a slider-crank mechanism with a rigid crank OA, a connecting sandwich beam AB and a slider block located at B. The sandwich beam is composed by a viscoelastic core ($h_c = 0.2$ mm)

<figure><img src="image_4.png" /><figcaption>Fig. 3. Flexible sandwich robot arm.</figcaption></figure>

<table><caption>Table 1 Mechanic and piezoelectric characteristics of the sandwich beam</caption><tr><td>Aluminum</td><td>ρ = 2690 kg/m³, ν = 0.345, E = 70.3 GPa</td></tr><tr><td>ISD112</td><td>ρ = 1600 kg/m³, ν = 0.5<br>E₀ = 1.5 MPa, E∞ = 69.95 MPa<br>α = 0.7915, τ = 1.4052 × 10⁻² ms</td></tr><tr><td>PZT5H</td><td>ρ = 7500 kg/m³, c₁₁ = c₃₃ = 126 GPa<br>c₁₃ = 84.1 GPa, e₃₁ = -6.5 C/m²<br>e₃₃ = 23.3 C/m², d₃₃ = 1.3 × 10⁻⁸ F/m</td></tr></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Deflected shapes of the sandwich robot arm: (a) repositioning sequence up to $\psi = \pi/2$ rad and (b) free vibration about $\psi = \pi/2$ rad.</figcaption></figure>

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

<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Time histories for x and z tip displacement components and (b) phase-space diagram for x component.</figcaption></figure>

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

constrained by symmetrical elastic faces ($h_{a1} = h_{b1} = 1$ mm) with piezoelectric patches ($h_{a2} = h_{b2} = 0.5$ mm). The dimensions and material properties of the slider-crank mechanism are given in Fig. 6 and Table 1, respectively. The crank is driven by a constant angular speed $\dot{\psi} = 0.15$ rad/ms. The points O, A and B are initially aligned with A between O and B. Moreover, the initial velocity and acceleration of the mechanism are calculated by using kinematics of rigid mechanism.

The beam AB is discretized by a regular mesh of 10 elements. The study is performed up to 60 ms with a fixed time step $\Delta t = 0.05$ ms. The fractional operator is approximated using the whole history in the Grünwald series. Moreover, a sensor voltage feedback control law is used.

<figure><img src="image_9.png" /><figcaption>Fig. 6. Slider-crank mechanism.</figcaption></figure>

The top piezoelectric patch works as an actuator and the bottom one as a sensor. The time derivative of the voltage measured in the sensor is used for the feedback signal and then applied to the top piezoelectric patch with a constant gain $K_d = -100$ ms, as $V_A = -K_d \dot{V}_S$.

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

<figure><img src="image_11.png" /><figcaption>Fig. 7. Deflection versus crank angle and phase-space diagram for the midpoint of the sandwich beam.</figcaption></figure>
The transverse deflection of the midpoint of the connect-
ing beam AB is measured perpendicularly to the straight
line connecting points A and B. In Fig. 7, deflection versus
crank angle and phase-space diagram are plotted in order
to illustrate the attenuation of the elastic deformation
due to the active/passive damping treatment.

**5. Conclusions**

A co-rotational approach for non-linear dynamic analy-
sis of flexible linkage mechanisms with active constrained
layer damping treatment is proposed in this work. An ori-
ginal co-rotational sandwich beam finite element, which
allows large displacements and rotations, is developed for
this purpose. The sandwich beam is composed of a visco-
elastic core, covered by elastic/piezoelectric laminated
faces. No electrical degrees of freedom are required in the
finite element formulation, since the electromechanical cou-
pling is taken into account by means of an augmentation of
the stiffness of the piezoelectric layers in the sensor case.
Concerning the viscoelastic damping, a four-parameter
fractional derivative model is used to describe the dissipa-
tive behavior of the core. For the numerical examples, an
incremental-iterative method based on the Newmark direct
integration scheme and the Newton-Raphson algorithm is
employed. Two particularities are associated with the
numerical resolution method. The first one concerns the
approximation of fractional operators by introducing an
external dissipative force which involves the storage of
the displacement history. Obviously, this approach reduces
the numerical costs related to the non-locality of the viscoe-
lactic behavior. The second one concerns the non linear
algorithm which requires the evaluation of the internal
force vector and tangent stiffness matrix associated with
the sandwich element. As described in the paper, these
two quantities are calculated analytically in the global
coordinate system. Finally, two numerical examples with

large displacements are proposed in order to test the co-
rotational active/passive sandwich beam element devel-
oped in this work. The results illustrate the effectiveness
of the proposed approach. An evaluation of the active con-
strained layer damping treatment performances in the con-
text of non-linear dynamics will be investigated in a future
work. Moreover, an experimental validation of the present
formulation would be interesting particularly to demon-
strate the effectiveness of the viscoelastic fractional deriva-
tive model.

**References**

[1] Bagley RL, Torvik PJ. Fractional calculus – a different approach to the analysis of viscoelastically damped structures. AIAA J 1983;21: 741–8.

[2] Crisfield MA. Non-linear finite element analysis of solids and structures, Volume 1: essentials. Chichester: Wiley; 1991.

[3] Elkaranshawy HA, Dokainish MA. Corotational finite element analysis of planar flexible multibody systems. Comput Struct 1995;54:881–90.

[4] A.C. Galucio, Atténuation des réponses transitoires par traitement hybride piézoélectrique/viscoélastique en utilisant un modèle à dérivées fractionnaires, PhD thesis; 2004 [in French].

[5] Galucio AC, Deü J-F, Ohayon R. A fractional derivative viscoelastic model for hybrid active/passive damping treatments in time domain – application to sandwich beams. J Intell Mater Syst Struct 2005;16: 33–45.

[6] Hsiao K-M, Jang J-Y. Dynamic analysis of planar flexible mecha-
nisms by corotational formulation. Comput Methods Appl Mech
Eng 1991;87:1–14.

[7] Pritz T. Analysis of four-parameter fractional derivative model of real solid materials. J Sound Vib 1996;195:103–15.

[8] Simo JC, Vu-Quoc L. On the dynamics of flexible beams under large overall motions – the plane case: Parts I and II. J Appl Mech 1986;53: 849–63.

[9] Trindade MA, Benjeddou A, Ohayon R. Finite element modelling of hybrid active–passive vibration damping of multilayer piezoelectric sandwich beams – part I: Formulation; part II: System analysis. Int J Numer Methods Eng 2001;51:835–64.

[10] Wasfy TM, Noor AK. Computational strategies for flexible multi-body systems. Appl Mech Rev 2003;56:553–632.