Time Domain Modeling and Simulation of Nonlinear Slender Viscoelastic Beams Associating Cosserat Theory and a Fractional Derivative Model

A broad class of engineering systems can be satisfactory modeled under the assumptions of small deformations and linear material properties. However, many mechanical systems used in modern applications, like structural elements typical of aerospace and petroleum industries, have been characterized by increased slen-derness and high static and dynamic loads. In such situations, it becomes indispensable to consider the nonlinear geometric effects and/or material nonlinear behavior. At the same time, in many cases involving dynamic loads, there comes the need for attenuation of vibration levels. In this context, this paper describes the development and validation of numerical models of viscoelastic slender beam-like structures undergoing large displacements. The numerical approach is based on the combination of the nonlinear Cosserat beam theory and a viscoelastic model based on Fractional Derivatives. Such combination enables to derive nonlinear equations of motion that, upon finite element discretization, can be used for predicting the dynamic behavior of the structure in the time domain, accounting for geometric nonlinearity and viscoelastic damping. The modeling methodology is illustrated and validated by numerical simulations, the results of which are compared to others available in the literature.


INTRODUCTION
Highly flexible structural components such as slender beams and cables can be frequently found in several modern engineering systems, such as electrical energy transmission lines, cable-stayed towers and bridges, and subsea oil risers.Most often, slender structural members are distinguished by large aspect ratios, defined as the ratio between the length and a characteristic dimension of crosssections, and are prone to undergo large displacements and rotations, under static or dynamic loads.These features have led to the need for developing adequate modeling procedures accounting for geometric nonlinearities.Moreover, in the cases where the materials do no present Hookean behavior, physical nonlinearity must also be accounted for.
The so-called Cosserat theory, intended for the modeling of slender unidimensional continua, has attracted the interest of many researchers lately.This theory was originally developed by the French brothers Eugene and François Cosserat in the beginning of the 20th century, and remained overlooked for many years, most probably due to the insufficiency of computational resources required for its practical implementation.More recently, the increasing need for the accurate modeling of nonlinear structural systems of practical interest, combined with the advances achieved in computer processing capacity, has stimulated the rediscovery of Cosserat Theory, allowing its application to numerous problems in different fields of knowledge.
After the seminal work of the Cosserat brothers in 1909, some of the more prominent contributions to the development and application of their theory should be mentioned.Antman (1995aAntman ( , 1995b) ) presents an application of the Cosserat theory to describe the behavior of rods.This author states that the Cosserat theory presents several advantages with respect to others used in structural mechanics, as it is not based on specific geometrical approximations or mechanical assumptions.The Cosserat theory has been extensively employed for the modeling of different kinds of physical systems, including vertical columns of oil wells (Tucker and Wang, 1999), cables for surgical applications (Pai, 2002), Micro Electro Mechanical Systems (MEMS) (Liu et al., 2004), DNA chains (Maddocks, 2004), blood flow in veins (Carapau, 2006).More recently, it was used to describe the three-dimensional dynamic behavior of oil drills (Alamo and Weber, 2006) and subsea oil risers accounting for fluid-structure interactions (Borges et al., 2010).
In most of the previous studies, energy dissipation (damping) has been disregarded or considered in excessively simplified manners.On the other hand, the accurate estimation of structural damping is of primary importance in numerous situations, such as in the design against fatigue failure of offshore umbilical cables or risers subjected to Vortex Induced Vibrations (VIV).
Damping can be associated with the inherent characteristics of the structural systems, such as material hysteresis and friction between solid parts.Also, when needed, damping can be introduced or increased voluntarily by employing dissipative materials, among which viscoelastic materials have been frequently used, given their operational and cost effectiveness.Moreover, the maturity achieved in the use of those materials has led to their availability as commercial products, which can be purchased "off the shelf" or customized according to specific needs (Rao, 2003).
The viscoelastic behavior, presented by various types of solid materials, such as glass and polymers, can be regarded as the combination of the elastic behavior exhibited by Hookean solids and the behavior of Newtonian (viscous) fluids (Nashif et al., 1985;Christensen, 1982).As a result, the mechanical response of those materials to external loads is not instantaneous and depends on the Latin American Journal of Solids and Structures 14 (2017) 153-173 load history, which means that those materials exhibit memory effect.Additionally, their response is lagged behind under the action of cyclic loads, leading to a stress-strain hysteresis loops, the dimensions of which determine the energy dissipated in each cycle.This behavior is the rationale for using viscoelastic materials as a means of increasing damping in structures, which can be done either by means of discrete devices (mounts and joints), or surface treatments (free or constrained layers) (Nashif et al., 1985;Mead, 1998).
Although the use of viscoelastic materials in engineering systems has achieved considerable maturity, there are still a number of issues involved in the modeling of the dynamic behavior of complex viscoelastic structures.Indeed, it is widely known that the dynamic behavior of viscoelastic materials depends heavily on a number of environmental and operational factors, mainly the vibration frequency and temperature (Christensen, 1982).This dependency can be accounted for by employing the concept of complex modulus in combination with the Frequency-Temperature Superposition Principle (FTSP).The former, intended for the characterization of the viscoelastic behavior in the frequency domain, consists in ascribing complex-valued, frequency-dependent moduli to the viscoelastic material, while the FTSP establishes the equivalence between the effects of the excitation frequency and the temperature on the material moduli, for the so-called thermoreologically simple materials.
The complex modulus approach, combined with the FTSP, enables the modeling of the dynamic behavior of rather complex structures containing viscoelastic materials in the frequency domain, provided that the dependency of the moduli on frequency and temperature is adequately handled within efficient numerical procedures, such as finite element discretization.Frequency domain response predictions is quite straightforward, although much research effort has been devoted to important aspects such as model reduction procedures aiming at reducing the computation cost involved in the resolution of large scale problems, and the modeling of uncertainty propagation.Such an approach has been addressed in a number of previous studies (De Lima et al., 2010a;De Lima et al., 2010b;Johnson et al., 1981;Cazenove et al., 2012;Moita et al., 2013).
However, the numerical prediction of time domain responses of viscoelastic systems can be preferable or indispensable in a number of situations, such as problems involving transient loads or nonlinear behavior.As compared to frequency domain modeling, the prediction of time domain responses involves higher computational costs, since the computation of time domain responses involves the resolution of convolution integrals for the appropriate consideration of the memory effect of viscoelastic materials.
Various strategies have been proposed for the calculation of time domain responses of viscoelastic structures in association with finite element discretization.The Anelastic Displacement Fields (ADF) model, proposed by Lesieutre and Bianchini (1995), separates the mechanical displacement field into elastic and anelastic parts.The latter is used to describe the part of the strain that is not instantaneously proportional to the stress, and is expressed as a series of anelastic fields, whereby the time evolution of each of them is governed by a first-order ordinary differential equation.The resulting equations of motion are expressed in the same form as the equations of motion of linear, viscously damped discrete systems, comprising ( ) where n is the number of structural degrees-of-freedom and A n is the number of anelastic fields adopted.It has been recog- nized that the major shortcomings of the ADF model is the increase of the number of degrees-offreedom and the necessity of separating structural vibration modes from those associated to the internal variables.
Another approach is the Golla-Hughes-McTavish (GHM) model (Golla and Hughes, 1985).According to this model, the complex modulus of the viscoelastic material is expanded in the Laplace domain as the sum of G n rational transfer functions.The time-domain equations of motion can be derived after using the inverse Laplace transform, where internal dissipative coordinates can be defined and related to their elastic counterparts.The GHM model generates a set of ( ) ear second-order differential equations of motion, and their resolution involve the same difficulties as those already described for the ADF model.The use of viscoelastic models based on the notion of fractional derivatives, pioneered by Bagley andTorvik (1983, 1985), is another promising modeling strategy.In a first phase, these models were used to predict frequency domain responses.Subsequently, time domain computations associated with finite element models began receiving greater attention.Schmidt and Gaul (2001) developed a three-dimensional finite element model incorporating a five-parameter fractional derivative constitutive equation, solved by means of the Grünwald-Letnikov discretization approach.Schmidt and Gaul (2006) improved the efficiency of this model through a novel discretization scheme for the fractional derivative operator.Galucio et al. (2004) implemented a viscoelastic sandwich beam finite element based on a four-parameter fractional constitutive model solved by the Grünwald-Letnikov discretization scheme.
The association of viscoelastic models based on fractional derivatives and finite element discretization is considered to be an efficient strategy for the simulation of time domains responses of viscoelastic systems.Certainly, it does not avoid the computation of convolution integrals but, making use of appropriate time integration schemes, can render computations satisfactorily affordable and accurate.
The study reported in this paper is motivated by the existence of a number of practical problems involving large amplitude vibrations of slender one-dimensional structures, to which viscoelastic behavior can be introduced, either through the basic constituent materials or devices applied to the original structure, aiming at increasing damping levels and, as a consequence, decreasing vibration amplitudes.
In the next sections, the foundations of the Cosserat theory is reviewed and the incorporation of the fractional derivative viscoelastic model into the discretized nonlinear equations of motion is described.Numerical simulations are described for the purpose of demonstrating the main features of the modeling strategy and evaluating the effectiveness of the viscoelastic behavior in attenuating vibrations of the particular types of structures considered herein.  .However, when shear is accounted for, this coincidence no longer occurs.However, it is widely accepted that, for slender beams, the shear effect can be neglected, so that the cross-sections remain perpendicular to the tangents to the centroid line (Wang et al., 2004).This assumption is adopted herein.

