Vibrations of a beam between stops: numerical simulations and comparison of several numerical schemes

Mathematics and Computers in Simulation 60 (2002) 45–83

Vibrations of a beam between stops: numerical simulations and comparison of several numerical schemes夽 Yves Dumont IREMIA, Université de La Réunion, 15 Avenue René Cassin, Saint-Denis, 97715 Ile de La Réunion, France Received 10 October 2001; accepted 24 November 2001

Abstract We present and compare different schemes used to simulate the vibrations of an elastic beam between two stops. The beam is clamped at one end to a device which may vibrates periodically. The other end is restricted to move between two stops. The contact is modeled by the normal compliance condition which describes nonlinear stops and approximates the Signorini condition. We use the method of lines to obtain numerical algorithms. © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: (Visco)elastic beam; Noise characteristic; Normal compliance; Constrained vibrations; Method of lines

1. Introduction The purpose of this paper is to present and compare different numerical schemes used to simulate the vibrations of an elastic beam between two flexible (rigid) stops. The model and existence of weak solutions have been studied in [9]. There is considerable interest in industry in dynamic vibrations of mechanical systems. In the automotive industry the noise and vibration characteristics of cars and car components are considered as an important factor in the product customer satisfaction. Indeed, appreciable effort has been made in designing automotive components to reduce unwanted or disturbing noise. Some of the unwanted noise is generated by the dynamic contact of parts and components when periodically forced. For instance, if the mounting of the components on the car engine is not perfect, the motion in the clearances leads to dynamic contact, which in turn may generate unwanted noise. Another important problem is to analyze the wear, resulting from impacts, of some components of pressurised water reactor in order to predict the life-time of this components (combustible bar, etc.). The problem we consider is one-dimensional and describes a viscoelastic or an elastic beam that is clamped at one end to a possibly vibrating device while the other end oscillates between two rigid or 夽

This paper is dedicated to my daughter Joanne. E-mail address: [email protected] (Y. Dumont). 0378-4754/02/$ – see front matter © 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 2 ) 0 0 0 0 3 - 4

46

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

flexible stops. Such a setting was considered by Moon and Shaw [12] (also see [11] and references therein) but in a simplified form. Indeed, they exhibited complicated behavior such as quasiperiodic or chaotic oscillations. Here, we consider the full problem and we expect to observe chaotic behavior. Using the well-known method of lines (MOL), we also consider different spatial and temporal discretizations associated to the elastic case. Through these discretizations, we are interested to compare the different schemes we obtain and to discuss which is the most appropriate for solving this type of problems. We also consider time integration formulas that are commonly used in this type of problems, i.e. backward differentiation formula (BDF) and the Hilber–Hughes–Taylor’s-␣ (HHT-␣) method. But, we also consider a new time-integration algorithm related to the Runge–Kutta time-integration family and presented in [13] as a possible alternative to the use of the HHT-␣ algorithm for the time integration of motion in structural dynamics. From our knowledge, this algorithm has never been used to integrate such type of problems and so it seems to be interesting the study its numerical behavior. The outline of the paper is as follows: the model is described in Section 2, following [9]. We consider that the contact is modeled with the normal compliance condition, which models flexible stops with spring-like nonlinear reaction; this model is mechanically and mathematically the most attractive, e.g. to approximate Signorini’s condition (see, e.g. [1,2,8–10] and references therein). In Section 3, we use the MOL in order to present different space discretizations available for the elastic beam’s case. We also obtain stiff (oscillatory) ordinary differential equations (ODEs) systems. In Section 4, we also present different time integrations of the previous semidiscretized equations. Finally, in Section 5, we give different simulations and compare the results obtained using the different algorithms. 2. The model In this section, we present the classical model and its variational formulation, following [9]. Then we quote the existence results obtained there. The mechanical setting, depicted below, is as follows. A viscoelastic or an elastic beam is clamped at its left end to a device that may shake or oscillate. Its right end is free to move vertically, but this motion is constrained by two obstacles: the stops.

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

47

The area-center of gravity of the beam in its (stress free) reference configuration coincides with the interval 0 ≤ x ≤ 1; g1 and g2 (g1 < 0 < g2 ) are the positions of the stops. We set ΩT = (0, 1)×(0, T ), for T > 0. Let u = u(x, t), (x, t) ∈ ΩT , represent the vertical displacement of the beam, and σ (x, t) = σT (x, t) the shear stress. Then, the equation of motion, in nondimensional form, is utt − σx = f

in ΩT ,

(2.1)

where f denotes the density (per unit length) of applied forces, and the subscripts x and t indicate partial derivatives. The material is assumed to be viscoelastic, of constant cross section and we consider the Kelvin–Voigt law of viscoelasticity for shear stress σ (x, t) = −k 2 uxxx (x, t) − duxxxt (x, t), where k 2 is the scaled elastic modulus and d is a nonnegative viscosity coefficient. When d = 0 the material is elastic. The initial conditions take the form u(x, 0) = u0 (x),

ut (x, 0) = v0 (x),

(2.2)

for 0 ≤ x ≤ 1. The beam is rigidly attached at its left end, u(0, t) = φ(t)

and

ux (0, t) = 0,

(2.3)

where φ = φ(t) represents the motion of the supporting device. In this paper we present some simulation results when the device oscillates periodically, i.e. when φ is a periodic function (also see [11]). At the free end we choose the normal compliance condition (see, e.g. [7–9]): the stops are assumed to be flexible, with resistance force proportional to the deflection σ (1, t) = −κ[(u(1, t) − g2 )+ − (g1 − u(1, t))+ ].

(2.4)

We thus assume that the stops behave as one-sided nonlinear springs, with spring constant κ, a positive constant. Remark 1. When κ → ∞, we obtain the non-penetration Signorini condition, which models rigid stops [4,9]. Thus (2.4) may be considered as its regularization. Finally, we assume that no moments act on the free end, k 2 uxx (1, t) + duxxt (1, t) = 0.

(2.5)

The classical problem with normal compliance is: find a function u such that (2.1)–(2.5) hold. Consider now the following spaces: H = L2 (0, 1) and V = {w ∈ H 2 (0, 1) : w(0) = w (0) = 0}, then V ⊂ H = H ⊂ V , where V is the topological dual of V , so {V , H, V } is a Gelfand triplet. We denote by (·, ·) the inner product in H and by ·, · the duality pairing between V and V . The associated norm on H is denoted by | · |H and the one on V by || · ||V . Here and below “ ” denotes the time derivative. Under a suitable change of variable, i.e v(x, t) = u(x, t)−φ(t)(1−x 2 ), w0 (x) = u0 (x)−φ(0)(1−x 2 ) and z0 (x) = v0 (x)−φ (0)(1−x 2 ), the following existence result for the problem with normal compliance was etablished in [9].

