Finite element formulation of viscoelastic sandwich beams using fractional derivative operators

This paper presents a ﬁnite element formulation for transient dynamic analysis of sandwich beams with embedded viscoelastic material using fractional derivative constitutive equations. The sandwich conﬁguration is composed of a viscoelastic core (based on Timoshenko theory) sandwiched between elastic faces (based on Euler– Bernoulli assumptions). The viscoelastic model used to describe the behavior of the core is a four-parameter fractional derivative model. Concerning the parameter identiﬁcation, a strategy to estimate the fractional order of the time derivative and the relaxation time is outlined. Curve-ﬁtting aspects are focused, showing a good agreement with experimental data. In order to implement the viscoelastic model into the ﬁnite element formulation, the Gru¨nwald deﬁnition of the fractional operator is employed. To solve the equation of motion, a direct time integration method based on the implicit Newmark scheme is used. One of the particularities of the proposed algorithm lies in the storage of displacement history only, reducing considerably the numerical efforts related to the non-locality of fractional operators. After validations, numerical applications are presented in order to analyze truncation effects (fading memory phenomena) and solution convergence aspects.


Introduction
Many investigations have demonstrated the potential of viscoelastic materials to improve the dynamics of lightly damped structures. There are numerous techniques to incorporate these materials into structures. The constrained layer passive damping treatment is already largely used to reduce structural vibrations, especially in conjunction with active control [2,13]. One of the crucial questions is how to quantify such a material damping if the viscoelastic solid has a weak frequency dependence on its dynamic properties over a broad frequency range. Classical linear viscoelastic models, using integer derivative operators, convolution integral or internal variables, become cumbersome due to the high quantity of parameters needed to describe the material behavior. In order to overcome these difficulties, fractional derivative operators acting on both, strain and stress can be employed.
Until the beginning of the 80s, the concept of fractional derivatives associated to viscoelasticity was regarded as a sort of curve-fitting method. Later, Bagley and Torvik [1] gave a physical justification of this concept in a thermodynamic framework. Their fractional model has become a reference in literature. Special interest is today dedicated to the implementation of fractional constitutive equations into FE formulations. In this context, the numerical methods in the time domain are generally associated with the Grünwald formalism for the fractional order derivative of the stress-strain relation in conjunction with a time discretization scheme. Padovan [7] derived several implicit, explicit and predictor-corrector type algorithms. Escobedo-Torres and Ricles [6] analyzed a numerical procedure based on the central difference method and its stability aspects. The FE formulation proposed by Enelund and Josefson [5] employs a fractional calculus involving convolution integral description with a singular kernel function of Mittag-Leffler type. Most of these approaches were restricted to single-degree-of-freedom systems and bar-type structures. Recently, Schmidt and Gaul [9] presented a 3D finite element implementation of a fractional model which requires the storage of both displacement and stress history.
This work presents a finite element formulation for transient dynamic analysis of sandwich beams with embedded viscoelastic material. The four-parameter fractional derivative model [1] is used to describe the frequency-dependence of the viscoelastic layer. To illustrate this modeling capability, the master curves of a viscoelastic material are presented, showing the agreement between the fractional model and experimental data. Moreover, a strategy to identify the model parameters is discussed. In particular, it is shown that the order of the fractional time derivative is related to the maximum value of the mechanical loss factor, independently of the relaxation time. For applications to structural dynamics, the finite element implementation of the fractional derivative constitutive equation is firstly proposed. Next, the sandwich beam finite element model is presented, assuming Euler-Bernoulli hypotheses for the elastic faces and where r and e are the stress and the strain, E o and E 1 are the relaxed and non-relaxed elastic moduli, and s is the relaxation time. Moreover, the Riemann-Liouville definition of the fractional operator is where C is the gamma function and a the fractional order of the time derivative 0 < a < 1 ð Þ . 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,8].
After calculating the Fourier transform of Eq. (1), one obtains the expression of the elastic complex modulus where r Ã and e Ã are the Fourier transforms of rðtÞ and eðtÞ, respectively. Its behavior in the frequency domain is described between two asymptotic values: the static modulus of elasticity E o ¼ E Ã ðx ! 0Þ and the high-frequency limit value of the dynamic modulus E 1 ¼ E Ã ðx ! 1Þ. The statements 0 < a < 1, s > 0 and E 1 > E o must hold to fulfill the second law of thermodynamics.
It is important to emphasize that the previous complex modulus has successfully been used to fit experimental data for a wide variety of materials.