Referring to
The motion of the centroid line ( ) , s t r can be described in the inertial basis as: The characterization of the motion involves both the velocity of the centroids ( ) andthe angular velocity of the cross-sections, ( ) , s t w .First of all, it is necessary to establish the variation of the local vector basis with respect to time, which, from pure kinematic considerations, is expressed as follows: Latin American Journal of Solids and Structures 14 (2017) 153-173   ( According to the Cosserat theory, the strains are classified in two groups: linear strains, ( ) , s t v , and angular strains, ( ) From Eq. ( 2), it follows that: Similarly, from Eq. ( 4) one gets: As stated previously, since the effect of shear is neglected, each cross-sections of the beam remains perpendicular to the direction defined by the tangent to the centroid curve, which enables to write: In this case, one has: where: Latin American Journal of Solids and Structures 14 (2017) 153-173   ( , , The cross-section rotations must be related to the linear strains in order to completely describe the kinematics of the Cosserat beam, and this can be done by using rotation parameters such as Euler angles, rotational vector (Euler vector), quaternion parameters and Cayley transformation (Dongsheng et al., 2007;Cao et al., 2008;Rubin and Brand, 2007).
In the work reported here, Euler vector and Euler angles are jointly used.First, for the parameterization based on Euler vector, the torsion angle of the cross-section, ( ) , s t j , and the linear strain vector, ( ) , are used as parameters.Accordingly, the following relations relating the two vector bases shown in Fig. 1 are obtained: ( ) ) v with the aid of Eq. (11.c), the following relations are obtained (for simplifi- cation, time and space dependences are omitted): (12.a)After obtaining the relationship between the moving basis and the inertial basis using Euler vector, in order to relate the linear strains with angles of rotation of the cross-section, one makes use of Euler angles.Accordingly, the orientation of the moving axes is obtained by performing three successive rotations of the inertial axes.The angles of rotation about the axes 1 2 3 , and e e e are defined as , Expanding the trigonometric relations of the transformation matrix obtained through the Euler angles in third order Taylor series, and equating the resulting expression to Eqs. (12.a)-(12.c), the following relations are obtained (Cao et al., 2006):

Equations of Motion
According to Antman (1995a), from the application of the Newton-Euler equations, the local dynamic behavior of a differential element of Cosserat beam with mass density ( ) s r and cross-section are ( ) A s is described by the following partial differential equations: where ( ) In order to obtain the equations of motion in terms of displacements ( ) Thus, the derivative of the vector of internal forces with respect to s is developed as follows: ( ) Similarly, the derivatives of the vector of internal moments and angular momenta are developed, respectively, as: ( ) Moreover, the second term on the right side of Equation (14.a) is developed as follows: ( ) ( ) ( ) It is assumed that the density of the material, ( ) s r , as well as the cross-sectional area, ( ) A s , depend only the spatial variable s, and that the center of mass of the differential element coincides with the cross-section centroid.Another hypothesis adopted is that the beam can undergo large displacements while maintaining small deformations.Accordingly, a constitutive model that characterizes homogeneous, isotropic materials, based on Kirchhoff constitutive relations is adopted.
Latin American Journal of Solids and Structures 14 (2017) 153-173 Based on these hypotheses, the following expressions are obtained for the elements of the main diagonals of matrices ( ) s K and ( ) s J , representing the mechanical stiffness of the structure: where ) are the principal second moments of area of cross-section, 11 J and 22 J represent the flexural rigidity, 33 J represents the torsional rigidity.The components 11 K and 22 K represent shear stiffness while 33 K is the axial stiffness.
Using Kirchhoff's constitutive relations, the contact forces and moments are expressed in terms of linear and angular deformations, respectively, as follows (Cao et al., 2006):