48

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Theorem 2. Assume that f ∈ L2 (ΩT ),

w0 ∈ V ,

z0 ∈ H,

φ ∈ H 2 (0, T ).

Then problem (2.1)–(2.5) has a weak solution for each T < ∞. When d > 0, the solution is unique and satisfies v ∈ L∞ (0, T ; V ),

v ∈ L2 (0, T ; V ) ∩ L∞ (0, T ; H ),

v

∈ L2 (0, T ; V ).

When d = 0, v ∈ L∞ (0, T ; V ),

v ∈ L∞ (0, T ; H ),

v

∈ L2 (0, T ; V ).

Remark 3. When d = 0, the uniqueness of the weak solution remains an open question. We now use the MOL to approximate numerically the solutions of the PDEs system (2.1)–(2.5): the MOL-approach permits various possibilities for the space and time discretizations. Thus, we first discretize in space to obtain a system of first or second order equations in time; then, we discretize in time by using well-suited methods for stiff (oscillatory) systems. 3. Space discretization We now present different space discretizations available when d ≥ 0 and d = 0: as in [3,4], we use the finite difference method, but it is possible to use finite elements too. We first consider the following change of variable v(x, t) = u(x, t) − φ(t)(1 − x 2 ), which makes the displacements at x = 0 vanish, for t > 0, but does not affect the contact conditions at x = 1. Thus, problem (2.1)–(2.5) can be written as 4 ∂ 2 v(x, t) ∂ 5 v(x, t) 2 ∂ v(x, t) + k + d = f (x, t) − φ

(t)(1 − x 2 ), ∂t 2 ∂x 4 ∂x 4 ∂t

v(x, 0) = u0 (x) − φ(0)(1 − x 2 ), v(0, t) = 0, vx (0, t) = 0, vt (0, t) = 0, vt (x, 0) = v0 (x) − φ (0)(1 − x 2 ), k 2 vxx (1, t) + dvxxt (1, t) = 2(k 2 φ(t) + dφ (t)). k 2 vxxx (1, t) + dvxxxt (1, t) = κ[(v(1, t) − g2 )+ − (g1 − v(1, t))+ ] Let t = T /N be the time step and let h = 1/M be the space step, then the grid points are given by (xi = ih, tn = n t), for i = 1, . . . , M − 1 and n = 0, . . . , N.

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

49

3.1. First space discretization Let d ≥ 0 and V the vector w = (v, u) = (v(x, t), vt (x, t)). Then, the equation in the interior of the beam can be written as ∂v(x, t) = u(x, t), ∂t ∂u(x, t) ∂ 4 v(x, t) ∂ 4 u(x, t) = −k 2 − d − φ

(t)(1 − x 2 ) + f (x, t). ∂t ∂x 4 ∂x 4 We discretize the previous equation with respect to the space variable using the classical fourth order finite difference scheme. We obtain a system of 2M − 4 ordinary differential equations (i = 2, . . . , M − 1):  dvi   = ui ,    dt    2    dui = − k (v − 4v + 6v − 4v + v ) i−2 i−1 i i+1 i+2 dt h4 (3.1)    d   − 4 (ui−2 − 4ui−1 + 6ui − 4ui+1 + ui+2 )   h     −φ

(t)(1 − xi2 ) + fi (t), where ui = u(xi , .), vi = v(xi , .) and fi = f (xi , .). Then, using the initial conditions v(0, t) = 0 and vt (0, t) = 0, we set v0 = 0 and u0 = 0. Moreover, from the assumption vx (0, t) = 0, we deduce  dv1   = u1 ,    dt    2    du1 = − k (7v − 4v + v ) 1 2 3 dt h4 (3.2)    d   − 4 (7u1 − 4u2 + u3 )   h     −(1 − h2 )φ

(t) + f (x1 , t) For i = M, recalling that k 2 vxxx (1, t) + duxxx (1, t) = −σ (1, t) and k 2 vxx (1, t) + duxx (1, t) = 2(k 2 φ(t) + dφ (t)), and using Taylor’s expansion at the virtual points x = 1 + h and x = 1 + 2h, we have  du k2    M = − 4 (2vM − 4vM−1 + 2vM−2 )   dt h    d (3.3) − (2uM − 4uM−1 + 2uM−2 )  h4      2σ 4   −f (xM , t) + + 2 (k 2 φ + dφ ) h h

50

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Then, taking i = M − 1 in (3.1) and using again the assumption k 2 vxx (1, t) + duxx (1, t) = 2(k 2 φ(t) + dφ (t)), we obtain  dvM−1     dt      duM−1       dt

= uM−1 , =

              

k2 (vM−3 − 4vM−2 + 5vM−1 − 2vM ) h4 d − 4 (uM−3 − 4uM−2 + 5uM−1 − 2uM ) h +f (xM−1 , t) − φ

(t)(1 − (1 − h)2 ) −

−

(3.4)

2 2 (k φ + dφ ) h2

Finally, from (3.1)–(3.4), we obtain the following ODE system y = Kd y + F1 + G1

(3.5)

where 

0

  7k 2 −  h4    0    4k 2   h4    0   Kd =  ..  .   .  .  .    ..  .    .  ..    0

1

0

0

0

0

7d − 4 h 0

4k 2 h4 0

4d h4 1

k2 − 4 h 0

d − 4 h 0

6d h4 0

4k 2 h4 0

4d k2 − h4 h4 1

4d h4 0

−

6k 2 h4 0

−

0

0

0

0

0

0

0

0

0

0

0

0

d h4 0

0

0

0

···

−

0 .. .

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

d h4

4k 2 h4

4d h4

5d h4

2k 2 h4

2d h4

0

0

0

1

···

−

k2 h4

−

···

0

0

···

0

0

−

2k 2 h4

−

is a 2M × 2M matrix, y is the column vector y = (v1 , vt,1 , v2 , vt,2 , . . . , vM−1 , vt,M−1 , vM , vt,M )T ,

2d h4

−

5k 2 h4

−

0

0

4k 2 h4

4d h4

−

.. .

.

2k 2 h4

−

2d h4

                  ,                

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

in R2M , F1 is a bounded function 

51



0

    −(1 − h2 )φ

(t) + f (h, t)       0     .   ..       F1 =  0       −φ

(t)(1 − (1 − h)2 ) − 2 (k 2 φ + dφ ) + f (1 − h, t)    h2       0     4 2

