Kelvin Viscoelasticity and Lagrange Multipliers Applied to the Simulation of Nonlinear Structural Vibration Control

This study proposes a new pure numerical way to model mass / spring / damper devices to control the vibration of truss structures developing large displacements. It avoids the solution of local differential equations present in traditional convolution approaches to solve viscoelasticity. The structure is modeled by the geometrically exact Finite Element Method based on positions. The introduction of the device's mass is made by means of Lagrange multipliers that imposes its movement along the straight line of a finite element. A pure numerical Kelvin/Voigt like rheological model capable of nonlinear large deformations is originally proposed here. It is numerically solved along time to accomplish the damping parameters of the device. Examples are solved to validate the formulation and to show the practical possibilities of the proposed technique

Latin American Journal of Solids and Structures 13 (2016) 964-991 structure.In particular, the formulations based on the Finite Element Method are very accepted for solving this kind of problem.Updated Lagrangian and co-rotational formulations can be seen, for example, in Bathe et al (1975), Simo and Vu-Quoc (1986), Armero (2006), Jelenic and Crisfield (2001) and Yoo et al. (2007).Total Lagrangian formulations based on positions and unconstrained vectors can be seen in Coda et al. (2013), Reis and Coda (2014) and Coda and Paccola (2008).The sole advantage of updated Lagrangian and co-rotational formulations is that they are extensions of linear formulations and, as a consequence, various well known literature solutions can be used to improve results.This advantage is also a disadvantage as in updated formulations finite rotations must be linearized and the current stress must be updated using the Kirchhoff-Jaumann formulae, while for the above mentioned total Lagrangian formulation none of these tricks are necessary.
Many lightweight structures, responsible for mechanical support of the designed object, consist of structural elements loaded predominantly by normal force.This kind of structure is usually called truss (Greco et al. (2006)).The finite element used to model trusses do not include bending stiffness bringing difficulties for modeling vibration control devices of the type mass / spring / damper in nonlinear analysis.In practice, this type of device has a mass which must slide over a surface (over a finite element) and its contact force is orthogonal to the member, resulting in bending stresses.Moreover, the dampers used in such devices can be represented by viscous constitutive models that, associated with their stiffness, result in viscoelastic models similar to the Kelvin/Voigt one.The treatment of these linear and nonlinear viscoelastic formulations is usually done by more or less complicated traditional convolution techniques or decaying functions as explained by Lemaitre and Chaboche (1990), Lemaitre (2001), Lima et al. (2015), Holzapfel (1996), Simo (1987), Huber and Tsakmakis (2000) and Petiteau et al (2013) used or not in vibration control (Marko et al. (2006) and Clough and Penzien (1985)).There are no significant differences regarding the traditional rheological viscoelastic approach of materials in all consulted bibliography.They use the same procedure, convolution, to solve different viscous materials, developing small or large strains.
In this study, it is presented: (i) a Total Lagrangian truss finite element for geometric nonlinear dynamic modeling of two-and three-dimensional structures.(ii) a technique based on Lagrange multipliers for the consideration of passive mass / spring / damper vibration control mechanisms and (iii) an innovative approach to the consideration of viscous damping in the form of the rheological model of Kelvin/Voigt adapted to the Green nonlinear strain derived from linear applications (Mesquita and Coda (2002), Mesquita and Coda (2003) and Mesquita and Coda (2007)).
The main contribution of this study is item (iii) that introduces a pure numerical way to consider viscous material behavior in dynamic or static applications.The proposed formulation avoids the analytical solution of time differential equations for each considered material law and the use of its solution as a residual in a convolution technique.Therefore, differently from the existent formulations, the proposed strategy can be easily extended to be used in any viscoelastic constitutive relation, linear or not.
The work is organized as follows.Section 2 presents the Energy Stationary Principle on which the formulation is based.It also presents each part of the total energy necessary to the development of the numerical process, detaching viscous damping and its numerical treatment, the kernel of the work.Section 3 presents the equilibrium equation and the organization of the numerical solution process.Section 4 presents the mass / spring / damper arrangement for large strain applications Latin American Journal of Solids and Structures 13 (2016) 964-991 making use of Lagrange multipliers.Representative examples are presented in section 5 to validate the proposed formulation and describing its practical possibilities.Conclusions are presented in section 6.