Finite Element Discretization
Lagrange equations are used to derive the differential equations of motion for a general Cosserat beam finite element.For this purpose, the kinetic and potential energy, as well as the virtual work of the nonconservative forces must be formulated, after performing interpolation of the displacement field within the element.As suggested by Cao et al. (2008), for this interpolation, non-linear shape functions obtained from the resolution of the equilibrium equations corresponding to the associated static problem are used.
This problem is represented by the following particularization of Eqs. ( 14 The ordinary differential equations of static equilibrium are solved by using a perturbation technique.While the literature offers several methods, the Frobenius method (Arfken et al., 2000) was selected in this paper.
Certainly, the shape functions can be arbitrarily chosen.However, aiming at achieving the numerical accuracy required for representing a system undergoing large displacements, at affordable computational cost, in this paper, the shape functions were approximated using terms up to the third order, as follows: More details about the derivation of the shape functions using the method of perturbation can be found in Cao et al. (2008) and Wang et al. (2004).

Finite Element Dynamic Analysis
The kinetic energy per unit length of the Cosserat beam finite element is expressed as: where ( ) s I is the mass inertia matrix, with respect to a set of barycentric orthogonal axes, and the strain energy per unit length can be expressed in terms of the vectors ( ) where ( ) s K and ( ) s J are defined in Eqs. ( 19) and ( 20).
At the element of a single finite element, the Lagrangian function can thus be expressed in terms of a set of generalized coordinates and velocities associated to the nodal degrees of freedom, ( ) ( ) e q  , as follows (from now on, superscripts (e) are used to indicate quantities defined at element level): while the expression of the virtual work of non-conservative forces is: It should be further clarified that the forces acting on the element are composed of three additive parts: the first comes from the interaction between neighboring elements, ( ) i e f ; the second is due to the action of external concentrated forces applied at the nodes, ( ) c e f ; and the third is associ- ated to the loads forces distributed along the element, ( ) f .As a result, the total virtual work done by the nonconservative forces is given by: Substituting Eqs. ( 28), ( 29) and (30) into the Lagrange equations, the equations of motion at element level are obtained as follows.It should be mentioned that, owing to the complexity of the algebraic manipulations involved, they have been conducted with the aid of specialized computer programs intended for symbolic manipulation: In the equation above, ( ) e M is the matrix of inertia, ( ) e K is the linear stiffness matrix, and g q is the nonlinear vector containing the cubic and quadratic terms in ( ) Once the matrices of mass, stiffness and nodal equivalent forces for each Cosserat beam element are determined, the mechanical connection among different elements must be enforced.This is made through the standard procedure of assembling element-level equations of motion into the Latin American Journal of Solids and Structures 14 (2017) 153-173 global equations of motion for the system.Also, geometrical boundary conditions must be enforced by a proper modification of the global equations of motion (Bathe, 1996).

VISCOELASTIC FRACTIONAL DERIVATIVE MODEL
In the study reported herein, a model based on fractional derivatives was adopted for modeling of the viscoelasticity behavior in the time domain, in association with finite element discretization, according to the procedure suggested by Galucio et al. (2004).
According to the model adopted, the one-dimensional viscoelastic constitutive equation is expressed as follows: where ( ) t s is the stress, ( ) t e is the strain, t is the relaxation time, and 0 E and E are, respec- tively, the low-frequency and high-frequency elastic modulus.Also, in Eq. ( 32), t D a é ù ⋅ ê ú ë û indicates the Riemann-Liouville fractional derivative operator, defined as: where 0 1 a < < is the fraction derivative order and ( ) G ⋅ designates the Gamma function.
The fractional derivative operator is approximated by using the Grünwald-Letnikov backwarddifferentiation approximation, as follows: ( ) where t D is the constant-time step used for discretization, t N is the number of time intervals such that t t N t = D , and 1 j A + are the so-called Grünwald coefficients, which are expressed as: and satisfy the following recurrence formula: Latin American Journal of Solids and Structures 14 (2017) 153-173 which is interpreted as an anelastic strain, using the Grünwald-Letnikov approximation given in Eq. ( 35), and recalling that 1 1, A = one obtains: ( ) where c t  32) and ( 37), the following discretized form of the constitutive equation is obtained: The discrete-time constitutive equation applicable to multiaxial stress/strain states can be obtained by replacing the scalar stress, strain and material moduli by the respective equivalent vectors and matrices, as follows: In the derivation of Eq. ( 39), it has been assumed that: ( ) ( ) in such a way that the material coefficient matrices can be expressed as: where matrix  E depends only on the Poisson ratio .
One way to include the viscoelastic behavior into the finite element model is by formulating the virtual work done by the non-conservative internal forces developed within a viscoelastic material.For this purpose, one resumes the development presented in the previous section, writing:  and the virtual work of the internal forces is obtained as follows: where: e e e e e e t j t t j t t j t and, Noting that the inclusion of the viscoelastic behavior does not affect the kinetic energy, upon using the Lagrange Equations, the following nonlinear differential equations of motion for a Cosserat viscoelastic beam at elementary level are obtained: Regarding the numerical resolution of the global equations of motion, a composite, implicit time integration scheme, previously proposed by Bathe and Baig (2005), has been chosen, since it has been found to provide accurate results with a constant time step.This latter is a constraint imposed by the Grünwald-Letnikov scheme adopted for the discretization fractional derivatives.Therefore, the numerical algorithm implemented consists in computing an intermediate response between time instants t and t t + D based on the trapezoidal rule, and using backward-differentiation formulae to evaluate the response at t t + D , from the response calculated at instant t and the intermediate response.Clearly, since the method is implicit, a method intended for solving nonlinear algebraic equations, such as Newton-Raphson, has also been used.
The flowchart presented in Fig. 2 illustrates the main operations involved in the computation of time domain responses associating Cosserat beam theory and fractional derivative model.