(k φ + dφ ) + f (1, t) h2 and G1 is a nonlinear function defined in the following way   0   ..     .   G1 =     0    2  − κ((vM (t) − g2 )+ − (g1 − vM (t))+ ) h 3.2. Second space discretization The previous space discretization can be rewritten in the following system of second order ODEs y

= k 2 Hy + d Hy + F2 + G2 , (3.6) where y = (u1 , . . . , uM ) is the vector displacement of the beam, H is an M × M matrix   −7 4 −1 0 0 0    4 −6 4 −1 0 0        .. .. .. .. . 0 . . .  0   ..  1  .. .. .. ..  .. . . . . . H= 4 . ,  h   . ..  .. .. .. ..  .. . . . . .       .   .. · · · −1 4 −5 2   0

···

···

−2

4

2

52

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83





−(1 − h2 )φ

(t) + f (h, t)

  ..     .     2

  −φ (t)(1 − xi ) + f (xi , t)    F2 =  ..     .      −φ

(t)(1 − (1 − h)2 ) − h22 (k 2 φ + dφ ) + f (1 − h, t)      4 2

(k φ + dφ ) + f (1, t) h2   0   ..     .   G2 =     0    2  − κ((uM (t) − g2 )+ − (g1 − uM (t))+ ) h 3.3. Third space discretization When d = 0, there is another way to obtain a space discretization of the problem: we rewrite it in a velocity-curvature formulation. This formulation was used in [3]. Denote by w the vector (vt , vxx ) = (vt (x, t), vxx (x, t)). Then, the equation in the interior of the beam becomes ∂w ∂2 = 2 (Aw) + F, ∂t ∂x where A=



0

−k 2

1

0

(3.7)



 and

F =

F1 F2



 =

f (x, t) − φ

(t)(1 − x 2 ) 0

.

We now give a spatial discretization of (3.7) using a classical second order finite difference formula. We thus obtain the following system of M ODEs dw i 1 = 2 A(wi+1 − 2w i + w i−1 ) + Fi , dt h where wi = (vt (xi , t), vxx (xi , t)), for 1 ≤ i ≤ M and t ∈ [0, T ]. The initial conditions are vt (x, 0) = v0 (x) − φ (0)(1 − x 2 ) and vxx (x, 0) = u

0 (x); thus, 0 = v0,i − φ (0)(1 − (ih)2 ) vt,i

and

0 vxx,i = u

0,i .

The boundary conditions vxx (1, t) = 0 and vt (0, t) = 0 are n vxx,M = 0,

n vt,0 = 0,

(3.8)

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

53

respectively. Next, we derive a boundary condition for vxx at x = 0. Recalling that vx (0, t) = 0 and using Taylor’s expansion at x = h and x = 2h yield vt (h, t) =

h2 h3 vxxt (0, t) + vxxxt (0, t) + O(h4 ), 2 6

vt (2h, t) =

4h2 8h3 vxxt (0, t) + vxxxt (0, t) + O(h4 ) 2 6

We deduce 1 d vxx,0 = 2 (8vt,1 − vt,2 ) + O(h2 ) dt 2h

(3.9)

We turn to the discretization of the contact condition. First, we use Taylor’s expansion of vxx at x = 1 − h and 1 − 2h vxx (1 − h, t) = vxx (1, t) − hvxxx (1, t) +

h2 h3 vxxxx (1, t) − vxxxxx (1, t) + O(h4 ), 2 6

vxx (1 − 2h, t) = vxx (1, t) − 2hvxxx (1, t) +

4h2 8h3 vxxxx (1, t) − vxxxxx (1, t) + O(h4 ). 2 6

Using the assumption vxx (1, t) = 2φ(t), we deduce d −k 2 3k 2 7k 2 vt,M = v (8v − v ) − + φ(t) xx,M−1 xx,M−2 xxx,M dt 2h2 h h2 Recalling that σ (1, t) = σM = −k 2 vxxx (1, t), we have −k 2 3 7k 2 d vt,M = (8v − v ) + + φ(t), σ xx,M−1 xx,M−2 M dt 2h2 h h2

(3.10)

where σM = −κ((vM − g2 )+ − (g1 − vM )+ ),

(3.11)

and vM is the displacement at the end of the beam. Then, from (3.8)–(3.10), we obtain the following system of ODEs y = M y + F 3 + G3

(3.12)

54

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

where



0

  −k 2    0    0    0  1   M = 2  ... h    .  ..    ..  .    ..  .  0

4

0

−1/2

0

0

0

0

···

0

2k 2

0

−k 2

0

0

0

···

−2

0

1

0

0

0

0

···

0

−k 2

0

2k 2

0

−k 2

0

···

1

0

−2

0

1

0

0

···

..

..

.

..

.

..

.

..

.

..

.

.. .

..

.

..

.

..

.

..

.

.

..

···

0

0

0

1

0

−2

0

···

0

0

0

0

k 2 /2

0

−4k 2

..

.

.

0

is a 2M × 2M matrix, y is the column vector y = (vxx,0 , vt,1 , vxx,1 , vt,2 , . . . , vxx,M−2 , vt,M−1 , vxx,M−1 , vt,M )T , in R2M , F3 is a bounded function 

0



    f (h, t) − φ

(t)(1 − h2 )       0       ..   .       0     2

  f (x , t) − φ (t)(1 − x ) i   i     0 F3 =       . ..         0       2  f (1 − h, t) − φ

(t)(1 − (1 − h)2 ) − 2k φ(t)    h2       0     2   7k f (1, t) + 2 φ(t) h

0



 0   0   0   0   ..  , .   ..  .        1  0

(3.13)

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

55

and G3 is the following nonlinear function   0   ..     .   G3 =  .   0    3κ  − ((v(1, t) − g2 )+ − (g1 − v(1, t))+ ) h We note that M is a block matrix with diagonal and off-diagonal matrices   0 −k 2 0 2k 2 , , −2 0 1 0 respectively, except for the first and last rows and the first and last columns. Remark 4. Note that we have to compute the displacement every time step, because the unknows are the velocity and the curvature. The previous discretization could be useful to study the bending of the beam without computing the displacement of the whole beam. Remark 5. Both discretizations in space are second order discretizations, but note that the second one gives only a system of M differential equations. For the rest of the paper, we focus on the case d = 0, i.e. the elastic beam. 4. Time discretization of problems (3.5), (3.6) and (3.12) We now propose and discuss different time integration formulas. From the previous space discretizations, we obtain very stiff (and/or high oscillatory) initial value problems