Parameter identification
In order to identify the four parameters E o , E 1 , s and a of the one-dimensional constitutive equation, various experimental procedures can be used. For isotropic viscoelastic materials, dynamic tests (transient or harmonic tests) allow the measurements of the complex Young modulus (by traction-compression tests) or the complex shear modulus (by torsion tests, for example). Note that the following identification procedure is restricted to the complex Young modulus since it is supposed to be proportional to the complex shear modulus.
From Eq. (3), we can extract its real and imaginary parts, providing the storage modulus and the loss modulus The mechanical loss factor, defined as gðxÞ ¼ E 00 ðxÞ= E 0 ðxÞ, becomes always non-negative due to the hypotheses previously adopted for a, s, E o and E 1 . The measured data used for building the master curves for the 3M ISD112 material at 27 C, in a frequency range between 20 Hz and 5000 Hz, are the real part of the complex shear modulus G 0 and the loss factor g ¼ G 00 =G 0 (see Fig. 1). These material data were supplied by the 3M company and previously employed in [12].
It should be noted that the Poisson ratio is supposed to be frequency-independent. This classical assumption, Assuming the asymptotical values of the storage modulus (E o and E 1 ) and the maximum of the loss factor g max to be known, we can evaluate a as follows which does not depend on s. The estimation of the relaxation time s can be performed by minimization of the gap between theoretical and experimental data of the complex modulus, using, for example, a least squares method.
The procedure described above can be useful for a first estimate of a. Its improvement should be the focus of a future investigation.
For comparison purposes, the master curves for the classical Zener model and the fractional one are illustrated in Fig. 1. For given values of E o and E 1 , we note that the fractional derivative model (represented by the solid line) is much closer to the measured data [12] than the standard solid model (dashed line).

Approximation for the fractional derivatives
The fractional operator d a =dt a , appearing in the constitutive equation (1), 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 t=N can be chosen equal to the time step increment Dt of the numerical scheme. The upper limit of the sum N t is strictly lower than N, and A jþ1 represents the Grünwald coefficients given either in terms of the gamma function or by the recurrence formulae Let us introduce the internal variable as a strain function such that the constitutive equation (1) can be rewritten as This variable change implies that Eq. (9) contains only one fractional derivative term instead of two as in (1). Using the Grünwald approximation (7) with t=N ¼ Dt and noting that A 1 ¼ 1, relation (9) takes the following discretized form where c is a dimensionless constant given by It should be stated that the Grünwald coefficients in Eq. (10), 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.

Kinematic assumptions
Consider a sandwich beam made of two elastic faces and a viscoelastic core. Euler-Bernoulli assumptions are considered for the faces, whereas the core is based on the Timoshenko theory. All layers are assumed perfectly bonded and in plane stress state. The kinematics of the sandwich beam can be clearly interpreted in Fig. 2. The mechanical displacement field within the ith layer can be written as where the subscript i ¼ a; b; c stands for upper, lower 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. Let us introduce the mean and relative axial displacements given by Using Euler-Bernoulli hypotheses for the faces (h f ¼ w 0 for f ¼ a; b with ðÁÞ 0 ¼ @ðÁÞ=@x) and displacement continuity conditions between layers (u xa ¼ u xc at z ¼ h c =2 and u xb ¼ u xc at z ¼ Àh c =2), 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ũ where " h andh are defined by Notice that these displacement fields have already been used in [13].

Strain-displacement relations
From (11) and (12), 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 the membrane strain i and the curvature j i of the ith layer, and the shear strain of the core c c , are defined by Without covering face layers (i.e., h a ¼ h b ¼ 0), the previous generalized strain quantities of the core correspond to those of a single Timoshenko beam.