Axial Loading
For the axial loading case, the dimensions of the beam illustrated in Fig. 3(a), are as follows: L=500 mm, b=50 mm, h=50 mm.This example was considered by Galucio et al. (2004).Although it does not involve transverse motion, it is considered to be useful for validation of the numerical implementation of the Cosserat theory in association with the viscoelastic model, by comparison with results available in the literature.The properties adopted for the viscoelastic material are: r =1000 kg/m 3 , 0 E = 7 MPa, E ¥ = 10 MPa, t = ms and a = 0.5.The external load is applied at the free-end of the clamped beam, acting along the positive direction of the longitudinal x axis, as illustrated in Fig. 3(a).According to Fig. 3(c) this load is given by 0 Spatial discretization is done with 10 finite elements, leading to a mesh with 11 nodes and 66 degrees of freedom before the of the boundary conditions.Total simulation time has been chosen as 400 ms, and integration was performed with a time increment of 1 ms.
Figure 4 shows the normalized axial displacement of the tip of the viscoelastic beam, calculated through the Cosserat/fractional derivative modeling approach, as compared to other modeling strategies.These latter are based on a three-layer sandwich beam finite element, containing a viscoelastic core between two elastic faces, developed by Galucio et al. (2004).This model degenerates to a Timoshenko beam element when the outer faces are neglected.As the responses provided by both models are very close to each other, it can be concluded that the viscoelastic damping effect is correctly modeled in association with the Cosserat beam theory.Moreover, as only small displacements take place, the linear and nonlinear models predict the virtually the same dynamic responses.To verify that this is not the case when the structure undergoes large displacements, transverse loading of higher magnitude, applied to a slender viscoelastic beam, is considered next.