y (t) = f (t, y) (4.1) y(0) = y0 or 



 y (t) = f (t, y, y )  y(0) = y0 (4.2)  

y (0) = z0 that can be rewritten in the following manner



y (t) + Ay(t) = F (t) + G(t, y), t > 0 y(0) = y0 or



y (t) + Ay(t) = F (t) + G(t, y),   y(0) = y0  

y (0) = z0

(4.3)

t >0 (4.4)

56

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

where y, y and y

are vectors in RM or in R2M while F and G define M (2M) linear and nonlinear components of the systems, respectively. Note carefully that the stiffness comes both from the space discretization and the normal compliance condition. From the spatial discretizations, we obtain systems that are (very) stiff and nonlinear (normal compliance condition). More than the inner frequencies of the beam, some frequencies, and possibly high frequencies, appear when the beam hit the stops. The presence of high frequencies in the solution components require methods that damp out these high frequencies through numerical dissipation, i.e. L-stable methods. There exists various integration formulas available to solve our systems, but we focus in particular on algorithms widely used in structural dynamics, such as Newmark’s method and the HHT-␣ method, or used for the same type of problems, i.e. BDF method [15]. We also propose another method which belongs to the family of Runge–Kutta methods [13] and, from the author’s knowledge, is less used for the time integration in structural dynamics. In the following section we also present the previous time integration formula adapted to the different spatial discretizations. 4.1. Time integration of problem (3.12) 4.1.1. The second order backward differentiation formula In [3], we only use the second order BDF (BDF2) method. It is well known that the BDF2 method or Gear’s method, is unconditionally stable, L-stable and thus suitable for stiff problems. Remark 6. Because M has only purely imaginary eigenvalues, it is not possible to use “classical” high order BDF methods, i.e. up to three, because they are only A(α)-stable. Using (4.1), the BDF2 method is defined as follows: 4 1 2t yn+1 = yn − yn−1 + f (yn+1 , tn+1 ) 3 3 3

(4.5)

Thus from (3.12), each time step we have to solve the following linearized system BV n+1 =

 4 n 1 n−1 2t  n+1 IV − IV , + F3 + Gn+1 3 3 3 3

where B = I − 2/3t M, I is the identity matrix and V n ∈ R2N is a vector defined by  n T n n n n vxx,0 vt,1 vxx,1 vt,2 vxx,2 ··· Vn = . n n n n n · · · vt,M−2 vxx,M−2 vt,M−1 vxx,M−1 vt,M n n Note that vt,0 and vxx,M are omitted, since both vanish. Next, F3n and Gn3 ∈ R2M are given by

F3n+1 = F3 (tn+1 ),       T 3κ  n+1 n+1 n+1 , vM − g2 + − g1 − vM + G3 = 0, . . . , 0, − h

(4.6)

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

57

n−1 n−1 n n where the last term of Gn+1 is computed from the known values vM , vM , vt,M and vt,M . This is the 3 linearization of the contact condition and therefore of the problem.

Remark 7. Note carefully, that the BDF2 method computes only the right-hand side of (3.12) at time tn+1 . We now turn to another important family used to integrate (non)linear ODEs systems: the Runge–Kutta family. Instead of using the “classical” formulae of Runge–Kutta, we present here a subclass of Runge– Kutta family introduced by Owren and Simonsen [13]. 4.1.2. A single diagonal implicit Runge–Kutta method Owren and Simonsen present in [13] a subclass of Runge–Kutta family (SDIRK-OS), which possess properties needed to integrate structural dynamics systems. Thus, they should become a very interesting alternative of other schemes [13]. An s-stage single implicit Runge–Kutta (SDIRK) method for the ODE system (4.1) can be written in the form   i s     Ki = f tn + ci t, yn + t aij Kj , (1 ≤ i ≤ s), bj Kj . (4.7) yn+1 = yn + t j =1

j =1

In the SDIRK methods, we have a11 = · · · = aii = · · · = ass = γ = 0 and we usually impose that  i j =1 aij = ci , which implies c1 = γ . In the litterature, the following convenient representation, called the Butcher array, of the method is used

The SDIRK methods have good stability properties and are relatively cheap to implement. Here, we focus on a SDIRK methods which are stiffly accurate and L-stable, such that the numerical dissipation can be controlled by a parameter. Stiff accuracy recquire that asj = bj ,

j = 1, . . . , s.

In [13], Owren and Simonsen present several schemes that present the previous properties. It exists only one three-stages, second order SDIRK method which is A-stable, stiffly accurate and L-stable. This method is given by the following Butcher array

(4.8)

58

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

where γ is a well choosen parameter and γ 2 − 2γ + 1/2 , b1 = 1 − γ − b 2 . σ The stability function of the two-parameter family of second order methods is b2 =

R(z) =

(4.9)

(3γ 2 − 3γ + 1/2)z2 − (3γ − 1)z + 1 . (1 − γ z)3

We deduce that the methods is L-stable for all γ in the interval [α, β], where α and β are the real roots of the polynomial −6γ 3 + 18γ 2 − 9γ + 1, i.e. α ≈ 0.1804 and β ≈ 2.1856. It is possible to choose σ and γ in such a way that the method becomes three order method for linear and nonlinear ODEs systems. Owren and Simonsen propose to choose σ =

γ 3 − 3γ 2 + 2γ − 1/3 , γ 2 − 2γ + 1/2

and γ the middle root of p(γ ) =

1 6

− 23 γ + 3γ 2 − γ 3 ,

that is γ is approximatively equal to 0.4359. Applying (4.8) to problem (3.12), every time step, we have the following system to solve     n+γ   n+γ  F3 G3 I − tγ M 0 0 K1     n+γ +σ   n+γ +σ   −σ t M I − tγ M   K2  =  F +G  0     3   3  K3 −b1 tγ M −b2 tγ M I − tγ M F3n+1 Gn+1 3 n+γ

n+γ

where, for example, F3 = F3 (tn + tγ ). Like in the BDF2 schemes, the nonlinear term in G3 computed from known values. Then we deduce Yn+1 = Yn + t (b1 K1 + b2 K2 + γ K3 ).

is

(4.10)