Sandwich beam element
The finite element formulation of the elastic three layers sandwich beam is described in this section. The generalized displacements d e ¼ ½" u wũ T are discretized with linear (axial displacement) and cubic (deflection) shape functions. They are related to the elementary degreesof-freedom vector q e ¼ ½" u 1 w 1 w 0 1ũ 1 j " u 2 w 2 w 0 2ũ 2 T by d e ¼ Hq e , where the interpolation matrix H is defined as follows Moreover, the strain operators are written in terms of their membrane, bending and shear components where the interpolation vectors related to axial displacements H x , transverse displacements H z and rotations H r , for faces f ¼ aðþÞ; bðÀÞ and core (c), are defined by Thus, for each layer (i ¼ a; b; c), axial and transverse displacements and rotations are and membrane, bending and shear strains can be expressed by Using the previous discretization, the governing equation is written as where F is the external load vector, M i (i ¼ a; b; c) the mass matrix of each layer and K f (f ¼ a; b) the stiffness matrix of the elastic faces and where superscript e stands for elementary quantities.
The element mass matrix of each layer is given by and the element stiffness matrix of the elastic faces (K e f ) is decomposed into a membrane and a bending term The third term on the left-hand side of (22) is defined by Z where B 1c ¼ B mc þ zB bc and B 5c ¼ B sc . The crucial step in this finite element formulation, which is addressed in the following section, is the computation of the integral expression (25) using the fractional derivative constitutive equation. We recall that if the core is elastic, such an expression corresponds to Z where G c and k c are the shear modulus and the shear correction factor, respectively.

Finite element formulation of the viscoelastic core
Taking into account the assumptions made in the beam theory, and using the hypothesis of proportionality between G and E (since m is supposed frequency-independent), we can rearrange Eq. (8) considering the discretized form of the aneslastic strain (10). Therefore, the stress components can be written in terms of axial and shear strains as (for i ¼ 1; 5) where the constants n i can be written in terms of E o (relaxed modulus) and m (Poisson ratio) as Recall that axial and shear strains are defined by (15) and (16). Their discretization is given by e nþ1 ic ¼ B ic q e nþ1 , with B 1c ¼ B mc þ zB bc and B 5c ¼ B sc . Let us extend this discretization to the internal variable " e (i.e., the anelastic strain) such that " e nþ1 ic ¼ B ic " q e nþ1 . Consequently, Eq. (27) becomes By abuse of language, we shall denote the discretized unknowns " q e nþ1 , associated with " e, by ''anelastic displacements''. These unknowns depend on the displacement memory and are updated using the following equation Substituting Eq. (28) into (25), one obtains Z so that the elementary semidiscrete equation of motion can be rewritten as follows where M e ¼ M e a þ M e b þ M e c and K e ¼ K e a þ K e b þ K e c and where the modified stiffness matrix " K e c and loading vector " F e nþ1 , arising from the viscoelastic behavior of the core, are given by " F e nþ1 ¼ Àc It is worthwhile to notice that all the time history dependent terms were shifted to the right-hand side of the governing equation. Remarks: (a) For an elastic material, " K e c and " F e nþ1 ¼ whereq e is the anelastic displacement vector, associated withẽ, and updated using Such a formulation can be derived from the fractional one, taking into account the fact that the only non-zero Grünwald coefficient of the series in (33) is A 2 ¼ À1. It could also be obtained directly from a backward-Euler discretization of the standard solid model. In this case, the constant is Eq. (34), (35), and (36) reduce to those derived from the integration algorithm of viscoelastic solids proposed by Taylor et al. [11]. This formulation, which is unconditionally stable and second-order accurate, requires the evaluation of the convolution integral arising from the solution of the viscoelastic constitutive equation, over a time step ½t n ; t nþ1 , in this particular case d ¼ 1 À expðÀDt=sÞ Dt=s

Algorithm implementation
The Newmark scheme is adopted here due to its versatility for implementation in structural dynamics. Some modifications are carried out in the classical algorithm to obtain a modified scheme (see the scheme below) that is suitable to achieve the transient responses of a sandwich beam with a viscoelastic core in fractional calculus [4]. The Newmark parameters b ¼ 1=4 and c ¼ 1=2 are chosen in order to obtain an unconditionally stable and second order accurate scheme.
1. Enter data a) integration parameters: b, c, Dt b) model parameters: E o , E 1 , s, a, m c) assembled matrices: M, K c , 3. Enter time step loop, assuming that at t n the state is completely known: q n ; _ q n ; € q n ; " q n a) Predict displacement and velocity Evaluate acceleration by solving the following linear system Correct displacement and velocity f) Evaluate the anelastic displacements history 4. Update time step and return to 3. n n þ 1: Recall that for a constant time step, the stiffness matrix is evaluated once. It is important to note that evaluation of the modified loading (step 3b in the algorithm) is relatively simple because only the anelastic displacement history needs to be stored at each time step (see step 3f). We shall observe in the following examples that it is not necessary to store the whole time history but only the most recent one. This approach of truncating the history is accurate and has low numerical costs when dynamical oscillating responses are considered. However, if creep processes are examined, the whole history has to be taken under consideration. Hence, in this case, the truncation of the history will result in noticable errors [10].