TOTAL MECHANICAL ENERGY AND STATIONARY PRINCIPLE
The principle of stationary energy indicates that the equilibrium of mechanical systems occurs when the variation of the total mechanical energy is zero.When it comes to dynamic applications, using the D'Alambert principle, the dynamic equilibrium is meant the equation of motion.
The principle of stationarity is usually employed when dissipative forces do not appear because, in this case, system energy level falls continuously over time.However, it is possible to apply the principle of stationarity considering a larger energy system than the system studied, in which the dissipated energy is part of it.Thus, the total mechanical energy to be considered in this study is written as.
Where U is the strain energy K is the kinetic energy, P is the potential of external forces (considered conservative) and Q is the dissipative potential resulting from the considered damping device.The total energy variation can be expressed as: In truss structure analysis it is usual to consider external forces applied to the nodes of the structure, as well as concentrated mass at nodes, otherwise it would be admitting the existence of bending stresses in structural elements, which is not possible by the mechanical definition of truss elements.Thus, the primary variables of the problem are nodal positions.
Therefore, to establish the adopted finite element formulation it is necessary to write equation (2) as a function of nodal position of the structure.

Kinetic Energy and External Forces Potential
Using index notation, the potential of external applied forces is given by: in which the external forces i F  are considered conservative, i.e., they do not depend on positions.In equation (3) i represents direction and  truss nodes on which the forces are applied or po- sitions measured.The variation of external forces potential, taken in relation to nodal positions, is written as: in which   assumes zero when node  is not the same as node  , otherwise it assumes one.
Latin American Journal of Solids and Structures 13 (2016) 964-991 One writes the kinetic energy for the analyzed structure as: where the over dot indicates time derivative, fe stands for finite element and nfe is the number of finite elements.To write the variation of kinetic energy one applies the principle of conservation of mass, i.e.: Remembering that, for trusses, mass is considered concentrated at nodes, one writes: in which, M  is the mass corresponding to node  , calculated as half the sum of the masses corresponding to finite elements connected to the node plus, when needed, the concentrated masses existent in the structure.

Strain Energy
In the previous item the mechanical energies that depend directly on the position of the nodes of the structure are presented.This section shows the strain energy that depends indirectly on the nodal positions.Firstly, one defines the adopted strain measure and its dependence of positions.
Then the specific strain energy is presented as a function of strain and the dependence of nodal position becomes clear.
Figure 1 shows a truss element before and after the change of position or movement.The initial position of the element is defined by the initial position of nodes, i.e., i X  in which 1, 2,3 i  represents direction and  nodes.The current position is defined by the unknown variables, i.e., the current nodal positions called i Y  .First of all it is important to inform that unknown nodal positions are determined by a predictor corrector numerical algorithm.Therefore, the following expressions are written containing these unknowns, but trial values are always known.The initial and current lengths 0 L and L are calculated by: The one-dimensional Green-Lagrange strain is chosen to describe the material behavior, i.e.: From the calculated Green's strain one can choose any Lagrangian constitutive relation to represent the material.In this work the uniaxial constitutive relation of Saint-Venant-Kirchhoff is chosen and is expressed by the specific strain energy and its derivative with respect to the Green strain as in which E is the elastic modulus that coincides with the Young modulus for infinitesimal (or small) strains and S is the longitudinal second Piola-Kirchhoff stress related to the Cauchy stress by as described by Ogden (1984) and Bonet et al. (2000).
Integrating the specific strain energy on the initial volume of the bar (Lagrangian quantities) one finds the accumulated strain energy on a finite element and structure, such as: It is worth noting that all quantities are assumed to be constant along the bar, facilitating the integration process.It is obvious that expression (12) depends on the current positions of the nodes of the structure.Then, one calculates the variation of strain energy regarding positions as: Using expression (12) one writes, and considering ( 9) and (10) one achieves: Latin American Journal of Solids and Structures 13 (2016) 964-991 Substituting equation ( 16) into ( 14) and using the energy conjugate definition (Ogden (1984) and Bonet et al. (2000)), that is, the derivative of strain energy regarding position (or displacement) results internal force, one writes: ( 1) Equation ( 17) represents the contribution of an element for the total internal force of the analyzed structure.The internal forces are calculated for all finite elements and are combined into a single vector of internal forces containing all nodes of the structure by simply respecting their numbering.
Considering all finite elements the variation of strain energy, equation ( 13), results: where the index fe is omitted as the sum that corresponds to all nodes of the structure has already been performed.It is of interest, as will be shown next, to calculate the hessian of the strain energy potential as: For one finite element one has:

Damping -Modified Kelvin/Voigt Model
Despite the constitutive model presented here is dedicated to the device mass / spring / damper it can be used in any finite element.The Kelvin constitutive model indicates that the dissipation of power is generated by a force proportional to the speed of strain.As the Green strain is adopted, Latin American Journal of Solids and Structures 13 (2016) 964-991 the proposed model is called as Modified Kelvin Model.Before describing the proposed model it is important to note that, in general, it is not possible to find an explicit expression for the dissipative potential Q presented in equation ( 1), however, it is possible to write its internal force and, therefore, the variation of dissipation Q  .A schematic representation of the proposed model can be seen in Figure 2.
, From Figure 2 the Piola-Kirchhof stress is given by: for which the first term is the derivative of the specific strain energy in relation to the Green strain, see equation ( 11).So the second term corresponds to the dissipative part of the potential energy.
For a finite element one writes the variation of the dissipative potential regarding nodal positions as: where q is the dissipative potential by unit of volume (not written) and all variables are considered constant for a truss element.Using equation ( 16) one defines the contribution of one element to the global dissipative force, as: Considering all elements the variation of the dissipative potential results: where the index fe is omitted as the expression includes all nodes of the analyzed structure.
Latin American Journal of Solids and Structures 13 (2016) 964-991 To complete the consideration of viscosity, it should be mentioned that the solution proposed here is pure numerical and hence the speed of strain is calculated by the finite difference method in solution of nonlinear equilibrium equation.This strategy is originally described here for nonlinear strain measure, but was inspired in Mesquita and Coda (2002), Mesquita and Coda (2003) and Mesquita and Coda (2007) where linear viscoelastic models are proposed.As far as the authors knowledge goes it is different from all nonlinear viscoelastic formulations present in literature, see for example, Lemaitre and Chaboche (1990), Lemaitre (2001), Lima et al. (2015), Holzapfel (1996), Simo (1987), Huber and Tsakmakis (2000), Petiteau et al (2013), Marko et al. (2006) and Clough and Penzien (1985).
To introduce the numerical approximation one assumes a time step t  and writes: Substituting equation ( 26) into (24) results: Considering the time step sufficiently small one assumes: Comparing equations ( 27) and ( 17),one writes (27) as: Therefore, the viscous internal force is numerically related to the elastic internal force.

NONLINEAR EQUILIBRIUM EQUATION AND SOLUTION PROCEDURE
Before introducing the Lagrange multipliers it is of interest to present the structure solution considering the dissipative term concluding the description of main original contribution of this work.Substituting equations ( 4), ( 7), ( 18), (25) into equation ( 2) one writes: As the position variation i Y   is arbitrary, equation ( 30) can be rewritten as: Or making the correspondence between the node number  and the movement direction i with a generic degree of freedom that is the nonlinear equation of dynamic equilibrium, or movement, of the analyzed structure.The adopted solution process combines: (i) the time integration of the dynamic equilibrium by the Newmark method, (ii) the original viscous force approximation by equation ( 29) and (iii) the Newton-Raphson nonlinear solver for equation (32).
Equation ( 32) is rewritten for a current instant ( ) where the mass proportional damping is also introduced, following Coda et al. ( 2013), Coda and Paccola (2014) and Coda and Paccola (2008), in order to not loss the generality.The Newmark approximations are used in the following form, where t  is a time step.
Substituting equations ( 35) and (36) in equation ( 34) one finds:  .Its solution is done by the Newton-Raphson procedure from a Taylor expansion truncated at first order, as: where In equation ( 41) estat H is the static hessian given, for each finite element, by equation ( 19).In this equation the approximation described by equation ( 29) is employed, which implies that only the current viscous force depends on the current position and has influence in the dynamic hessian, see equation ( 37).
From equation ( 40) results the linear system of equation used to calculate a correction for the trial current position, i.e.,

 
  Equation ( 44) is used to correct velocity in equation ( 36) and the iteration continues until, in which TOL is a pre defined tolerance.It is worth noting that t Q , t R and int t F  are updated only at the end of the iterative procedure, i.e., at the beginning of the next time step.For the first time step the acceleration is calculated by: Next section shows how the solution procedure is completed considering the Lagrange multiplier to include vibration control devices.