Now it is possible to control the numerical dissipation or the damping process in higher modes through the parameter γ . Higher is the value of γ , higher is the damping. The difficulty is to find the appropriate value for γ , but this depends mainly on the problem. Thus, we have now a second order method which permits to control the damping. Moreover, the method is self-starting because of its one-step nature. Remark 8. In the forthcoming simulations we only choose γ ∈ [α, 1] in order to avoid analytic extension of the solution beyond the current step. Remark 9. Owren and Simonsen present other “modified” SDIRK methods and in particular a four-stages three order method. Unfortunately, this method does not introduce sufficient damping to prevent a very bad behavior of the solution. This is why we only consider the three-stages method. Remark 10. The right-hand side is computed at times tn + γ t, tn + (γ + σ )t and tn+1 , which increase the time of computations.

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

59

4.2. Time integration of problem (3.5) We now apply the previous methods to integrate (3.5). 4.2.1. The second order backward differentiation formula At each time step, we have the following linearized system to solve Md Y n =

4 n−1 1 n−2 2t n Y − Y + (F1 + Gn1 ), 3 3 3

(4.11)

where Md = I −

2t BKd , 3

I is the identity matrix and Y n ∈ R2N is a vector defined by  n T n n y1 yt,1 y2n yt,2 ··· ··· ··· Yn = . n n n n · · · · · · yM−1 yt,M−1 yM yt,M n Note that y0n and yt,0 are omitted, since both vanish. Next, F n and Gn ∈ R2M are defined by

F1n = F1 (tn ),       T 2κ  n n n , G1 = 0, . . . , 0, − yM − g2 + − g1 − yM + h n−1 n−1 n n , vM , vt,M and vt,M . where the last term of Gn1 is computed from the known values vM

4.2.2. A single diagonal implicit Runge–Kutta method We also apply the SDIRK method (4.7)–(4.9) described above. Every time-step, we have the following system to solve:     n+γ   n+γ  F1 I − tγ Md 0 0 K1 G1        −σ t Md   K2  =  F n+γ +σ  +  n+γ +σ  I − tγ Md 0     1   G1 Gn+1 K3 −b1 tγ Md −b2 tγ Md I − tγ Md 1 F1n+1 Then we compute Yn+1 = Yn + t (b1 K1 + b2 K2 + γ K3 ),

(4.12)

where γ , σ are given parameters and the coefficients b1 , b2 are computed through (4.9). 4.3. Time integration of problem (3.6) 4.3.1. The HHT-α method The previous BDF and SDIRK-OS methods are of second order which leads to think to use an new second order method.

60

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

In order to solve problem (3.6), we first consider HHT-␣ methods [6], which is defined in the following way: let (yn , zn , an ) be an approximation of (y(tn ), y (tn ), y

(tn )) in (4.4), then  an+1 = (1 + α)k 2 Hyn+1 − αk 2 Hyn + Fn+1 + Gn+1 ,    yn+1 = yn + tzn + (t)2 (βan+1 + (1/2 − β)an ), (4.13)    zn+1 = zn + t (γ an+1 + (1 − γ )an ), 0 ≤ n ≤ N − 1 where α, β and γ are parameters. The previous parameters permit to control the stability, the accuracy and the damping. By choosing particular values of γ and β, we reduce the HHT-α family to a one parameter family, like the previous SDIRK method. In particular, it is known that the HHT-␣ algorithm [6] is as follows: • quadratically convergent if γ = 1/2 − α and β = 1/4(1 − α)2 and • unconditionaly stable for linear problems if α ∈ [−1/3, 0]. The parameter α permits to control the damping of the higher modes in the solution. When α = 0, we recover the well-known Newmark’s trapezoidal rule, which give a second order, non-dissipative and energy preserving algorithm. Unfortunately, this second order Newmark’s scheme does not give satisfactory results because of the presence of high-frequency oscillations: the trapezoid method cannot simultaneously give quadratic accuracy and damps the high-frequency effects. Remark 11. One order Newmark’s schemes, i.e. α = 0 and γ > 1/2, introduce some numerical damping in order to stabilize the integration. Unfortunately, they introduce amplitude decay and period elongation which affect the low modes of the solutions. In [6], it was showed that the HHT method has better dissipation properties than the Newmark’s method. 4.3.2. A single diagonal implicit Runge–Kutta method (SDIRK) Runge–Kutta methods can be a very interesting alternative of the use of HHT-␣ method: in [13], Owren and Simonsen showed that the SDIRK-OS methods have superior stability properties than the HHT-␣ family; furthermore, the SDIRK-OS methods present relative period errors smaller than the HHT. If we apply the SDIRK method (4.7)–(4.9) to (3.6), we have, for i = 1, . . . , s, to solve the systems   i−1  (I − (tγ )2 k 2 H)ki = F2 (tn + ci t) + k 2 H yn + ci ty n + (t)2 Aij kj  , j =1

with respect to ki . Note that Aij =

i 

aik akj .

k=1

Then we compute the approximations yn+1 and yn+1 to y(tn+1 ) and y (tn+1 ), respectively, from the formulae

yn+1 = yn +

ty n

+ (t)

2

s  j =1

Bj kj ,

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

yn+1 = yn + t

s

s 

61

bj kj ,

j =1

where Bj = k=j akj bk . Thus, every time step, we have the following system to solve:     n+γ F2 0 0 I − (tγ )2 k 2 H K1        K2  =  F n+γ +σ  −σ (t)2 k 2 H I − (tγ )2 k 2 H 0     2 2 2 2 2 2 2 K3 −b1 (tγ ) k H −b2 (tγ ) k H I − (tγ ) k H F2n+1 and then we compute

Yn+1 = Yn + tY n + (t)2 (B1 K1 + B2 K2 + B3 K2 ),

Yn+1 = Yn + t (b1 K1 + b2 K2 + γ K3 ).





n+γ

G2



  n+γ +σ  +G    2  Gn+1 2

(4.14)

Now, as before, it is possible to control the numerical damping int the higher modes through the parameter γ which belongs to [α, β], the real roots of the polynomial 24γ 4 −72γ 3 +48γ 2 −12γ +1, i.e. α ≈ 0.1804 and β ≈ 2.1856. Remark 12. In the forthcoming simulations we only choose γ ∈ [α, 1] in order to avoid analytic extension of the solution beyond the current step.

5. The simulations All the simulations presented in this section deal only with the elastic case which is the more difficult and thus the more interesting from √ the numerical √ point of view. Along the simulation, we consider two values for the stiffness, i.e. k = 300 and k = 8000. These values are related to large and small values for the Young modulus, which correspond to more or less flexible beams. 5.1. Solution without contact Let us first consider the no-contact case. Thus, our problem becomes  utt (x, t) + k 2 uxxxx (x, t) = f (x, t), 0 < x < 1, t > 0       v(0, t) = φ(t),       vx (0, t) = 0,    v(x, 0) = u0 (x),     vt (x, 0) = v0 (x),       vxx (1, t) = 0,      vxxx (1, t) = 0.