Energy conservartion aspects
It is well known that the unconditional stability of the Newmark scheme leads to the conservation of the total energy of the system during the computation of its motion, such that where T, U and U d are the kinetic energy, the stored elastic energy and the left-hand component of the dissipated energy and where the powers P and " P, associated to the external and modified loads, are given by P ¼ _ q T F and " P ¼ _ q T " F. In the following, their corresponding energies are written W and W d , respectively.
Integrating (37) over a time step ½t n ; t nþ1 , we can rewrite the dissipated energies as These two quantities are such that the total dissipated energy is defined by D ¼ U d À W d . It should be worthwhile to note that the approximation for the work of the dissipation forces over one time step (40) results from the integration operator of the Newmark algorithm.

5
Results and analysis

Validation
Two examples taken from literature are chosen in this section in order to validate the fractional derivative algorithm and the sandwich beam model. The first one was extracted from [5] in which the fractional viscoelastic model is based on a convolution integral description with a singular memory kernel of Mittag-Leffler type. In the second example, our solution is compared to the hybrid Laplace transform/finite element method proposed by [3] for a Timoshenko beam with a classical Zener model. In this case, the relaxation modulus is expressed in form of Prony series.

Viscoelastic bar
Consider a viscoelastic fixed-free bar of length L ¼ 500 mm, width b ¼ 50 mm and thickness h ¼ 50 mm, discretized with 10 finite elements. The material data of the viscoelastic material (fictitious polymer) is q ¼ 1000 kg/m 3 , E o ¼ 7 MPa, E 1 ¼ 10 MPa and s ¼ 20 ms. The structure is subjected to a unit step load FðtÞ ¼ 1 HðtÞ N at its free end (where HðtÞ is the Heaviside function). We note T the observation time and Dt the time step.
The dynamic responses at the free end of the bar are plotted in Fig. 3. These results are in good agreement with those obtained by Enelund et al. [5], thus validating our numerical approach in a traction-compression case. In Fig. 3 (a), we note the influence of the fractional order of the time derivative. When a tends to 1 (solid line), the solution with the fractional derivative algorithm is quite close to that one obtained with the integration method proposed by Taylor et al. [11] for the classical Zener model (line with the sign + at each data point). The case where a ¼ 1 could also be computed by taking only the two first terms A 1 and A 2 in the Grünwald series. To carry out this calculation, we have chosen a time step slightly smaller than 0:5Dt c , where Dt c is the critical time step defined in [5].
In Fig. 3 (b), the time step adopted is greater than Dt c . In this way, contrary to the explicit time integration method used by Enelund et al. [5], our numerical procedure avoids time step limitations due to stability considerations.

Viscoelastic beam
An example with a viscoelastic simply supported beam [3] was chosen to validate our formulation for a Timoshenko beam (sandwich beam without external faces). The length, width and thickness of the beam are L ¼ 10 m, b ¼ 2 m and h ¼ 0:5 m. The shear correction factor adopted is defined by k c ¼ 10ð1 þ mÞ=ð12 þ 11mÞ. The beam is modeled by 10 finite elements, with a uniform step function load of 10 HðtÞ N/m acting all over the beam. The material properties are q ¼ 500 kg/m 3 , m ¼ 0:3, E o ¼ 19:6 MPa, E 1 ¼ 98 MPa and s ¼ 2:24 s. Figure 4 shows the timedependent displacements at the center of the beam for a standard solid model (a ¼ 1) and for a fractional derivative one (a ¼ 0:75 and a ¼ 0:5). For a ¼ 1, our solution is in good agreement with that one taken as reference [3]. This demonstrates that our sandwich model without faces degenerates into a standard Timoshenko beam. As expected, for a decreasing fractional order of the time derivative, damping and the time required to reach the quasi-static long time solution increase considerably.