INTRODUCTION OF THE LAGRANGE MULTIPLIER
Despite the damping described above can be used in any bar of the structure,in this section it is specially adapted to simulate the passive mass / spring / damper system of vibration control shown in Figure 3.This figure shows a mass / spring /damper system named z for which the stiffness is 2K , the damping parameter is 2 and the mass is m .
As the finite element simple bar (truss) cannot resist transverse loads it is necessary to constrain the device mass movement to a straight line between nodes k and r of the system.The two equations to impose the necessary constraint for a three-dimensional device are: and The introduction of these constraints on the equilibrium equations is made by adding in the total energy, equation (1), the sum of the following Lagrangian potentials where z varies from 1 to the number of installed devices nz , i.e.: Latin American Journal of Solids and Structures 13 (2016) 964-991 The search of the stationarity is now made including the above restriction, as: Remembering that the Lagrange multipliers are also free variables, it is written where forces ˆi F  are actives when nodes  belong to the devices z and z j  are the associated La- grange multipliers.Values of ˆi F  and z j  are given at appendix.
Using index notation over equation ( 51) and taking into account equation ( 52), equation ( 30) is rewritten as: Equation ( 54) is very similar to equation (31), as the active coordinates of ˆi F  coincide with the coordinates of (int) i F  .Equation (55) coincides with the imposed constraints, see equations ( 47) and (48).Thus, the additional number of equations is   where d is the dimension of the con- sidered space (2D or 3D) and nz is the number of devices.
The solution process is identical to that described in the previous section where both the force vector (equations (A1) to (A11)) and the Hessian matrix will be changed.In the three-dimensional case, for a generic device z , one writes the terms to be added into the hessian matrix, as: or, developing the involved derivatives, one finds: Latin American Journal of Solids and Structures 13 (2016) 964-991 where the last two rows and columns correspond to the additional degrees of freedom per device.
The other values are added to existent positions of the original hessian matrix.It is worth noting that in the solution process vectors Y  and Y   also contain the updating of Lagrange multipliers, but these are not considered to calculate the acceptable error / Y X    .

NUMERICAL SIMULATIONS
Some examples are shown to confirm the various aspects of the formulation.They are organized from very simple cases to more complicated ones.

Mass / Spring / Damper Device Validation
The example in figure 4  For this level of strain it is also expected that the quasi-static problem, including Kelvin/Voigt viscosity and discarding inertial forces, approaches the linear problem.Thus, the horizontal displacement as a function of time is given by this is the viscous damping equation however, for this simple problem, solutions are similar.Adopting 4MPa s    , figure 5 compares the achieved solutions for various time steps, revealing the convergence of the proposed numerical viscoelastic approach to the analytical solution.When the load is small the numerical response approaches better the linear response, for larger loads it was expected very different results, however, it is a good surprise that the damping behavior continues presenting the same pattern.For the purposes of this study the results are more than satisfactory, moreover the small difference indicates that it is possible to search better constitutive models to better fit the real device to be used.

Optimization of a Mass / Spring / Damper Device
In this and the next example it is used a hypothetical plane structure for representing a building with 12 floors, 36m in height and 5m in base.Its is used 49 truss finite elements with density , elastic modulus of 210GPa 

E
and cross sectional area of The structure is constrained at the base by two fixed supports and the vibration control device is fixed at the top floor as indicated in figure 6. Figure 7   The natural frequencies of the structure are calculated without and with the mass/spring device (not damped) by equation  where 0 H is the hessian matrix for the initial position that coincides with the linear stiffness matrix of the structure.The first frequency without the mass/spring device is 3.1253Hz when the device is coupled the first frequency is 0.963 .Hz Applying sinusoidal horizontal base movement in the same frequency of structures (with and without the mass/spring device) and with amplitude of 5cm , as both systems are not damped, the resonance phenomenon takes place.The structure without the device reaches a top displacement amplitude of 6m in less than 5s , while the structure with the mass/spring (not damped) device reaches for the same period of time an amplitude less than 2 .
m Thus, the introduction of the mass/spring system without damping would have vibration control effect only for external actions of short duration.Now it is introduced the modified Kelvin/Voigt damping    in the vibration control device as shown in equation ( 41).In order to optimize the device one varies  from .From previous analysis, it was found that with the damping variation there was a small variation in the natural frequency of the structure.Due to this reason it is necessary to vary the excitement frequencies from 0.963Hz through 1.04Hz with a step of in order to find the worse excitation frequency for each damping parameter  .Figure 8 shows the maximum displacement of the top of the structure as a function of the damping parameter  for the worse possible excitation frequency.Figure 9 shows the worse frequency for each damping parameter  . . Figure 10 shows the displacement of the top of the structure as a function of time for these values of damping and frequency.As can be seen, in addition to smaller amplitude, the presence of damping eliminates, as expected, the resonance profile.If a longer time is presented the maximum displacement stabilizes at 0, 7m as shown in figure 8.The proposed device works very well and its use is not limited to small displacements as the formulation presents a geometrically exact description.One should note that the mass of the device developed horizontal and vertical movements keeping the alignment with nodes 1 and 3 of figure 7.