62

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Under suitable assumptions on the initial conditions, it is possible to give a series solution of the previous problem. Indeed, the solution is given by the following eigenfunction expansion v(x, t) = φ(t) +

∞ 

Xi (x)Ti (t),

i=1

where the eigenfunctions Xi are given by sin λi + sinh λi ( cosh λi x − cos λi x) + ( sin λi x − sinh λi x), cos λi + cosh λi

Xi (x) =

where λi , i = 1, 2, . . . , are positive solutions of the nonlinear equation cos λi +

1 =0 cosh λi

The coefficients Ti verify Ti

(t)

+

λ4i k 2 Ti (t)

= 1 ( 0 Xi2 dx)2

with 1

Ti (0) =  1 ( 0 Xi2 dx)2 Ti (0)

1

1



1



1

(f (x, t) − φ

(t))Xi (x) dx,

0

(u0 (x) − φ(0))Xi (x) dx

0



= 1 ( 0 Xi2 dx)2

1

(v0 (x) − φ (0))Xi (x) dx

0

We now consider a test case in order to show stability and accuracy of our algorithms. Let    f (x, t) = k 2 λ41 − 4 X1 (x) sin 2t,       φ(t) = 0,       u0 (x) = 0,    v0 (x) = 2X1 (x),     vx (0, t) = 0,       v (1, t) = 0,    xx   vxxx (1, t) = 0. Thus the exact solution u is given by u(x, t) = X1 (x) sin 2t and λ1 = 1.875104.

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

63

Taking k = 8000, h = 0.04, t = 0.001, we obtain the following figures representing the error between the exact and the numerical solution computed with some of the previous schemes.

5.2. Numerical simulations Here, we describe some of our numerical simulations, showing the types of behavior the solutions exhibit. For the simulations, we first consider that the motion of the left end is given by φ(t) = E sin (at), where a is the frequency of the driving device and f (x, t) = 0. First, we give some comparison results between the HHT-␣ scheme (4.13) and the two BDF2 schemes (4.5) and (4.11) . We also present some comparison between the SDIRK schemes (4.10), (4.12) and (4.14). Then, we compare hte HHT-␣ scheme (4.13) and the SDIRK scheme (4.14). Finally, we present other simulations and an example of a possible chaotic behavior. For each simulation, we give a table that contains all the parameters used for the computations and we present the noise characteristics of the system. We depict the (fast) Fourier transforms (FFT) of the motion of the beam’s end, u(1, t), that shows the frequency distribution of the vibrations. It is seen that higher

64

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

modes are excited by the contact and the system sometimes becomes rather ‘noisy’. The free vibration frequencies are given by 2 ωn = kβn+1 ,

n = 0, 1, . . .

where βn are ordinate solutions of cosh β cos β = −1. Thus, a straightforward calculation gives the first natural frequencies ω0 = 3.5148k (rad/s) and ω1 = 22.3480k (rad/s) √ √ When k = 300, ω0 = 60.878 (rad/s) and ω1 = 387.0787 (rad/s). When k = 8000, ω0 = 314.373 (rad/s) and ω1 = 1998, 865 (rad/s). 5.2.1. Comparison between the HHT scheme and the two BDF2 schemes For the use of the BDF2 schemes (4.5) and (4.11), we need to compute the values of V 0 and V 1 from the initial conditions u0 and v0 . The θ-method, with θ = 1/2, is used to compute Vi1 , for 0 ≤ i ≤ M. Then, we march in time using the scheme described in the previous section (Table 1). In Fig. 1a we present low driving frequency (a = 10) solutions obtained with the three methods. There is contact between the beam’s end and the stops over considerable periods of time (Fig. 1a). The results are similar for the displacement at the end of the beam, and there is little deflection at the stops (Fig. 1c). Note that the deflection is smaller with the first BDF scheme. The phase portraits show some differences for the velocity (Fig. 1b): the end oscillates rapidly during the free motion between the stops, but the phase portrait obtained with the HHT method is more noisy. Finally, the FFT depicted in Fig. 1d shows that the excited frequencies are odd multiples of the driving frequency a, i.e. 3a, 5a, . . . (Table 2). √ In Fig. 2, we consider the same simulation but with a lower value for the elastic modulus, i.e. k = 300: the results obtained for the displacement at the end of the beam are similar (Fig. 2a) and the phase portraits are irregular. At contact, there is no defection (Fig. 2c) and as can be shown the contact is more complex than before in the sense that the number of impact is higher. Moreover, the schemes give different results at contact (Fig. 2c): the results obtained by the two BDF schemes are near similar whereas the HHT scheme gives a (very) different results. Note carefully that this difference (at impact) is not very significant in the time interval [0, 2.85]: we know that the BDF and the HHT schemes introduce numerical damping which lead to period elongation and phase lag in the lower modes. This could explain the difference in the displacement we observed in Fig. 2c. The FFT, Fig. 2d, shows that odd multiples of the driving frequency are excited: 3a, 5a, . . . , with decaying amplitudes peaks (Table 3). We consider the same elastic modulus but with a smaller clearance. Here, the duration of free flights is very short comparing with the duration of contacts. The results obtained by both methods are similar with a small deflection at contact (Fig. 3c), but the FFT shows a wide band of excited frequencies (Table 4). Table 1 k √

8000

a

h

t

g1

g2

κ

E

α

10

0.04

10−4

−0.1

0.1

4 × 106

0.2

−0.1

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

65

Fig. 1.

Table 2 k √ 300

a

h

t

g1

g2

κ

E

α

10

0.04

10−4

−0.1

0.1

5 × 106

0.2

−0.1

Table 3 k √ 300

a

h

t

g1

g2

κ

E

α

10

0.02

2 × 10−5

−0.001

0.001

2 × 107

0.2

−0.1

66

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 2.

We now consider a “resonance case”: we choose a driving frequency near the first natural frequency of the beam. It is know that near an inner frequency, it is possible to obtain a chaotic behavior. For both methods the results look similar, but after some contacts (Fig. 4c), the computed displacement is different for all methods. The FFT shows that a large band of frequencies are excited, but the odd multiples of the driving frequency have higher amplitudes (Table 5).

Table 4 k √

300

a

h

t

g1

g2

κ

E

α

60

0.02

10−5

−0.001

0.001

3 × 107

0.2

−0.3

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

67

Fig. 3.

Finally, we perform the same computations with a bigger clearance. The phase portrait (Fig. 5b) indicates that the motion becomes noisy. Moreover, the displacement near the stops are very complicate and depends on the scheme used. Consider now the CPU-time of the previous computations: Table