Results
In this section, we present the example of a cantilevered sandwich beam with viscoelastic core and symmetrical faces, discretized along its length with 5 finite elements. To study the influence of computational model parameters on the structural behavior of the sandwich beam, we present some numerical results.
Firstly, the influence of the time step (Dt ¼ 0:50; 0:25; and; 0:01 ms) on the dynamic tip displacement response of the structure is shown in Fig. 6. Taking the whole history of the anelastic displacements (i.e., N t ¼ 500, 1000, and 2500 terms), we note that the solutions are quite close. This result shows that the solution converges for a decreasing time step, proving the consistency of the present algorithm.
From a numerical point of view, an investigation on the total energy conservation of the system is proposed in Fig. 7 for Dt ¼ 0:5 ms and 500 terms in the Grünwald series. As already discussed in paragraph 4.1, the energy balance is strictly zero (see dotted line in Fig. 7 (a)) since the integration scheme has no supplementary numerical damping. It can be observed that the work done by the external force W increases quickly up to 4 ms and remains constant after. This value of 4 ms is that in the end of the triangular impulse (see Fig. 5) applied to the beam. The sum of the kinetic and elastic energies T þ U, which also reaches its maximum value at 4 ms, decreases with time. At about 200 ms, this sum vanishes, and the total dissipated energy U d À W d remains constant and equals W.
In Fig. 7 (b), the two components of the total dissipated energy are plotted. We note that the major contribution to the total energy dissipated in the system arises from the time history dependent term W d , except in the beginning, when U d is important. This can be explained by the action of the external load, since U d depends on the displacement variation over the most recent time step.
It should be pointed out that U and U d have always the same shape (see Eq. (38)), however, for a given time step, their amplitudes are different. In fact, the amplitude of U d increases when the time step decreases since " K c is inversely proportional to Dt a .
In order to demonstrate that it is not necessary to take a large number of terms in the Grünwald series, a more detailed study of the numerical parameters is carried out below. In Fig. 8 (a), we plot the error of total dissipated energy versus the number of terms in the Grünwald expansion of the fractional operator approximation for various values of Dt. The error is expressed by where D ref is the total dissipated energy associated to the reference solution. The so called reference solution is a priori a very accurate solution with a fine time discretization (Dt ¼ 0:1 ms, i.e. 2500 time steps) and the whole history in the Grünwald expansion (2500 terms).
The basic idea of this strategy consists of choosing, for a given time step, the truncation corresponding to the first local minimum of each curve in Fig. 8 (a). For example, taking a rough time step (Dt ¼ 1 ms or 250 time steps) and 13 terms in the history, an error of about 7% is observed. This error slightly increases when more than 13 terms are retained. Taking 500 time steps (Dt ¼ 0:5 ms), we obtain an error lower than 2% with only 26 terms in the series. Moreover, the first local minimum in Fig. 8 (a) is also the global one and the curves are plotted up to 100 terms in the Grünwald expansion since the error values after this point are quasi-constant even through some small oscillations appear.
The above results lead us to think that there is a truncation time corresponding to each of these minima, where the recent history needs to be stored. In this way, Fig. 8 (b) shows that we must only store the time history over a truncation time in order to obtain the most accurate solution. In this example, the truncation time is quasiconstant (% 13 ms) for various values of Dt. Moreover, refining the time step means decreasing the error but increasing the numerical cost. This might justify to work with an implicit scheme.

Conclusion
A finite element formulation for transient dynamic analysis of sandwich beams with embedded viscoelastic material is proposed. The three-layer sandwich model is built assuming Euler-Bernoulli hypotheses for the elastic faces and Timoshenko ones for the viscoelastic core. The four-parameter fractional derivative model is used to describe the frequency-dependence of the viscoelastic layer, and a strategy to identify the model parameters is discussed. The master curves of a viscoelastic material shows good agreement between the fractional model and experimental data.
A numerical scheme to obtain the dynamic response of the viscoelastic sandwich beam is presented. This direct time integration method is based on the implicit Newmark algorithm in conjunction with the Grünwald approximation for the fractional order derivative of the constitutive equations. The time dependent terms, arising from the viscoelastic model, are shifted to the right-hand side of the equations of motion, modifying in this way the transient excitation. Moreover, only the anelastic displacement history is stored in the algorithm calculation. This implies a substantial reduction of the numerical efforts related to the non-locality of fractional operators. Numerical applications are presented and analyzed. They show the effectiveness of the proposed numerical treatment of fractional order for a suitable truncation of the Grünwald approximation (fading memory phenomena) and for small enough time increments.