F t F H t =
, is applied at the free-end of the clamped beam, as shown in Fig. 3(b), and takes, successfully, the following values of amplitude: 0.01 N, 0.5 N, 1 N, 2.5 N and 5 N.This pattern has been chosen to demonstrate increasing nonlinear behavior with increasing load.The discretization mesh employed in this case consists of 20 elements, 21 nodes, and 126 degrees of freedom (before enforcing boundary conditions).A simulation time of 10s was adopted, and numerical integration was performed with a time increment of 2 ms.
In Fig. 5 the results, in terms of transverse displacement at the free end of the beam, are compared to counterparts obtained from a linear viscoelastic beam finite element model based on Timoshenko beam theory, also implemented by the authors.
It can be clearly seen that as the magnitude of the load increases, thus inducing increasing geometric nonlinearities, the response predicted by the linear Timoshenko beam theory and the nonlin-earCosserat beam theory progressively deviate from each other.Differences are observed not only in the vibration frequency and amplitude but also in the damping level.

CONCLUSIONS
This paper proposed a modeling strategy for viscoelastic slender structures consisting in the association of the nonlinear Cosserat beam theory with a fractional differential viscoelastic model.The formulations previously developed for each part of the modeling procedure were combined in such a way to enable the simulation of the dynamic behavior of structures including nonlinear geometric and viscoelastic dissipation phenomena.Numerical simulations, the results of which were compared to others drawn from the literature, indicate that, at least under the considered circumstances, the developed model can accurately predict transient responses.It was also found that the Cosserat-viscoelastic beam finite element model can capture increasing levels of nonlinearity, which affect the amplitude, frequency and damping of the calculated vibration responses.
Fig.1, according to the Cosserat theory, the motion of a beam element is spatially defined in terms of the motion of the curve connecting the centroids of the cross-sections, defined by Latin American Journal of Solids and Structures 14 e e and a set of orthogonal unity vectors attached to the moving cross section, s is the curvilinear coordinate that indicates the position of the cross section along the centroid curve, and t indicates time.

Figure 1 :For
Figure 1:Schematic of a Cosserat beam element.
in Taylor series up to the third order in terms of 1 v , 2 v , j and eliminating the term 3 respectively, the vectors of internal forces (normal and shear), internal torques (flexural and torsional) and angular momentum per unit length; ( ) are related to the vector of angular deformation ( ) s u .The angular momen- tum ( ) , s t h is related to the rotational inertia and the angular velocity vector ( ) . (18), the development of Eqs.(23.a) and (23.b) leads to the following nonlinear differential equations:

Figure 2 :
Figure 2: Procedural flowchart for the computation of time domain responses associating Cosserat beam theory and fractional derivative model.

Figure 3 :
Figure 3: Structures considered in the numerical simulations: (a) clamped-free subjected to axial load; and (b) clamped-free beam subjected to transversal load, (c) Time variation of the external loads applied to cases (a) and (b).

Figure 4 :
Figure 4: Normalized axial displacement of the tip of the viscoelastic beam shown in Fig. 3(a).