Structure under Earthquake-2D Model
In this example, the behavior of the previous studied structure is analyzed when it is submitted to the effect of a real earthquake, see for instance Reis and Coda (2014), Coda and Paccola (2014) and <http://peer.berkeley.edu/peer_ground_motion_database/spectras/21713/unscaled_searches/84357/edit.>.The imposed base movements in horizontal and vertical directions are depicted in figure 11.
In addition to the mass / spring / damper device, the sliding base device is created by changing the configuration of the first floor, see Figure 12.As in the previous example 50 truss elements with density   An optimization process is used to find the elastic modulus ( 0.5GPa  E ) of the cross diagonals of the first floor limiting the base angle  to be less than 0 15 (see figure 12).The same damping constant floor and the horizontal acceleration of the top floor are evaluated.To perform this analysis four situations are considered: (i) without any vibration control device (ii) with vibration control mass / spring / damper (iii) with the sliding base device Figure 12 and (iv) with the combination of the two devices.
Figures 13 to 21 show results along time for the above described cases.As can be seen the sliding base device allows larger displacement amplitude, but reduces the frequency of the horizontal movement of the top floor.Regarding the acceleration of top floor, the device mass / spring / damper is the best.Regarding the behavior of the normal force on the right column of the second floor, both the sliding base device and the mass / spring / damper present equivalent performances.Now the combination of devices, case (iv), is analyzed and results are shown in figures 22, 23 and 24.As shown in Figures 22, 23 and 24 the structure with two devices presents intermediate displacement in relation to the individual devices, moreover one can observe the beating profile.Furthermore, it was found that in this case the oscillations present larger periods.On the other hand, it should be noted that in terms of acceleration, this combination presents values lower than the massspring-damper system, thus offering more comfort to users.The device combination present an intermediate maximum normal force amplitude.Thus, the combination of devices offers a better comfort for users and improves the strength of the structure when compared to the sole mass / spring / damper device.Moreover, this example shows that the proposed viscoelastic rheological model is successfully implemented, is very stable and useful for general analysis.

3D Water Thank Tower
To confirm the possibilities and usefulness of the proposed formulation in practical analysis, in this example a 3D application is carried out with the intention of giving obvious and consistent results.The example consists of the analysis of a light structure that supports a water tank subjected to a horizontal impact with the composition of forces Due to the Lagrange multiplier technique, the water tank is constrained to move along the straight line between nodes 1 and 3. Figure 28 shows some selected positions of the tower movement.The proposed examples reveal that the mass / spring / damping device is working properly for 2D and 3D applications.Example 1 shows that the rheological model and the Lagrange Multiplier technique are appropriate to static and dynamic application and unify approaches.Examples 2, 3 and 4 show that the formulation is able to model large displacement situations using large strain mass / spring / damping devices.Furthermore, the success of numerically integrating the Kelvin-Voigt rheological model opens the possibility of developing more realistic models for nonlinear dynamic damping.

CONCLUSIONS
This study presents a total Lagrangian formulation to perform nonlinear dynamic analysis of plane and space trusses.For this total Lagrangian formulation it is proposed an original Lagrangian multiplier and a modified Kelvin viscoelastic constitutive model approach to consider sliding mass / spring / damper devices for vibration control.The more important contribution is the proposed modified Kelvin viscoelastic constitutive model that make possible assembling damping for static and dynamic analysis from phenomenological viscosity at large strain using the Green strain measure.The main advantage of this formulation, when compared to the classical ones, is the possibility of its extension to comprise more complex viscous behavior without the necessity of analytically solving differential equations.An interesting surprise of technical interest is the same damping behavior of the modified Kelvin model for small and large displacements and strains.The success of numerically integrating the Kelvin-Voigt rheological model opens the possibility of developing more realistic models for nonlinear dynamic damping and its extension to 3D viscoplastic analysis.

Appendix
For each device z one calculates the conjugate forces, ˆi F  and z j  , as: When a 2D truss is considered one disregards the components 3 F  , 2 z  and the terms 3 Y  .

Figure 5 :
Figure 5: Quasi-static results for the device.
shows a detail of the device and the fixing strategy.The supporting bars have the same properties of other structural bars, but bars 1-2 and 2-3 of the device have elastic modulus 1GPa 

Figure 8 :Figure 9 :
Figure 8: Maximum displacement amplitude at the top as a function of  .

Figure 10 :
Figure 10: Top displacement as a function of time for the optimum damping.

Figure 11 :
Figure 11: Base movements in horizontal and vertical directions.
diagonal bars of the first floor.In this example the response of the structure to the excitation presented in Figure 11 is studied.The horizontal displacement of the top floor, the normal force on the right column of the second Latin American Journal of Solids and Structures 13 (2016) 964-991

Table 1 :
Maximum amplitude by the static displacement.