A dynamic viscoelastic problem: Experimental and numerical results of a finite vibrating plate

Abstract: In this paper, we numerically study a dynamic viscoelastic problem. The variational formulation is written as a linear parabolic variational equation for the velocity field. An existence and uniqueness result is recalled. Then, fully discrete approximations are introduced using the implicit Euler scheme and the finite element method, for which some a priori error estimates are derived, leading to the linear convergence of the algorithm under suitable additional regularity conditions. Finally, some oneand three-dimensional numerical simulations are presented to show the accuracy of the algorithm and the behaviour of the solution, including a comparison with an experimental study.


ABOUT THE AUTHORS
A. Segade received his PhD degree in Industrial Engineering from the University of Vigo (Spain) in 2008. Currently, he is an associate professor at this university.

PUBLIC INTEREST STATEMENT
Dynamic problems involving viscoelastic materials appear usually in industry. These materials have been utilized in many engineering applications since they can be customized to meet a desired performance while maintaining low cost. An important issue concerning such materials is that they may exhibit time-dependent and inelastic deformations. In this work, we introduce an efficient algorithm for solving a linear case, proving error estimates which allow to show its convergence. Moreover, we compare the numerical results with those obtained experimentally in a vibrating plate.
These viscoelastic materials have been utilized in many engineering applications since they can be customized to meet a desired performance while maintaining low cost. An important issue concerning such materials is that they may exhibit time-dependent and inelastic deformations. The viscoelastic strain component consists of a recoverable-reversible part (elastic strain) and a recoverable-dissipative deformation part (inelastic strain).
In this paper, we assume that the material is linear for the sake of simplicity. The variational formulation is written in terms of the velocity field, leading to a linear parabolic variational equation. Then, its numerical analysis is performed, providing a priori error estimates for the fully discrete approximations. Finally, some numerical simulations are presented in one-and three-dimensional examples to demonstrate the accuracy of the approximation and the behaviour of the solution. In order to check the applicability of the proposed model, an experimental study has been designed, related to the vibration of a metallic plate, comparing the measured results with those obtained from the numerical implementation. This work has application in the modelling of vibrating systems in order to better understand their behaviour and to optimize damping structures.
The outline of this paper is as follows. In Section 2, we describe the mathematical problem and derive its variational formulation. An existence and uniqueness result is recalled. Then, fully discrete approximations are introduced in Section 3 by using the finite element method for the spatial approximation and the implicit Euler scheme for the discretization of the time derivatives. An error estimate result is proved from which the linear convergence is deduced under suitable regularity assumptions. Finally, in Section 4, some one-and three-dimensional numerical examples are shown to demonstrate the accuracy of the algorithm and the behaviour of the solution.

Mechanical problem and its variational formulation
Denote by d , d = 1, 2, 3, the space of second-order symmetric tensors on ℝ d and by " ⋅ " and ‖ ⋅ ‖ the inner product and the Euclidean norms on ℝ d and d .
Let Ω ⊂ ℝ d denote a domain occupied by a viscoelastic body with a Lipschitz boundary Γ = Ω decomposed into two measurable parts Γ D and Γ F , such that meas (Γ D ) < 0. Let [0, T], T > 0, be the time interval of interest. The body is being acted upon by a volume force with density f 0 , it is clamped on Γ D and surface tractions with density f F act on Γ F . Moreover, let = ( i ) d i=1 be the outward unit normal vector. Let x ∈ Ω and t ∈ [0, T] be the spatial and time variables, respectively. In order to simplify the writing, in some places, we do not indicate the dependence of the functions on x and t. Moreover, a dot above a variable represents its first derivative with respect to the time variable and two dots indicate its derivative of second order.
the stress tensor and the linearized strain tensor, respectively. We recall that The body is assumed linearly viscoelastic and it satisfies the following constitutive law (see, for instance, Duvaut & Lions, 1976), where  = (a ijkl ) and  = (b ijkl ) are the fourth-order viscous and elastic tensors, respectively.
We turn now to describe the boundary conditions.
On the boundary part Γ D , we assume that the body is clamped and thus the displacement field vanishes there (and so u = 0 on Γ D × (0, T)). Moreover, since the density of traction forces f F is applied on the boundary part Γ F , it follows that = f F on Γ F × (0, T).
The mechanical problem of the dynamic deformation of a viscoelastic body is then written as follows.
Problem P. Find a displacement field u:Ω × [0, T] → ℝ d and a stress field Here, > 0 is the density of the material (which is assumed constant for simplicity), and u 0 and v 0 are initial conditions for the displacement and velocity fields, respectively. Moreover, Div represents the divergence operator for tensor-valued functions.
In order to obtain the variational formulation of Problem P, let us denote by H = [L 2 (Ω)] d , and define the variational spaces V and Q as follows, where for a Hilbert space X, let (⋅, ⋅) X and ‖ ⋅ ‖ X be the scalar product and norm in X, respectively.
We will make the following assumptions on the problem data.
The following regularity is assumed on the density of volume forces and tractions: Finally, we assume that the initial displacement and velocity satisfy Moreover, we denote by V ′ the dual space of V. We identify H with its dual and consider the Gelfand triple We use the notation < ⋅, ⋅ > V � ×V to denote the duality product and, in particular, we have Using Riesz theorem, from (8) we can define the element f (t) ∈ V � given by and then f ∈ C([0, T];V � ).
Plugging (1) into (2) and using the previous boundary conditions, applying Green's formula, we derive the following variational formulation of Problem P, written in terms of the velocity field v(t) =u(t).
Problem VP. Find a velocity field v:[0, T] → V such that v(0) = v 0 and for a.e. t ∈ (0, T) and for all w ∈ V, where the displacement field u(t) is given by The existence of a unique solution can be obtained proceeding as in Kuttler, Shillor and Fernández (2006), where the contact with a deformable obstacle is also considered and the damage and adhesion effects are included. Thus, we have the following.
Theorem 2.1 Assume (6)-(9) hold. Then, there exists a unique solution v to Problem VP. Furthermore, the solution satisfies