BDF2 (first discretization) BDF2 (second discretization) HHT method

1

2

3

4

5

15.48 14.01 15.47

15.11 13.82 15.13

64.86 60.38 44.13

63.26 60.13 44.32

32.09 30.43 21.97

68

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 4.

Except for some computations, the three methods need the same CPU-time to perform the computations. Moreover, for some computations, the results are not similar, and thus, the choice of a scheme depends mainly of the problem. From my point of view, the use of the HHT method is preferable in that sense that we could control the numerical dissipation trough the parameter α. We do not have control on the dissipation with the BDF scheme. We compare now the previous results with the “modified” SDIRK method presented before and for which control of the damping is available. Table 5 k √

300

a

h

t

g1

g2

κ

E

α

60

0.02

2 × 10−5

−0.01

0.01

2 × 107

0.2

−0.3

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

69

Fig. 5.

5.3. Comparison between the SDIRK schemes (4.10), (4.12) and (4.14) From our knowledge, it is the first time that the SDIRK-OS scheme is used for such type of problems. In the following simulations, we only consider γ = 1: the schemes (4.10) and (4.12) give bad results when γ = 0.4359. We also consider the same simulations as in the previous section (Table 6). Table 6 k √ 300

a

h

t

g1

g2

κ

E

γ

10

0.04

10−4

−0.1

0.1

5 × 106

0.2

1

70

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 6.

The results obtained by both methods are almost similar (Fig. 6a), except that the phase portrait obtained by the third spatial discretization is smoother (Fig. 6b). The FFT shows again that the odd multiples of a are excited, in particular the first four ones (Table 7). The displacements are similar and there is only a small difference at impact (Fig. 7c). There is a little deflection at impact and the solutions “oscillate” which make the appearance of numerical frequencies. Thus a wide band of frequencies are excited (Fig. 7d and Table 8).

Table 7 k √

300

a

h

t

g1

g2

κ

E

γ

10

0.04

5 × 10−5

−0.001

0.001

9 × 106

0.2

1

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

71

Fig. 7.

Here, we consider again a driving frequency near the first inner mode of the beam. The displacements are irregular (Fig. 8a and b) and at impact the numerical solutions are clearly different (Fig. 8c). In fact the solutions obtained with the schemes (4.12) and (4.14) are exactly similar. We suggest that the scheme (4.10) is very sensitiv to the damping parameter. A large band of frequencies are excited and the odd mutiples of a present higher amplitudes (Table 9). We consider the same “physical” parameters with a bigger clearance. A difference between the numerical solutions of (4.12) and (4.14) (which are again similar) with (4.10) appears clearly (Fig. 9). We Table 8 k √ 300

a

h

t

g1

g2

κ

E

γ

60

0.02

10−5

−0.001

0.001

3 × 107

0.2

1

72

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 8.

suggest that this (important) difference depends both on the space discretization, the damping parameter and the clearance. Finally, it seems that at “resonance”, the numerical solutions depend strongly on the space discretization. Consider now the CPU-time of the previous computations: Table

SDIRK(4.10) SDIRK(4.12) SDIRK(4.14)

6

7

8

9

73.81 123.86 42.011

140.47 88.31 50.69

556.98 1018.99 304.59

559.12 1195.82 259.28

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

73

Table 9 k √ 300

a

h

t

g1

g2

κ

E

γ

60

0.02

2 × 10−5

−0.01

0.01

2 × 107

0.2

1

It can be shown that the CPU-time depends on the SDIRK schemes we used and thus depends on the space discretization. 5.4. Comparison between the HHT scheme (4.13) and the SDIRK scheme (4.14) We now compare the schemes (4.13) and (4.14), and we consider a new parameter: the number of impacts. The “impact” is also defined in the following way: there is an impact at time tn+1 if at time tn the

Fig. 9.

74

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Table 10 k √

300

a

h

t

g1

g2

κ

E

α

γ

10

0.04

10−4

−0.1

0.1

5 × 106

0.2

−0.1

0.4359

beam is not in contact with one of the stops. We also count the number of impacts of the beam against the stops showing that the use of the HHT scheme or the SDIRK scheme give very different results. In the following simulations, we consider σ = 0.4359 such that the SDIRK scheme becomes third order accurate (Table 10). These two methods give different results (Fig. 10c) at impacts. This difference comes both from the phase lag error, the period elongation introduced by both methods and the relatively significant

Fig. 10.

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

75

Fig. 11.

clearance. Note, the the FFTs indicate the same excited frequencies (see Fig. 11 and Table 11). The results are similar because the clearance is very small and, thus, the duration of free flights is very short between the two stops (Table 12). Though the frequency of the vibrating device increases, the results obtained by both methods are similar, except for some impacts (see Fig. 12c and Table 13).

Table 11 k √ 300

a

h

t

g1

g2

κ

E

α

γ

10

0.025

25 × 10−5

−0.001

0.001

2 × 107

0.2

−0.1

0.4359

76

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Table 12 k √

a 300

60

h 0.02

t 10

−5

g1 −0.001

g2

κ

0.001

6 × 10

7

E

α

γ

0.2

−0.1

0.4359

Fig. 12.

Table 13 k √

300

a

h

t

g1

g2

κ

E

α

γ

60

0.02

2 × 10−5

−0.01

0.01

4 × 107

0.2

−0.1

0.4359

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

77

Fig. 13.

The clearance is larger, the free flight duration too. Thus, the phase lag error and the period elongation become significant. These would explain the conclusion: the SDIRK and the HHT schemes give approximatively the same results. It is clear that the choice of the damping parameter is crucial and plays a very important role to control the undesired high frequencies. In other words, too much damping can avoid some impacts. Consider now the CPU-time (Fig. 13): Table

HHT method SDIRK method

10

11

12

13

15.59 43.24

51.74 192.93

43.27 157.60

21.80 119.37

78

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

The HHT-␣ is clearly the faster algorithm. In other words, the SDIRK-OS is of three order and present good stability properties; moreover it permits to implement embedded scheme. Let us now consider the following table which gives the number of impacts of the beam against the stops: Table

HHT method SDIRK method

10

11

12

13

3904 4494

291 266

688 824

2062 1130

