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

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 Grunwald approximation of fractional derivatives, and the Newton-Raphson algorithm.


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

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.

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 Þh i ðx; tÞ ð 1aÞ u zi ðx; z; tÞ ¼ wðx; tÞ ð 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 h i are the axial displacement of the center line and the fiber rotation of each layer; and w is the transverse displacement ( Fig. 1).
Let us introduce the mean and relative axial displacements given by u ¼ ðu a þ u b Þ=2 andũ ¼ u a À u b . Euler-Bernoulli hypotheses for the faces lead to h k ¼ w 0 for k ¼ a; b with ðÁÞ 0 ¼ oðÁÞ=ox. 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 0 and the above defined variables u andũ as From (1) and (2), and taking the hypothesis of plane stress state into account, we can write the axial strain of the ith layer e 1i and the shear strain of the core e 5c as follows: where i and j i are the membrane strain and curvature of the ith layer, and c 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 / kj ðx; z; tÞ ¼ / kj ðx; tÞ þ ðz À z kj Þ V kj ðx; tÞ h kj ð4Þ The quantities / þ kj and / À 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ðÀÞ) 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

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 where r i and e 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: c 11 ¼ E 0 =ð1 À m 2 Þ and c 55 ¼ E 0 =ð2 þ 2mÞ. Furthermore, E 0 and E 1 are the relaxed and non-relaxed elastic moduli, s is the relaxation time, a is the fractional order of the time derivative, and D a denotes the operator of fractional derivation of ath order. The statements 0 < a < 1, s > 0 and E 1 > 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 1 . The Fourier transform of Eq. (6) leads to the calculation of the complex modulus of elasticity, which is given by The fractional operator D a , 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 a and easy to implement numerically. The finite difference approximation of the Grü nwald definition is given by where Dt is the time step increment of the numerical scheme (the function f n is approximated by f ðt n Þ, with t n ¼ nDt), N t 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

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 e 1 and the transverse electrical field E 3 . The three-dimensional constitutive equations can be reduced to where r 1 and D 3 are the axial stress and the transverse electrical displacement. Modified elastic, piezoelectric and dielectric constants are respectively given by For elastic faces, the piezoelectric constants vanish. Moreover, if the material is isotropic, c 11 ¼ E=ð1 À m 2 Þ, where E and m are the elastic modulus and the Poisson's ratio.

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: dU ¼ dU a þ dU b þ dU 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.

Viscoelastic core
For the core, the variation of the strain energy is given by In order to take into account the fractional derivative viscoelastic behavior, let us introduce the internal variable e ic such that The constitutive equation (6) can then be rewritten as 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: where c a is a dimensionless constant given by c a ¼ s a =ðs a þ Dt a Þ.
Replacing this expression in Eq. (12), the axial and shear stresses in the core at time t nþ1 are given by 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 a is zero, using strain relations (3) and integrating over the cross section, Eq. (11) gives where A c ¼ bh c and I c ¼ bh 3 c =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.

Piezoelectric laminated faces
For the piezoelectric/elastic laminated faces ðk ¼ a; bÞ, the virtual internal energy is where dU M k and dU E k are the virtual mechanical and dielectrical internal energies, and dU ME k and dU EM k 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 where A k is the extensional stiffness, D k is the bending stiffness, and B k is the bending-extensional coupling stiffness defined as 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 dU E k ¼ Àb It is easy to see in Eqs. (20a) and (20b) that the electromechanical coupling is taken into account by means of the piezoelectric constant e kj 31 . 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 kth lamina corresponds to an increase of the internal force vector or stiffness matrix of the beam due to the direct piezoelectric effect [5].

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 Fig. 2. Co-rotational sandwich beam coordinate system. 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.

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 Using the notations defined in Fig. 2, the displacements and rotations for nodes 1 and 2 in the local coordinate system ðx;ẑÞ are given, in terms of global quantities, bŷ where ' 0 and ' 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) Using previous equations and geometrical considerations, it can be shown that the variations of ' and / are given by where the following notations are introduced: The virtual local displacements are obtained through differentiation of the components of the local degrees of freedom vector. Hence, the transformation matrix T, which connects local and global virtual displacement vectors dq ¼ Tdq ð28Þ is given by

Internal force vector
The local and global internal force vectors are calculated by using the internal virtual work in both local and global systems From the internal virtual energy in the local coordinate system dU ¼ dU a þ dU b þ dU c and using linear and cubic shape functions, the local internal force vector can be evaluated analytically. It can be noted that the resultant forces N and e N are associated with mean and relative axial displacements.
Since Eq. (29) must hold for any arbitrary dq, the global internal force vector is simply given by

Tangent stiffness matrix
The global tangent stiffness matrix K T is obtained by differentiation of Eq. (31) where K M and K G are material and geometrical stiffness matrices. The global material stiffness matrix K M is computed by where the 6 · 6 local stiffness matrix is such that Using previous definitions forq and b F, this matrix can be evaluated analytically. For example, if the three layers are elastic homogeneous, we have with the following components:

Resolution algorithm
An incremental-iterative method based on the Newmark direct integration scheme (with b ¼ 1=4 and c ¼ 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 of motion is satisfied: where R is the residual vector, M is the mass matrix, F the internal force vector, G the mechanical external load, e G the electrical force associated with the actuator configuration for the faces and 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 Dq is computed using the linearized equation Sðq nþ1 ÞDq ¼ ÀR nþ1 , where the iteration matrix is defined by S ¼ K T þ M=ðbDt 2 Þ.

Flexible robot arm
This simulation is concerned with the repositioning of a sandwich flexible beam rotating horizontally about a verti-cal 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 frequencyindependent.
The robot arm is first repositioned to an angle of p=2 rad from its initial position. This is achieved by prescribing the rotation angle wð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 Dt ¼ 0:05 ms. The sequence of motion during the repositioning stage is depicted in Fig. 4a. Once wðtÞ ¼ p=2 rad for all time t P 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).

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Þ   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 _ w ¼ 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 Dt ¼ 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.
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  The transverse deflection of the midpoint of the connecting 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.

Conclusions
A co-rotational approach for non-linear dynamic analysis of flexible linkage mechanisms with active constrained layer damping treatment is proposed in this work. An original co-rotational sandwich beam finite element, which allows large displacements and rotations, is developed for this purpose. The sandwich beam is composed of a viscoelastic core, covered by elastic/piezoelectric laminated faces. No electrical degrees of freedom are required in the finite element formulation, since the electromechanical coupling 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 dissipative 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 viscoelactic 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 corotational active/passive sandwich beam element developed in this work. The results illustrate the effectiveness of the proposed approach. An evaluation of the active constrained layer damping treatment performances in the context of non-linear dynamics will be investigated in a future work. Moreover, an experimental validation of the present formulation would be interesting particularly to demonstrate the effectiveness of the viscoelastic fractional derivative model.