Fully discrete approximations and an a priori error analysis
In this section, we introduce a finite element algorithm for approximating solutions to variational problem VP. Its discretization is done in two steps. First, we consider the finite element spaces where Ω is assumed to be a polyhedral domain,  h denotes a triangulation of Ω compatible with the partition of the boundary Γ = Ω into Γ D and Γ F , and P q (T), q = 0, 1, represents the space of polynomials of global degree less or equal to q in T. Here, h > 0 denotes the spatial discretization parameter.
Secondly, the time derivatives are discretized using a uniform partition of the time interval [0, T], denoted by 0 = t 0 < t 1 < ⋯ < t N = T, and let k be the time step size, k = T∕N. Moreover, for a continuous function f(t) we denote f n = f (t n ) and, for the sequence {z n } N n=0 , we denote by z n = (z n − z n−1 )∕k its corresponding divided differences.
Using the well-known implicit Euler scheme, the fully discrete approximation of Problem VP is the following. Using assumptions (6)-(9) and the classical Lax-Milgram lemma, it is easy to prove that Problem VP hk has a unique discrete solution v hk ⊂ V h .

Problem VP
Our aim in this section is to derive some a priori error estimates for the numerical errors u n − u hk n and v n − v hk n . For the sake of simplicity in the calculations developed below, we assume that f ∈ C([0, T]; V h ) and so, f h n = f n . It is straightforward to extend the results presented in the rest of this section to a more general situation.
Thus, we have the following.
Theorem 3.1 Let assumptions (6)-(9) hold. If we denote by v and v hk the respective solutions to Problems VP and VP hk , then there exists a positive constant C > 0, independent of the discretization parameters h and k, such that, for all ( Proof We subtract variational Equation (10), at time t = t n and for w = w h ∈ V h , and discrete variational Equation (14) to get, for all w h ∈ V h , Therefore, we find that, for all w h ∈ V h , Keeping in mind assumptions (6) and (7), it follows that where u n = (u n − u n−1 )∕k, and since where v n = (v n − v n−1 )∕k, using several times Cauchy's inequality ab ≤ a 2 + 1 4 b 2 , a, b, ∈ ℝ with > 0, and assumptions (6)-(9), we have, for all w h ∈ V h , Thus, by induction we find that, for all Now, taking into account that We will use now the following lemma which constitutes a discrete version of Gronwall's inequality (see Barboteu, Fernández, & Hoarau-Mantel, 2005;Campo, Fernández, Kuttler, Shillor, & Viaño, 2006). where C = C 1 (1 + C 1 Te 2C 1 T ) and T = Nk.
Finally, if we define and G n the remaining terms on the right-hand side of the previous estimates, applying Lemma 3.2 we derive the a priori error estimates (17).

✷
We note that from estimates (17) we can derive the convergence order under suitable additional regularity conditions. For instance, if we assume that the continuous solution has the additional regularity: then we have the following result.  (19), it follows the linear convergence of the solution obtained by Problem VP hk ; that is, there exists a positive constant C, independent of the discretization parameters h and k, such that Notice that this linear convergence is based on some well-known results concerning the approximation by the finite element method (see, for instance, Ciarlet, 1993), the discretization of the time derivatives and the following result (see Barboteu et al., 2005;Campo et al., 2006 for details),

Numerical results
In this final section, we present the numerical scheme which has been implemented in MATLAB (1D) and the commercial code ANSYS (3D) to obtain the solutions to Problem VP hk and then we show some numerical examples to demonstrate its accuracy and its behaviour.

Numerical scheme
Using the finite element spaces defined in (12) and (13), let n = 1, 2, … , N and given u hk n−1 and v hk n−1 ∈ V h , the discrete velocity field v hk n , at time t = t n , is then obtained from Equation (14). That is, we solve the following problem, for all w h ∈ V h , Then, the discrete displacement field u hk n is updated from Equation (15) as follows, We point out that this problem leads to a linear system which is solved using classical Cholesky's method. In the one-dimensional simulations, this numerical scheme was implemented on a 3.2 Ghz PC using MATLAB, and a typical run (h = k = 0.01) took about 0.303 s of CPU time. For the three-dimensional analysis, we used the commercial software ANSYS with a 2.86 Ghz PC, where the CPU time was 43 min and 35 s due to the large size of the finite element mesh and a time step k = 0.001.

A first example: Numerical convergence
As an academical example, in order to show the accuracy of the approximations, we consider the following simpler problem. where the volume force f 0 is given by We note that Problem P ex corresponds to Problem P with the following data: and initial conditions Moreover, the exact solution to Problem P ex can be easily calculated and it has the form, for (x, t) ∈ (0, 1) × (0, 1), To show the numerical convergence and the asymptotic behaviour of the algorithm, the numerical errors given by are calculated and presented (multiplied by 10 2 ) in Table 1 for several values of the discretization parameters h and k. Finally, the evolution of the error depending on the parameter h + k is plotted in Figure 1. We notice that the convergence of the algorithm is clearly observed, and the linear convergence, stated in Corollary 3.3, is achieved.

A second example: Comparison with experimental results
In the second example, we consider a three-dimensional setting. Then, the aim is to simulate a finite vibrating plate made of aluminium. The plate occupies the domain Ω = (0, 35) × (0, 2) × (0, 250) (where the unit is the millimetre), and we study its deformations during 0.35 s (i.e. T = 0.35 s). We assume no body forces acting in Ω and that the body is clamped on the part   Figure 2). An initial displacement is prescribed moving point x = (17.5, 0, 210) 3.2 mm in the Y-direction. We note that no surface forces are then applied in this example and that we simulate the action of the sensor, used to obtain the experimental measures, with an additional mass of 13 g located at its position (see Figure 3).
The following data were used in the simulations: where the initial displacement u 0 is the deformation of the plate corresponding to the prescribed initial displacement.  The elastic tensor  is defined using Young's modulus E and Poisson's ratio ; i.e.
where = ( ij ) 3 i,j=1 and ( ij ) 3 i,j=1 denote a three-dimensional symmetric tensor and the classical Kronecker delta, respectively. Moreover, since we are simulating a metallic plate made of   aluminium, we take values E = 72,000 MPa and = 0.3. Moreover, its density is assumed to be 2,770 Kg∕m 3 .
For the definition of the viscous part , we employ a damping coefficient with respect to the elastic part, and so  = . In these simulations, we use the value = 1.0954 × 10 −4 .
Using a finite element mesh composed of 35,611 nodes and 6,804 elements (hexahedra), which is shown in Figure 4, and time step k = 10 −3 , Problem VP hk is now solved using the well-known commercial software ANSYS, which we checked it uses a similar implementation as described in the previous section (the unique difference is that it is based on a displacement-type formulation, as a result of plugging expression v hk n = (u hk n − u hk n−1 )∕k into the discrete variational Equation (14)).
Therefore, in Figure 5, the deformation of the plate is plotted at final time. As expected, due to the vibration process, a bending is produced in the Y-direction.
Moreover, the viscoelastic stresses (von Mises norm) are plotted in Figure 6. Obviously, due to the clamping conditions and since no external forces are applied, they concentrate around the clamped area.
From these results, we find a numerical frequency hk = 182.131 rad/s and a numerical coefficient of relative damping z hk = 0.0099754.
Finally, our aim is to compare the above numerical results with those obtained in an experimental study. Then, a metallic plate is subjected to a prescribed initial displacement as we did with the numerical simulation (see Figure 7). As can be seen, the plate made of aluminium is clamped on a part of its right-hand side, and there is a sensor, which we simulated numerically as an added mass, and a prescribed initial displacement on the left-hand part.
Thus, we perform some experimental studies obtaining (in mean values) a natural frequency = 180.453 rad/s and a coefficient of relative damping z = 0.0098830. As we can see, both values are similar with those obtained numerically. As an example, the results obtained in one of these experiments are plotted in Figure 8. There we show the evolution in time of the voltage measured by the sensor located at point x = (17.5, 0, 210). As expected, there is a vibration of the plate and, due to the viscoelastic (dissipative) behaviour, there is a decrease in its amplitude. Finally, in Figure 9, we plot the displacement in the Y-direction at point x = (17.5, 0, 210) obtained experimentally (blue curve) and numerically (orange curve). As can be seen, both almost coincide.