The table explain, in part, the “noisy” phase-portrait we obtain with the HHT method. In other words, we show that the number of impacts could depend on the choice of the method and the damping parameter. This is of great importance if we want to evaluate the impact forces and thus the damage process. 5.5. Other simulations We now consider other interesting cases computed using the HHT-␣ scheme. 5.5.1. One stop We first consider one of the above simulation but with only one stop. The results and in particular the phase portrait show that the motion is periodic. Moreover, the FFT shows us that multiples of the driven frequency a are excited, but the amplitudes of the even multiples are much higher than the odd ones. In particular, the frequency 6a, which is near the first mode of the cantilever beam, has high amplitude because the duration of the free flight of the beam is longer with respect to the duration of contact (see Table 14 and Fig. 14). 5.5.2. Simulation with a loaded force We now consider that the supporting device is not moving and we take the following right-hand side [14] f (x, t) = E cos (wt)δ(x − xd ), where δ(x − xd )

=1

if x = xd

=0

if x = xd

We use the same parameters as in [14] (see Table 15). Here, the force acts only at xd = 0.23. In Fig. 15a, it can be seen that the contact time is relatively long in comparison with the free flight. The motion is very irregular, Fig. 15b. The FFT show that only some frequencies are excited. Finally, in Fig. 15c, we Table 14 k √

300

a

h

t

g1

κ

E

α

10

0.04

5 × 10−5

−0.001

7 × 106

0.2

−0.01

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

79

Fig. 14.

show the shear stress at the end of the beam, i.e. σ (1, t), with respect to the time, which indicates that the motion is periodic (Table 16). We now consider that the force acts at xd = 0.7: the contact time is longer and the displacement is more regular, Fig. 16a and b. In Fig. 16c, we show the shear stress at the end of the beam, i.e. σ (1, t),

Table 15 k2

h

t

g1

g2

κ

w

xd

E

α

282.84

0.04

5 × 10−5

−0.001

0.001

3 × 107

6

0.23

67.024

−0.1

80

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 15.

with respect to the time. Finally, the FFT shows that multiples of w are excited. The motion seems to be periodic. Remark 13. The computations were done on a PIII-800 Mhz Personal Computer, using Matlab for computing and graphing.

Table 16 k2

h

t

g1

g2

κ

w

xd

E

α

282.84

0.04

5 × 10−5

−0.001

0.001

3 × 107

6

0.70

67.024

−0.1

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

81

Fig. 16.

5.5.3. Chaotic behavior It is well known that impact systems could exhibit some chaotic behavior for a wide variety of parameters. We also present a (possible) chaotic behavior obtained with one of the previous models. We use the HHT-␣ algorithm (Table 17). The displacement may be considered as chaotic, because the phase portrait (u(1, t +2π/w), u (1, t +2π/w)) in Fig. 17b and the stress, σ (1, t), in Fig. 17c are very irregular. Moreover, in Fig. 17d, we show that a wide band of frequencies is excited. Table 17 k √ 300

w

h

t

g1

g2

κ

E

α

60

0.025

2 × 10−5

−0.1

0.1

5 × 107

0.2

−0.1

82

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

Fig. 17.

6. Conclusions We compare different algorithms used to simulate the vibrations of an elastic beam between (flexible) stops. These algorithms are obtained using the method of lines. We perform and present some simulations. The HHT-␣ scheme is clearly the most interesting from different point of view: it is fast, simple to implement and it permits to damp out spurious high frequencies. The numerical simulations exhibit some complicated behavior. It seems that the stops introduce a varied and complicated type of oscillations, some of them seem to be irregular and might be chaotic. Moreover, the FFT of the solutions shows that multiples of the driving frequency are excited in the process. Sometimes only the odd or even multiples, in other cases all multiples appear with different amplitudes. At high driving frequencies or for a low stiffness, it is seen that a wide range of frequencies

Y. Dumont / Mathematics and Computers in Simulation 60 (2002) 45–83

83

is excited, indicating a complicated motion. Finally, we show that the number of impacts of the beam against the stops could depend on the shemes we used and even on the damping parameter. This could be of great importance if we want to give estimates of the contact forces, which could be used to evaluate the “fatigue-life” of the system. The next step is also to perform a numerical analysis of the previous schemes, and if possible, to compare the numerical simulations to experimental results. Another interesting study would be to consider various positions of the stops along the beam or to consider that the stops are moving on a slider [5]. Finally, a numerical comparison of the normal compliance condition with the Signorini condition coupled with the Newton’s impact law [4,14] could be of great interest. References [1] K.T. Andrews, M. Shillor, S. Wright, On the dynamic vibrations of an elastic beam in frictional contact with a rigid obstacle, J. Elasticity 42 (1996) 1–30. [2] K.T. Andrews, K.L. Kuttler, M. Shillor, On the dynamic behaviour of a thermoviscoelastic body in frictional contact with a rigid obstacle, Eur. J. Appl. Math. 8 (1997) 417–436. [3] Y. Dumont, D. Goeleven, M. Rochdi, M. Shillor, Simulations of Beam’s vibrations between stops, Preprint 99. [4] Y. Dumont, Some remarks on a vibro-impact scheme, submitted for publication. [5] Y. Dumont, K. L. Kuttler, M. Shillor, Analysis and simulations of vibrations of a beam with a slider, submitted for publication. [6] H.M. Hilber, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Eng. Struct. Dyn. 5, 1977, 283-292 [7] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity, SIAM, Philadelphia, 1988. [8] A. Klarbring, A. Mikelic, M. Shillor, Frictional contact problems with normal compliance, Int. J. Eng. Sci. 26 (8) (1988) 811–832. [9] K.L. Kuttler, M. Shillor, Vibrations of a beam between two stops, Dyn. Contin. Discrete Impuls. Syst. Ser. B: Appl. Algorithms 8 (1) (2001) 93–110. [10] J.A.C. Martins, J.T. Oden, A numerical analysis of a class of problems in elastodynamics with friction, Comput. Methods Appl. Mech. Eng. 40 (1983) 327–360. [11] F.C. Moon, Chaotic and Fractal Dynamics, Wiley, New York, 1992. [12] F.C. Moon, S.W. Shaw, Chaotic vibration of a beam with nonlinear boundary conditions, J. Nonlinear Mech. 18 (1983) 465–477. [13] B. Owren, H.H. Simonsen, Alternative integration methods for problems in structural dynamics, Comput. Methods Appl. Mech. Eng. 122 (1995) 1–10. [14] L. Paoli, Analyse numérique de vibrations avec contraintes unilatérales, Ph.D. thesis, Université Claude Bernard, Lyon, France, 1993. [15] K.L. Kuttler, A. Park, M. Shillor, W. Zhang, Unilateral dynamic contact of two beams, Math. Comput. Modelling 34 (3/4) (2001) 365–384.