A viscoelastic Timoshenko beam with Coulomb law of friction

Applied Mathematics and Computation 218 (2012) 7078–7099

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A viscoelastic Timoshenko beam with Coulomb law of friction Jeongho Ahn Department of Mathematics and Statistics, Arkansas State University, P.O. Box 70, State University, AR 72467, United States

a r t i c l e

i n f o

Keywords: Viscoelastic Timoshenko beams Signorini contact conditions Coulomb law of friction Energy balance

a b s t r a c t This work focuses on the analysis and numerical simulations of dynamic frictional contact between a Timoshenko beam and a stationary rigid obstacle. The beam is assumed to be viscoelastic, and is clamped at its left end while it is free at its right end. When the beam moves horizontally and its right end contacts the rigid obstacle, contact forces arise and then the vertical motion of its right end is considered with friction. Thus the right end of the beam satisfies two contact conditions: the Signorini unilateral contact condition and Coulomb law of dry friction. Moreover, the slip rate dependence of the coefficient of friction is taken into account. The existence of weak solutions to the problem is shown by using finite time stepping and the necessary a priori estimates that allow us to vanish the time step in the limiting process. The energy balance in the system is investigated both theoretically and numerically. A fully discrete numerical scheme for the variational formulation of the problem is implemented and numerical results are presented. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction In this work, we consider a system which models dynamic frictional contact of a long slab. The system consists of a viscoelastic (Kelvin-Voigt type) rod with the Signorini contact conditions for the horizontal displacement and a viscoelastic Timoshenko beam with the classical Coulomb friction law for the vertical displacement. Our aim is to extend the results in the literature concerning beams in frictional contact to more realistic setting. Both the rod and the beam which have their left ends clamped and their right ends free are assumed to be deformed horizontally before their contact and vertically after contact. When the rod is stretched out by applying tensile traction to its right end, it may come in contact with a stationary rigid obstacle. Once there are contact forces, the viscoelastic Timoshenko is activated and starts moving upward or downward. Thus, the Timoshenko beam which is described by two couple partial differential equations is in frictional contact at its right end. The similar physical setting can be found in the papers [6,7,19]. Note that general d-dimensional viscoelastic bodies (1 6 d 6 3 in applications) have been considered in [19]. Since the horizontal motion of the rod has already been studied in the papers [14,15,26] for the purely elastic case and in the paper [17] for the viscoelastic case, we mostly deal with the frictional contact problem of the viscoelastic Timoshenko beam. Andrew et al. [7] considered the horizontal motion of viscoelastic rods and the vertical motion of viscoelastic Euler–Bernoulli beams with two couples of contact conditions; one is Signorini contact conditions and Coulomb law of friction and the other is normal compliance conditions and Coulomb law of friction. Indeed, they focused on showing the existence of solutions for the vertical displacement of Euler–Bernoulli beams with frictional contact, based on the theory of variational inequalities. In the paper [6], dynamic frictionless contact of viscoelastic Timoshenko beams has been studied. By adding frictional effects to the model in [6], we use the time discretization to prove the existence of weak solutions. Unlike frictionless contact E-mail address: [email protected] URL: http://www.mathstat.astate.edu/~jahn/ 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.12.069

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7079

problems, differential operators for frictional contact problems cannot be self-adjoint, which causes one of the main mathematical difficulties. This implies that improving the regularities of solutions for frictional contact problems is considerably harder than for frictionless contact problems. Investigating an energy balance (or conservation of energy) is an important issue in dynamic contact problems. The previous work on energy balances has conducted for frictionless contact (e.g., see the papers [1,2,5,6]). It was proved that energy balances for boundary thin obstacle problems (e.g., see [5,6]) and the thick obstacle problems with even adhesion (e.g., see [2]) under the assumption that contact forces are bounded in a suitable space. Jarušek and Eck [13] showed the existence of solutions, assuming that the coefficient of friction is sufficiently small but can be dependent of solutions. Cocou [9] proved the existence of solutions with the constant coefficient of friction, using a nonlocal version of the Coulomb friction law in viscoelasticity and regularizing the normal contact forces. We note that there are many ways to regularize contact forces; for example, one can use the convolution to average the displacement and the velocity, which enables one to take an average operator on the normal contact forces (e.g., see [18]). Extending the results of [9], Kuttler and Shillor [19] introduced a coefficient of friction which is dependent of the relative slip between a viscoelastic body and a moving rigid (or deformable) foundation. This case is applicable to many physical settings. In order to show a possible energy balance with Coulomb law for friction, the friction coefficient is assume to be time dependent but given. Thus, in Section 6 we assume that it is a function of the slip rate of the right end. While the important one-dimensional problem with slip rate dependent coefficient of friction has been studied by Inonescu and Paumier [12] and Paumier and Renard [22], the energy balance was not investigated in these papers. Applying the result of Stewart’s paper [27] and imposing the slip rate dependent friction, we can prove an energy balance theoretically. Our numerical results support the energy balance as well. Numerical methods for frictional contact problems over d-dimensional domain with 1 6 d 6 3 in the quasi or static case were studied in the paper [24]. Numerical methods for the dynamic one-dimensional problem or similar physical setting to our frictional problem are presented in the papers [10,16], using the finite difference method over both the time space and the spacial domain; especially, the numerical algorithm of [16] is based on a Crank–Nicolson scheme. Unlike those papers, we will use the time discretization and the finite element method on the spacial domain. Our numerical schemes are useful to support numerical evidences of energy balance. The remaining sections of this paper are organized as follows. Section 2 illustrates our dynamic frictional contact problem based on the physical reference configuration. Section 3 provides the mathematical notation and background. In Section 4 we derive a variational formulation which plays an important role in establishing numerical formulations in the time space [0, T]. In Section 5 we set up the numerical formulations, employing the midpoint rule. In Section 6, a possible energy balance with Coulomb friction law is discussed. In Section 7 the fully numerical schemes with Galerkin approximations are introduced. In order to compute the next step solutions of a nonlinear equation at each time step, we apply Newton’s method. In Section 8, numerical results (simulations) are presented. We conclude this paper in Section 9. 2. Problem statement We begin by describing the setting of the problem for the viscoelastic rod; in the reference configuration it has the initial length l, and its left end is fixed, and as a result of the horizontal traction fh it stretches horizontally. When its right end touches a stationary rigid obstacle such as a wall, which is situated at x = u, the beam moves upward or downward, producing tangential stresses (shear stresses) between its right end and the rigid obstacle. For the sake of a simplicity, we normalize the material density in the equation of motion. Then the horizontal displacement v = v(t, x) of the rod satisfies the following third order partial differential equation

v€ ¼ c1 v xx þ a1 v_ xx þ fh ðt; xÞ

in ð0; T  ð0; lÞ;

where T > 0 is a fixed end time and fh denotes applied forces in the horizontal direction and c1 and a1 are the coefficient of elasticity and viscosity, respectively. Here the subscripts denote the derivatives with respect to the subscribed variable x and (_) and (__) denote the first and second derivative with respect to time t, respectively. Since the left end of the rod is fixed and its right end is free, the essential boundary condition v(t, 0) = 0 is imposed at the left end x = 0 and the natural boundary condition rN ðtÞ ¼ c1 v x ðt; lÞ þ a1 v_ ðt; lÞ at its right end x = l, where rN are normal contact forces (pressures). When the beam moves horizontally and touches a wall but does not penetrate the rigid obstacle, normal contact forces arise and rN(t) < 0 and thus the Signorini contact conditions are applied at x = l. Note that the wall is located on the right of the end of the beam, i.e., x = u > l. Thus the horizontal motion of the viscoelastic rod with the Signorini contact conditions is established in the strong sense:

v€ ¼ c1 v xx þ a1 v_ xx þ fh ðt; xÞ in ð0; T  ð0; lÞ; v ðt; 0Þ ¼ 0 on ð0; T; rN ðtÞ ¼ c1 v x ðt; lÞ þ a1 v_ x ðt; lÞ on ð0; T; 0 P rN ðtÞ ? v ðt; lÞ  u 6 0 on ð0; T; v ð0; xÞ ¼ v 0 ðxÞ; v_ ð0; xÞ ¼ v_ 0 ðxÞ in ð0; lÞ;

ð2:1Þ ð2:2Þ ð2:3Þ ð2:4Þ ð2:5Þ

7080

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

where two initial data are given in (2.5). Indeed, the Signorini contact conditions are written by the complementarity conditions (2.4). The meaning of the notation 0 P a \ b 6 0 is that both are nonpositive and either a or b is zero. In more general, for vectors a, b the complementarity conditions 0 P a \ b 6 0 means that a, b 6 0 componentwise and the dot product a  b = 0. In this paper, bold face characters are used to denote vectors or vector-valued functions. According to our physical reference configuration (see Fig. 2.1) two quantities rN and v(t, l)  u are nonpositive in (2.4) unlike the normal complementarity conditions 0 6 a \ b P 0. The variational formulation which is equivalent to the PDE systems (2.1)–(2.5) will be provided in Section 4, in order to consider the existence of solutions. We now turn into the vertical displacement of the viscoelastic Timoshenko beam which is described by the following two coupled partial differential equations:

_ in ð0; T  ð0; lÞ; qI€h ¼ EIhxx þ KAGðux  hÞ þ a2 EIh_ xx þ a2 KAGðu_ x  hÞ _ qAu€ ¼ KAGðuxx  hx Þ þ a2 KAGðu_ x  hx Þ þ fv ðt; xÞ in ð0; T  ð0; lÞ;

ð2:6Þ ð2:7Þ

where fv denotes the vertical applied force, and a2 is the coefficient of viscosity and u = u(t, x) and h = h(t, x) represent the vertical displacement and the rotatory angle of the beam, respectively. The Timoshenko beam is assumed to be uniform. Here q is the density of the beam and G the modulus of elasticity of shear, and K the geometric dependent distribution of shear stress. EI denotes the flexural rigidity of the beam, where I is the second moment of inertia and E Young modulus. From now on, we shall put all coefficients q = A = K = G = E = I = 1 for our convenience. The Timoshenko beam is also assumed to be clamped at x = 0 and be free at x = l. Therefore the essential boundary conditions of the Timoshenko beam imposed at the left end x = 0 are u(t, 0) = h(t, 0) = 0 and the natural boundary conditions at its right end x = l are _ lÞ  u_ x ðt; lÞÞ for t 2 (0, T], where rT is the tangential stress. In order hx ðt; lÞ þ a2 h_ x ðt; lÞ ¼ 0 and rT ðtÞ ¼ hðt; lÞ  ux ðt; lÞ þ a2 ðhðt; to see how the tangential stresses rT work at x = l, we introduce the Coulomb law of dry friction:

jrT ðtÞj 6 lðtÞjrN ðtÞj on ð0; T; 

ð2:8Þ

_ lÞ ¼ krT ðtÞ for some k P 0; jrT ðtÞj ¼ lðtÞjrN ðtÞj ) uðt; _ lÞ ¼ 0; jrT ðtÞj < lðtÞjrN ðtÞj ) uðt;

ð2:9Þ

where l() P 0 is called the coefficient of friction. In this paper, we consider the coefficient l of friction as a function in terms of time t, unlike a positive constant l considered in the paper [7]. The Coulomb friction law (2.8) and (2.9) can be understood from the physical point of view; if the magnitude of tangential stress is the same as the coefficient of friction times the magnitude of normal contact force, i.e., l()jrN()j, the tangential stress will be in the opposite direction to the sliding velocity and if the magnitude of the tangential stress is strictly than l()jrN()j, the right end of the beam will stick to the obstacle. We note that there are other ways to formulate Coulomb law for frictional contact; for example, maximum dissipation principle (see Erdmann’s paper [11]). Thus we are led to consider the following equations of the viscoelastic Timoshenko beam with Coulomb’s law of dry friction:

€h ¼ hxx þ ux  h þ a2 ðh_ xx þ u_ x  hÞ _ in ð0; T  ð0; lÞ; _ _ € u ¼ ðuxx  hx Þ þ a2 ðux  hx Þ þ fv ðt; xÞ in ð0; T  ð0; lÞ;

ð2:11Þ

0 ¼ uðt; 0Þ ¼ hðt; 0Þ on ð0; T; 0 ¼ hx ðt; lÞ þ a2 h_ x ðt; lÞ on ð0; T;

ð2:13Þ

_ lÞ  u_ x ðt; lÞÞ on ð0; T; rT ðtÞ ¼ hðt; lÞ  ux ðt; lÞ þ a2 ðhðt; jrT ðtÞj 6 lðtÞjrN ðtÞj on ð0; T;

ð2:15Þ

ð2:10Þ ð2:12Þ

Fig. 2.1. Frictional contact of a viscoelastic beam.

ð2:14Þ

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099



_ lÞ ¼ krT ðtÞ for some k P 0; jrT ðtÞj ¼ lðtÞjrN ðtÞj ) uðt; _ lÞ ¼ 0; jrT ðtÞj < lðtÞjrN ðtÞj ) uðt;

hð0; xÞ ¼ h0 ðxÞ;

7081

ð2:16Þ

_ xÞ ¼ h_ 0 ðxÞ in ð0; lÞ; hð0;

ð2:17Þ

_ uð0; xÞ ¼ u_ 0 ðxÞ in ð0; lÞ;

ð2:18Þ

0

uð0; xÞ ¼ u ðxÞ;

where the initial data are given in (2.17) and (2.18). The mathematical analysis of this PDEs (2.10)–(2.18) will be also carried out through Sections 4 and 5. 3. Preliminaries and notations In this Section, the basic mathematical backgrounds and notations are presented to consider the existence of solutions. We introduce the Sobolev space H1 ð0; l; R2 Þ which consists of the vector-valued functions u from [0, l] to R2 . If u ¼ ðu; hÞ 2 H1 ð0; l; R2 Þ, its norm is defined as

Z

kukH1 ð0;l;R2 Þ ¼

!1=2

l

2

2

jux j þ juj dx

;

0

where juj2 = jhj2 + juj2 and juxj2 = jhxj2 + juxj2. If u 2 L2 ð0; l; R2 Þ, its norm is defined to be

Z

kukL2 ð0;l;R2 Þ ¼

!1=2

l

juj2 dx

:

0

Based on the Gelfand triples (see [29]):

V  H ¼ H  V  ; let w = (#, w)T 2 V and V ¼ H1cf ð0; l; R2 Þ and H ¼ L2 ð0; l; R2 Þ in which the subscripts ‘‘c’’ and ‘‘f’’ stand for ‘‘clamped’’ and ‘‘free’’, respectively. A duality pairing between X and X⁄ is denoted by h; iX  X , where X is a Banach space and X⁄ is its dual space. For our convenient notations, an inner product (, )H on space H will be used instead of (, )HH and if dual spaces are clear, h, i will be used instead of h; iX  X . The duality pairing between V and V⁄ is defined to be

h-; wi ¼

Z

l

ðm# þ -wÞdx;

0

_ Þ; uð; _ xÞÞT 2 V we define a linear where - = (m, -)T 2 V⁄. For the vector valued functions u = (h(t, ), u(t, ))T 2 V and u_ ¼ ðhðt;   elasticity operator Le : V ! V and a linear viscosity operator Lv : V ! V with appropriate boundary conditions at x = l as

hLe u; wi :¼

Z

l

½hx #x þ ðux  hÞðwx  #Þdx 0

_ wi :¼ a2 hLv u;

Z

l

and

_ ½h_ x #x þ ðu_ x  hÞðw x  #Þdx:

0

We note that hLe u; wi ¼ ðu; wÞV and hLe u; ui ¼ kuk2V . If the functions u satisfy the following inequality kuðtÞ  uðsÞkH1 ð0;lÞ 6 Cjt  sjp for some constant C > 0 with s – t 2 [0, T], the functions u are called to be Hölder continuous with  exponent 0< p 6 1. The Hölder spaces that contain all Hölder continuous functions u : ½0; T ! H1L ð0; lÞ are denoted by C p 0; T; H1L ð0; lÞ , where H1L ð0; lÞ ¼ fu 2 H1 ð0; lÞjuð0Þ ¼ 0g and ‘‘L’’ stands for the left end of beams. For any real numbers r Sobolev spaces Hr ðRÞ with the fractional index can be defined, using the Fourier transform (e.g., see [23, Section 5.1] or [28, Section 4.1]);

^ 2 L2 ðRÞg; Hr ðRÞ ¼ fu 2 S  ðRÞjhnir u

ð3:1Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where SðRÞ is the space containing all smooth functions which decrease rapidly at infinity and hni ¼ 1 þ jnj2 . S  ðRÞ is the a linear mapping from SðRÞ to R and is called the space of tempered distributions. Here the Fourier transform of u 2 L1 ðRÞ is defined by

^ ðnÞ ¼ F ½uðnÞ :¼ u

Z

1 1=2

ð2pÞ

1

uðsÞeins ds

1

and the inverse of the Fourier transform of u is also defined by

^ ðsÞ :¼ uðsÞ ¼ F 1 ½u

Z

1 ð2pÞ

1=2

1

1

^ ðnÞeins dn: u

7082

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Then it is well known (see [28, Section 4.1]) that

ðu; v ÞHr ðRÞ ¼ ðKr u; Kr v ÞL2 ðRÞ ; r

ð3:2Þ r



1

r

^ Þ. where the operator K on S ðRÞ is defined by K u ¼ F ðhni u Let z(t, ) = (z1(t, ), z2(t, ))T 2 V with all t 2 [0, T]. Since Coulomb law for friction is applied to only Eq. (2.11) of motion for the displacement u, we consider the restriction operator b : V ! H1L ð0; lÞ with bz = z2. Using the restriction operator b, we define the friction functional J : V ! Rþ for almost all t 2 [0, T] to be

J ðzðtÞÞ :¼ J ðbzðt; ÞÞ ¼ lðtÞjrN ðtÞjjz2 ðt; lÞj:

ð3:3Þ

As we shall see in next Sections 4 and 5, l; rN 2 Mð½0; TÞ, where Mð½0; TÞ is a space of Borel measures. We also note that J can be identified with a nonnegative Borel measure on [0, T] by

J ðzðBÞÞ ¼

Z

J ðzðtÞÞdt;

B

where B is any Borel set in the time interval [0, T]. Similarly, the action for the duality pairing hJ ðzðtÞÞ; wi will be taken by hJ ðbzÞ; bwi. Therefore, it is easy to see that for almost all t 2 [0, T]

hJ ðzðtÞÞ; wi :¼ hJ ðbzÞ; bwi ¼ lðtÞjrN ðtÞkz2 ðt; lÞjwðlÞ:

ð3:4Þ

Since J is not differentiable in terms of z in the strong sense, we use the regularization function J n to approximate J : for each integer n P 1 the functional J n is defined by

hJ n ðzðtÞÞ; wi ¼ lðtÞjrN ðtÞjwn ðz2 ðt; lÞÞwðlÞ; where the approximation of jsj is chosen by wn ðsÞ ¼ given by

 0  lðtÞjrN ðtÞjz2 ðt; lÞwðlÞ J n ðzðtÞÞ; w ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz2 ðt; lÞÞ2 þ ð1=nÞ2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ ð1=nÞ2 with s 2 R. Thus the directional (Gateux) derivative of J n is

for almost all t 2 ½0; T;

ð3:5Þ

where (0 ) is interpreted as a directional (Gateaux) derivative. So one can see that Coulomb law (2.15) and (2.16) for the friction can be handled by J 0n with sufficiently large n P 1. We note that there are many functions to regularize the friction functional J ; for example, the different regularized function is used in [8]. Letting f = (0, fv)T, we define the energy functional of the viscoelastic Timoshenko beams to be

_ EðtÞ ¼ E½uðtÞ; uðtÞ :¼

1 _ uÞ _ H þ hLe u; uiV  V Þ  hf; ui; ððu; 2

ð3:6Þ

_ ui _ is called the kinetic energy and hLe u; ui the elastic energy and hf, ui the potential energy. Note that (, )H can be where hu; extended into h, i on V⁄  V. It shall be seen in Section 6 that energy dissipates due to not only viscosity but also Coulomb friction. 4. The results for existence of solutions We recall the frictionless contact problem of the viscoelastic rod (2.1)–(2.5) which describes its horizontal motion. Then, the weak solutions v : ½0; T ! H1L ð0; lÞ that we want to seek satisfy the following conditions:

v€ ðtÞ ¼ Av ðtÞ  Bv_ ðtÞ þ f h þ rN ðt; xÞ; 0 P rN ðtÞ ? v ðt; lÞ  u 6 0 on ð0; T; v ð0Þ ¼ v 0 ; v_ ð0Þ ¼ v_ 0 ; H1L ð0; lÞ

ð4:1Þ ð4:2Þ ð4:3Þ

1

where A; B : ! H ð0; lÞ are self-adjoint operators and rN(t, x) = rN(t)d(x  l) with the Dirac delta function d. When we prove the existence of solutions, the contact force rN has to be understood in the sense of measures. Since this contact problem (4.1)–(4.3) has already been studied in several important papers [14,15,26], we shall skip proving the existence of solutions. One can find that the idea of the paper [4,6] is able to be applied to show the existence of solutions. Therefore in this work, our main consideration will be taken into the frictional viscoelastic problem of a Timoshenko beam (2.10)–(2.18). The contact forces rN form a connection between the weak formulations (4.1) and (4.9). Based on the assumption that there are solutions (v, rN) satisfying Eqs. (4.1)–(4.3), we summarize the main results for the existence of solutions in Theorem 2. Its proof will be presented in the next Section 5. First of all, we will see in the next Lemma 1 that the dynamic frictional contact problem (2.10)–(2.16) is equivalent to the following regularized variational problem: for sufficiently large n P 1

€ n ðtÞ þ Le un ðtÞ þ Lv u_ n ðtÞ þ J 0n ðu_ n ðtÞÞ  fðtÞ; wiV  V ¼ 0 8w 2 V: un 2 V : h u n

ð4:4Þ

We assume implicitly that u is the approximation of the solution u and for almost all t 2 [0, T] un ðtÞ; u_ n ðtÞ 2 V and € n ðtÞ 2 V  . One can be aware that a variational inequality equivalent to the regularized nonlinear equation (4.4) can be u derived easily.

7083

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Lemma 1. u = (h, u)T is a solution of the formulations 2.10, 2.11, 2.12, 2.13, 2.14, 2.15, 2.16 if and only if there is a sequence un such that un weakly converges to the solution u in V and satisfies the variational equality (4.4) with sufficiently large n P 1. Proof. Consider the equations of motion (2.10)–(2.16) as the following vector valued equation:

" # " # €h _ hxx þ ux  h þ a2 ðh_ xx þ u_ x  hÞ ¼ : € u uxx  hx þ a2 ðu_ x  h_ x Þ þ fv ðt; xÞ

ð4:5Þ

Suppose that the solution u of Eqs. (2.10)–(2.16) is sufficiently smooth. Then taking any w = (w, #)T 2 V and applying integration parts, it follows from (4.5) that

Z

l

€ wÞdx ¼  ð€h# þ u

Z

0

l

ðhx #x þ ðux  hÞðwx  #ÞÞdx þ

0

Z

l

fv wdx  a2

Z

0

l

_ ðh_ x #x þ ðu_ x  hÞðw x  #ÞÞdx  rT ðÞwðlÞ:

0

Thus we can have

€ ; wi ¼ hLe u  Lv u_ þ f; wi  rT ðtÞwðlÞ: hu n

_n

ð4:6Þ 

€n

⁄

_ u € Þ in V  V  V as n"1. In order to Now we choose the sequences ðu ; u ; u Þ 2 V  V  V are convergent weakly to ðu; u; approximate the tangential stresses rT which satisfy Coulomb law for friction (2.15) and (2.16), we apply (3.5):



 lðtÞjrN ðtÞju_ n ðt; lÞwðlÞ J 0n ðu_ n ðtÞÞ; w ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðu_ n ðt; lÞÞ2 þ ð1=nÞ2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi To see it, consider the first condition of (2.16). If ju_ n ðt; lÞj= ðu_ n ðt; lÞÞ2 þ ð1=nÞ2 ¼ 1 with sufficiently large n P 1, then we can qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ lÞ ¼ krT ðtÞ with k ¼ ðu_ n ðt; lÞÞ2 þ ð1=nÞ2 =ðlðtÞjrN ðtÞjÞ P 0. This k may be understood get only ju_ n ðt; lÞj P 1=n and thus uðt; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi in the sense of measures. For the second condition of (2.16) let ju_ n ðt; lÞj= ðu_ n ðt; lÞÞ2 þ ð1=nÞ2 < 1 with sufficiently large n P 1. _ lÞ ¼ 0. Indeed, we will see that u_ n ðt; lÞ converges to uðt; _ lÞ pointwise almost everyThen we get ju_ n ðt; lÞj 6 1=n and thus uðt; where in the next Section 5. Therefore we obtain the regularized variational problem (4.4). Conversely, we suppose that there is a sequence un with sufficiently large n P 1 such that un satisfies the variational formulation (4.4) and un converges weakly to u in V. Choosing all w 2 V such that w(l) = 0 and applying integration by parts with the essential boundary conditions (2.12), we can get the following variational formulation: for all w = (#, w)T 2 V

 n  € þ Le un þ Lv u_ n þ J 0n ðu_ n ðtÞÞ  f; w ¼ 0¼ u

Z

l

€h#dx þ

Z

l

€ wdx  u

Z

l

ðhx #x þ ðux  hÞðwx  #ÞÞdx 0 0 " # Z l Z l Z l € _ # h  ðhxx þ ux  h þ a2 ðh_ xx þ u_ x  hÞÞ _ _ fv wdx ¼ dx:  a2 ðhx #x þ ðu_ x  hÞðwx  #ÞÞdx þ _ € _ w u  ðuxx  hx þ a2 ðux  hx Þ þ fv ðt; xÞÞ 0 0 0 0

Notice that an approximation un is replaced by u in the second line of the previous equation. Thus the solution u solves the equations of motion (2.10) and (2.11). We also need to show that the natural boundary conditions (2.13) and (2.14) and Coulomb friction conditions (2.15) and (2.16) hold. Again, we choose sufficiently large n P 1 to get the weak convergence of the € n Þ 2 V  V  V  . Since J n satisfies the friction conditions (2.15) and (2.16) for the large n, for all sequences ðun ; u_ n ; u T w = (#, w) 2 V we can obtain

 n  lðtÞjrN ðtÞju_ n ðt; lÞwðlÞ € þ Le un þ Lv u_ n þ J 0n ðu_ n ðtÞÞ  f; w ¼ hu € n þ Le un þ Lv u_ n  f; wiV 0 V þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0¼ u ðu_ n ðt; lÞÞ2 þ ð1=nÞ2 € þ Le u þ Lv u_  f; wiV 0 V þ rT ðtÞwðlÞ ¼ hu

ð4:7Þ

and rT satisfies the conditions (2.15) and (2.16). Similarly using integration by parts, it follows that

€ þ Le u þ Lv u_  f; wiV 0 V hu

# Z l" € _ # h  ðhxx þ ux  h þ a2 ðh_ xx þ u_ x  hÞÞ dx þ ðux ðt; lÞ  hðt; lÞ ¼ _ € _ w u  ðuxx  hx þ a2 ðux  hx Þ þ fv ðt; xÞÞ 0 _ lÞÞÞwðlÞ þ ðhx ðt; lÞ þ a2 h_ x ðt; lÞÞwðlÞ: þ a2 ðu_ x ðt; lÞ  hðt;

Therefore from (4.7) and (4.8) we have

 n  € þ Le un þ Lv u_ n þ J 0n ðu_ n ðtÞÞ  f; w 0¼ u _ lÞÞ þ rT ðtÞÞwðlÞ þ ðhx ðt; lÞ þ a2 h_ x ðt; lÞÞwðlÞ: ¼ ðux ðt; lÞ  hðt; lÞ þ a2 ðu_ x ðt; lÞ  hðt; Thus taking w(x) = x, we can obtain

_ lÞÞ þ rT ðtÞ: 0 ¼ hx ðt; lÞ þ a2 h_ x ðt; lÞ ¼ ux ðt; lÞ  hðt; lÞ þ aðu_ x ðt; lÞ  hðt; The proof is complete. h

ð4:8Þ

7084

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Now we want to seek u : [0, T] ? V such that un with sufficiently large n P 1 is the approximation of the solution u and satisfies the following regularized nonlinear equations;

€ n ðtÞ ¼ Le un ðtÞ  Lv u_ n ðtÞ  J 0n ðu_ n ðtÞÞ þ f; u n

0

u ð0Þ ¼ u ;

ð4:9Þ

u_ n ð0Þ ¼ u_ 0 ;

ð4:10Þ n

_n

where Eq. (4.9) has to be considered in the sense of distributions and u (0) and u ð0Þ are the same as the initial data u(0) and _ uð0Þ. Theorem 2. Assume that the contact forces rN which satisfy Eqs. (4.1) and (4.3) exist in Mð0; T; H1 ð0; lÞÞ, the space of Borel measure on the time space [0, T] and the coefficient of friction, l for Coulomb law is given in Mð½0; TÞ and f 2 L1(0, T; L2(0, l)). If the initial data u0 ; h0 2 H1L ð0; lÞ and u_ 0 ; h_ 0 2 L2 ð0; lÞ, then there exist solutions (h, u) such that

  h 2 L1 0; T; H1L ð0; lÞ \ W 1;1 ð0; T; L2 ð0; lÞÞ \ C 1=2 ð0; T; L2 ð0; lÞÞ;   u 2 L1 0; T; H1L ð0; lÞ \ W 1;1 ð0; T; L2 ð0; lÞÞ \ C 1=2 ð0; T; L2 ð0; lÞÞ \ Cð½0; T  ð0; TÞ;   u_ 2 L2 0; T; H1L ð0; lÞ ; h_ 2 L2 ð0; T; H1L ð0; lÞÞ; € 2 L2 ð0; T; H1 ð0; lÞÞ; u

€h 2 L2 ð0; T; H1 ð0; lÞÞ:

5. Convergence with time discretization As we mention in Section 1, there are several papers to show the existence of solutions for that kind of the dynamic frictionless contact problem (2.1)–(2.5). Time discretization also enables us to approach the dynamic frictionless contact problem both theoretically and numerically. We start by dividing the time interval [0, T] into n subintervals, not necessarily of the uniform spacing. In our approach, we denote the initial time 0 by t0 and the end time T by tP so that

0 ¼ t0 < t 1 < t 2 <    < t n1 < t P ¼ T; where the time step size ht = tk  tk1 for 1 6 k 6 P and each time step tk = kht. For the sake of the simple notations the numerical approximations v_ ðtk ; xÞ, v(tk, x), rN(tk) at each time step tk are denoted by v_ k , vk, rkN , u_ k , and uk, respectively. Then the numerical formulation which corresponds to the continuous formulation (2.1)–(2.5) can be established as follows:

1 kþ1 1 1 ðv_  v_ kþ1 Þ ¼  Aðv kþ1 þ v k Þ  Bðv_ kþ1 þ v_ k Þ þ fhk þ rkN ; h 2 2 1 kþ1 1 kþ1 k k _ _ ðv þ v Þ ¼ ðv  v Þ; 2 ht kþ1 k  u ? rN 6 0: 0Pv

ð5:1Þ ð5:2Þ ð5:3Þ

In the numerical scheme, the midpoint rule is used for the elasticity and viscosity and the implicit Euler method is used for the Signorini contact conditions (complementarity conditions). Readers may refer the papers (e.g., [1,3,4,6]) for proving the convergence of numerical trajectories based on the numerical scheme. Therefore we will focus on showing the existence of solutions for the frictional contact problem (2.10)–(2.18). _ k ; xÞ and u(tk, x) at each time step tk are denoted by u_ k and uk, respectively. Using a Similarly the numerical solutions uðt small regularized parameter ht > 0 with ht = 1/n, for a given lk = l(tk) we define the discrete friction functional J ht at each time step tk to be

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 hJ 1=ht ðzk Þ; wi ¼ lk rkN ðzk Þ2 þ ht wðlÞ Thus we apply Eq. (4.9) to set up the following numerical formulations over the time space (0, T]:

 1  1  ht ;kþ1 1  1  ht ;kþ1 k u_  u_ ht ;k ¼  Le uht ;kþ1 þ uht ;k  Lv u_ ht ;kþ1 þ u_ ht ;k  J 01=ht þ u_ ht ;k þ f ; u_ ht 2 2 2

ð5:4Þ

1  ht ;kþ1 1  ht ;kþ1 u_ þ u_ ht ;k ¼ u  uht ;k ; 2 ht

ð5:5Þ

where uht ;k and u_ ht ;k are the discrete numerical solutions dependent of the sufficiently small parameter ht and the body force f is defined by the step function f(t) = fk for t 2 [tk, tk+1). Since there is no vertical motion of the viscoelastic Timoshenko beam before contact,

(5.4) and (5.5) will not be considered until there occurs contact forces. So we assume that there is a time step tk such that rkN P 0 and vk = 0 for any t > tk and uht ;t ¼ 0 for any t < tk. For our simple notations, uht ;k and u_ ht ;k will be denoted by uk and u_ k , respectively. The additional approximate equation sets up in (5.5). The main numerical scheme employed here is the midpoint rule which is used for the elasticity and viscosity, and the regularized friction functional J n . Making a connection between the numerical formulations (5.1)–(5.3) for the horizontal motion of the rod and the numerical

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7085

formulations (5.4) and (5.5) for the vertical motion of the Timoshenko beam, we construct the approximation of contact forces, ðrN Þht ðt; Þ and the approximation of the coefficient of friction, lht ðtÞ and such that ðrN Þht ðt; Þ ¼ rkN and lht ðtÞ ¼ lk for t 2 [tk, tk+1). Then we define their product to be

lht ðrN Þht ðt; xÞ ¼ ht

P1 X





lk rkN dðt  kht Þdðx  lÞ:

ð5:6Þ

k¼0

As we will see the following Lemma 3, the energy functional with the time discretization which plays an important role in showing the boundedness of numerical trajectories in suitable spaces is defined by

E k :¼ Eðt k Þ ¼

1 k k ðu_ ; u_ Þ þ hLuk ; uk i : 2

ð5:7Þ

The numerical trajectories are constructed as follows; we use the linear continuous interpolant uht ¼ ðhht ; uht ÞT to approxi_ Þ such that mate the displacement u(t, ) and the constant interpolant u_ ht ¼ ðh_ ht ; u_ ht ÞT to approximate the velocity uðt; _ Þ ¼ u_ k on the half open interval (tk1, tk]. uht ðtk1 ; Þ ¼ uk1 and uht ðtk ; Þ ¼ uk on the subinterval [tk1, tk] and uðt;  Lemma 3. Suppose that the numerical solutions v k ; v_ k ; rkN 2 H1L ð0; lÞ  L2 ð0; lÞ  H1 ð0; lÞ at each time step tk satisfy 5.1, 5.2, 5.3 and the numerical solutions ðuk ; u_ k Þ at each time step tk satisfy Eqs. (5.4) and (5.5) with sufficiently small ht > 0. If ðu0 ; u_ 0 Þ 2 V  H and f : [0, T] ? V0 is approximated by the step function f(t) = fk for t 2 [tk, tk+1), then ðuk ; u_ k Þ is uniformly bounded in V  H, independent of the time step size ht > 0. Proof. We use Eq. (5.5) to obtain the following; for a sufficiently small time step size ht > 0

 1 1 1 1 kþ1 k ðu_ ðu_ kþ1  u_ k ; u_ kþ1 þ u_ k ÞH ¼ hLe ðukþ1 þ uk Þ  Lv ðu_ kþ1 þ u_ k Þ; ukþ1  uk i þ hJ 01=ht þ u_ k Þ þ f ; ukþ1  uk i 2ht 2ht ht 2 1 1 ðukþ1 þ uk ; ukþ1  uk ÞV  ðu_ kþ1 þ u_ k ; u_ kþ1 þ u_ k ÞV ¼ 2ht 4  

 1 1 kþ1 k 0 k _ _ J 1=ht ðu þ u Þ þ f ; u_ kþ1 þ u_ k þ 2 2

From (3.5) we have

  1  kþ1 2 1  kþ1 2 1 lk jrkN jðu_ kþ1 ðlÞ þ u_ k ðlÞÞ2 1 k þ ðf ; u_ kþ1 þ u_ k ÞH ku_ kH þ kukþ1 k2V  ku_ kH þ kuk k2V ¼  ku_ kþ1 þ u_ k k2V  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ht 2ht 4 2 2 2 2 ðu_ kþ1 ðlÞ þ u_ k ðlÞÞ þ 4ht 6

1 k kþ1 þ u_ k ÞH : ðf ; u_ 2

Using (5.7), it follows that

E kþ1 þ

ht kþ1 ht lk jrkN jðu_ kþ1 ðlÞ þ u_ k ðlÞÞ2 ht k 6 E k þ ðf ; u_ kþ1 þ u_ k ÞH : þ u_ k k2V þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ku_ 4 2 2 2 2 kþ1 k ðu_ ðlÞ þ u_ ðlÞÞ þ 4ht

ð5:8Þ

When we apply the trapezoidal rule and the Cauchy–Schwartz inequality to (5.8), the telescoping series allows us to have

Ek þ



 Z tk k1 k1 li riN ðu_ iþ1 ðlÞ þ u_ i ðlÞÞ2 ht X ht X 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 E0 þ Tkfk2L1 ð0;T;HÞ þ ku_ iþ1 þ u_ i k2V þ ku_ ht ðsÞk2H ds : 2 4 i¼1 2 i¼1 2 0 ðu_ iþ1 ðlÞ þ u_ i ðlÞÞ2 þ 4ht

ð5:9Þ

When we apply Grownwall inequality to the last term in (5.9), we can see easily that maxkP1kukkV 6 M1 < 1 and maxkP1 ku_ k kH 6 M 2 < 1, independent of ht > 0. This proof is complete. h Thus by the previous Lemma 3 we can see that uht is uniformly bounded in C(0, T; V) and u_ ht is in L1(0, T; H), which implies that uht and hht are uniformly bounded in Cð0; T; H1L ð0; lÞÞ and u_ ht and h_ ht are in L1(0, T; L2(0, l)). We note that throughout this paper the constants C may be different in each occurrence. Lemma 4. Assume that the numerical  trajectory u_ ht satisfies the numerical formulations (5.4) and (5.5). Then h_ ht and u_ ht are uniformly bounded in L2 0; T; H1L ð0; lÞ for sufficiently small ht > 0. Proof. Replacing k by P from (5.9), we can get

C P ht

P1 X k¼1

ku_ kþ1 þ u_ k k2V ¼ ht

P1 P1 X X ðu_ kþ1 þ u_ k ; u_ kþ1 þ u_ k ÞV ¼ ht ðku_ kþ1 k2V þ ku_ k k2V þ 2ðu_ kþ1 ; u_ k ÞV Þ: k¼1

i¼1

ð5:10Þ

7086

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Now the piecewise constant interpolant u_ ht can be regarded as a linear combination of B splines of degree 0, i.e., R P u_ ht ¼ Pk¼1 u_ k Bk , where Bk = 1 for t 2 (tk1, tk] and Bk = 0 for t 2 [0, T]n(tk1, tk]. Note that for each 1 6 k 6 P, ðtk1 ;tk  ku_ k k2V dt ¼ R 2 k kþ1 k ht ku_ kV and ðtk1 ;tk  ðu_ ; u_ ÞV dt ¼ 0. Therefore it follows from (5.10) that for any sufficiently small ht > 0

Z P1 P1 X X C P ht ðku_ kþ1 k2V þ ku_ k k2V þ ðu_ kþ1 ; u_ k ÞV Þ ¼ k¼1

Z

¼

0

T

ku_ ht k2V dt þ

Z

Tht

ht

i¼1

!

ðtk ;t kþ1 

ku_ kþ1 k2V dt þ

!

Z ðt k1 ;tk 

ku_ k k2V dt

ku_ ht k2V dt :

Thus u_ ht is uniformly bounded in L2(0, T; V), which completes the proof. h In order to justify the convergence of the numerical trajectories, we need to consider the coefficient of frictions, l() as measures. The next Lemma 5 is necessary to arrive at the final step for the convergence theory. Lemma 5. There exists a constant C > 0 such that for sufficiently small ht > 0

Z

T

0

Z

l

lht ðrN Þht ðt; xÞdxdt 6 C:

0



RT Rl P Pn1 k

k

k k

Proof. Recall the approximate formula (5.6). Since 0 0 lht ðrN Þht ðt; xÞxdxdt ¼ lht P1 k¼0 l rN , we claim that ht k¼0 l rN 6 C. Putting the vector valued function w(x) :¼ x for any x 2 (0, l) we obtain





  lk rkN ðu_ kþ1 ðlÞ þ u_ k ðlÞÞl 1 kþ1 : ðu_ J 01=ht þ u_ k Þ ; w ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðu_ kþ1 ðlÞ þ u_ k ðlÞÞ2 þ 4ht

Modifying the numerical formulation (5.4) properly, we have

lht

P1 X k¼0

P1 P1 P1 P1 X X lk jrkN jðu_ kþ1 ðlÞ þ u_ k ðlÞÞ X ht X k qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi¼ ðv kþ1  v k ; wÞH þ hLe ðukþ1 þ uk Þ; wi þ ht hLv ðukþ1  uk Þ; wi  ht ðf ; wÞH 2 2 2 k¼0 k¼0 k¼0 ðu_ kþ1 ðlÞ þ u_ k ðlÞÞ þ 4ht k¼0

6 Cðkv P kH þ kv 0 kH Þ þ CT max kv 0 kH þ CðkuP kV þ ku0 kV Þ þ kfv kL2 ð0;T;L2 ð0;lÞÞ þ CT: 16k6P

ð5:11Þ

Note that (5.6) can be expressed in another way:

Z

T

0

Z

l

lht ðrN Þht ðt; xÞdxdt ¼ ht

0

n1 X

lk jrkN j:

k¼0

Since for sufficiently small ht > 0

lk jrkN jðu_ kþ1 ðlÞ þ u_ k ðlÞÞ ; lk jrkN j ’ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

ðu_ kþ1 ðlÞ þ u_ k ðlÞÞ2 þ 4ht

it follows from (5.11) that

Z 0

T

Z 0

l

lht ðrN Þht ðt; xÞdxdt 6 C;

as desired. h From the frictionless viscoelastic rod problem (2.1)–(2.5), the contact forces rN have been considered as measures and we can see that there is a subsequence, say ðrN Þht as measures such that ðrN Þht * rN as ht ; 0. Since nonnegative Borel measures mht can be identified with ðrN Þht , we set up them as

mht ðBÞ ¼

Z

B

ðrN Þht dxdt;

ð5:12Þ

where B is any Borel set in [0, T]  [0, l]. Thus it follows that ðrN Þht , we define the measure

.ht ðBÞ ¼

Z

B

mht * m as ht ; 0. Letting lht be a measure corresponding to

lht dmht :

By the previous Lemma 5, as ht ; 0,

.ht ð½0; T  ½0; lÞ ¼

Z 0

T

Z 0

l

lht dmht ¼

Z 0

T

Z 0

l

lht ðrN Þht ðt; xÞdxdt 6 C:

ð5:13Þ

7087

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Thus .ht * . as measures which implies that lht * l, since lht can be identified with .ht . Hence the coefficients of friction are measures in terms of the contact forces. Since . is absolutely continuous with respect to m, by Radon–Nikodym theorem (see [20, VII, 2] or [25, Section 11.6]) there is a unique l almost everywhere. This also makes a sense in the physical point of view, that is, if there is no contact, friction will never happen. The following weak ⁄ convergences can be shown, by applying Alaoglu theorem. Since hht and uht are uniformly bounded     in L1 0; T; H1L ð0; lÞ wL1 0; T; H1L ð0; lÞ , there are subsequences, say hht and uht such that hht * h and uht * u in     L1 0; T; H1L ð0; lÞ as ht ; 0. Similarly, since h_ ht and u_ ht are uniformly bounded in L2 0; T; H1L ð0; lÞ from Lemma 4, there are sub  sequences, say h_ ht and u_ ht such that h_ ht * h_ and u_ ht * u_ in L2 0; T; H1L ð0; lÞ as ht ; 0. Note that one can show that u_ 2 L2 ð0; T; VÞ, using other approaches, the Galerkin method or the penalty method. € ht ¼ ðu_ ht ðt þ h; Þ  u_ ht ðt; ÞÞ=ht for t 2 (tk1, tk]. Then using the similar The approximation of acceleration is defined by u € ht are uniformly bounded in L2(0, T; H1(0, l)), independently idea that we did in Lemma 5 we can see easily that € hht and u of any sufficiently small ht P 0. We introduce the following useful theorem which can be founded in Lion [21]. Theorem 6. Let p P 1 and q > 1. If W # U # Y with compact inclusion map W ? U and continuous inclusion map U ? Y and

_ Lq ð0;T;YÞ 6 Mg; S ¼ fu 2 Lp ð0; T; WÞju_ 2 Lp ð0; T; YÞ; kukLp ð0;T;WÞ þ kuk then S is precompact in Lp(0, T; U), i.e., each bounded sequence in S has a convergent subsequence in Lp(0, T; U). Applying Theorem 6, there is a subsequence, say u_ ht such that u_ ht ! u_ in L2(0, T; L2(0, l)) as ht ; 0. Therefore since _ lÞ in L2(0, T), we can assume that there is a subsequence u_ ht such that u_ ht ð; lÞ ! uð; _ lÞ pointwise a.e. on [0, T] u_ ht ð; lÞ ! uð; _ lÞj pointwise a.e. as ht ; 0. Then we can show easily that the solutions as ht ; 0. Thus we observe that wht ðu_ ht ð; lÞÞ ! juð; _ lÞ which are converged by such a subsequence u_ ht ð; lÞ as ht ; 0 satisfy the Coulomb friction law (2.15) and (2.16). uð; In order to show the existence of solutions, we need to consider their compactness. It is not hard to prove that uht belongs   to the Hölder space C 1=2 0; T; H1L ð0; lÞ using the estimate (5.9) and the following identity

uht ðtÞ ¼ uht ðsÞ þ

1 2

Z

t

v h ðs  ht Þ þ v h ðsÞds t

t

for s; t 2 ½0; T

s

which is easily derived from the numerical equation (5.5). The Sobolev embedding theorem allows us to see that H1L ð0; lÞ is compactly embedded in C[0, l]. Thus by the Arzela–Ascoli theorem, there is a subsequence, say uht such that uht ! u in C(0, T; C[0, T]) = C[0, T]  C[0, l]. On the other hand, if we impose strong conditions on the approximations uht , we can obtain the better Hölder spaces and improve their regularities. To see it, we define the characteristic function vQ which enables to agree uht on any set Q and vanish outside of Q.   Lemma 7. Let the time step size be ht > 0. Suppose that the time discrete approximate solutions uht 2 Hb 0; T; H1L ð0; lÞ for b P p + 1 with 0 < p < 1. Then uht are uniformly Hölder continuous from ½0; T ! H1L ð0; lÞ with the exponent p. Proof. We choose s, t such that 0 6 s < t 6 T and jt  sj < ht. Then we want to show that kuht ðtÞ  uht ðsÞkH1 ð0;lÞ 6 Cjt  sjp . Since L jeist  eissj 6 Cjt  sjpjsjp with 0 < p < 1, it follows from the definition of Sobolev space with fractional index (3.2) that

Z  kuht ðt; xÞ  uht ðs; xÞk2H1 ð0;lÞ ¼ C   L

1

1

2  is t iss  ð s ; xÞðe  e Þd s uc ht 

6 Cjt  sj2p 6 Cjt  sj2p 6 Cjt  sj2p 6 Cjt  sj2p

Z Z Z Z

H1L ð0;lÞ

1

1

ds 6 ð1 þ jsj2 Þbp

Z

1 1

1 1 þ jsj2

1

2  p b b  ð s ; xÞj s j h s i h s i d s uc ht 

H1L ð0;lÞ

p b b k uc ht ðs; xÞkH1 ð0;lÞ jsj hsi hsi ds L Z 1 2 2b k uc jsj2p hsi2b ds ht ðs; xÞkH1 ð0;lÞ hsi ds

1 1 1 1 1 1

L

2 2b k uc ht ðs; xÞkH1 ð0;lÞ hsi ds L

2 2b k uc ht ðs; xÞkH1 ð0;lÞ hsi ds

Note that for (b  p) P 1 1

1

2

1

1

Z

Z  6 Cjt  sj2p  

ds ¼ p:

L

1

Z

1

1

Z

1

1

ð1 þ jsj2 Þp ð1 þ jsj2 Þb 1

ds

ð1 þ jsj2 Þbp

ds

7088

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Thus, restricting the whole domain into [0, T], for b P p + 1 we obtain

kuht ðt; xÞ  uht ðs; xÞkH1 ð0;lÞ 6 Cjt  sjp

Z

1

p

L

1

2 2b k uc ht ðs; xÞ  v½0;T kH1 ð0;lÞ hsi ds

1=2

¼ Cjt  sjp kuht kHb ð0;T;H1 ð0;lÞÞ ¼ Cjt  sjp : L

L

This completes the proof. h In the following Lemma 8, the big O notation is used; we write f = O(g), provided that there is a constant C > 0 such that jf(x)j 6 Cjg(x)j for a.e. x 2 (0, l). Lemma 8. Assume that the first weak derivative of uht satisfies ðuht Þx ð; xÞ ¼ Oððx  lÞð2q#1Þ=2# Þ over the spacial domain (0,l) for 0 < q# 6 1/2 with 0 < # < 1. Then the time discrete approximate solutions uht are uniformly Hölder continuous from ½0; T ! H#L ð0; lÞ with exponent p = 1  #. Proof. We choose s, t 2 [0, T], as we did in Lemma 7. Then we want to show that kuht ðtÞ  uht ðsÞkH# ð0;lÞ 6 Cjt  sjp with L p = 1  #. It follows that

kuht ðtÞ  uht ðsÞk2H# ð0;lÞ ¼ C L

¼C

Z

1

1 Z 1 1

¼C

Z

1

1

2 2# c j uc ht ðt; nÞ  u ht ðs; nÞj hni vð0;lÞ dn 22# 2# 2# c c j uc j uc ht ðt; nÞ  u ht ðs; nÞj ht ðt; nÞ  u ht ðs; nÞj hni dn

Z





1 1

22#

2# 2# c ðuht ðt; nÞ  uht ðs; nÞÞeinx dx

j uc ht ðt; nÞ  u ht ðs; nÞj hni dn:

Restricting the whole domain into (0, l) and using the Hölder’s inequality, we can obtain

kuht ðtÞ  uht ðsÞk2H# ð0;lÞ 6 Ckuht ðt; xÞ  uht ðs; xÞkL22# 2 ð0;lÞ

Z

L

Z  ¼ C 

s

t

22# Z 

v h ðs; xÞds t

2

L ð0;lÞ

1

1

1

1

2# 2# c j uc ht ðt; nÞ  u ht ðs; nÞj hni dn

2# 2# c j uc ht ðt; nÞ  u ht ðs; nÞj hni dn:

Then by the energy bounds, we have

kuht ðtÞ  uht ðsÞk2Hh ð0;lÞ 6 C # jt  sj22# L

Z

1 1

2# 2# c j uc ht ðt; nÞ  u ht ðs; nÞj hni dn:

Note that C# is a constant independent of s, t, and ht. Since 1/(jnj2 + 1) 6 1, we have

kuht ðtÞ  uht ðsÞk2H# ð0;lÞ ¼ C # jt  sj22# L

6 C # jt  sj22#

Z

1

1

Z

1

1

6 C # jt  sj22#

Z

1

1

2# 2 2#2 c j uc dn ht ðt; nÞ  u ht ðs; nÞj hni hni

2# 2 c j uc ht ðt; nÞ  u ht ðs; nÞj hni

!1#

1 jnj2 þ 1

dn

2# 2 22# c j uc kðuht ðtÞ  uht ðsÞÞ# k2H1 ð0;lÞ : ht ðt; nÞ  u ht ðs; nÞj v½0;l hni dn ¼ C # jt  sj L

Since uht is the linear interpolant in the time space [0, T] and jðuht Þx ð; xÞj 6 Cjx  ljð2q#1Þ=2# with q > 0 over (0, l), using the inequality kukH1 ð0;lÞ 6 Ckux kL2 ð0;lÞ we can obtain L

kuht ðtÞ  uht ðsÞk2H# ð0;lÞ 6 C # jt  sj22# L

Z 0

l

1 ðx  lÞ

dx ¼ C # jt  sj22# 2q#þ1

Z 0

l

1 dx: x2q#þ1

In order to bound the improper integral, we need to impose the condition 0 6 2q# + 1 < 1. Therefore if ðuht Þx ð; xÞ ¼ Oððx  lÞð2q#1Þ=2# Þ with 0 < q# 6 1/2, it follows that kuht ðtÞ  uht ðsÞkH# ð0;lÞ 6 Cjt  sjp , where p = 1  #. This completes the L proof. h

6. Energy balance In this section, we consider a possible energy balance with Coulomb friction. Friction process is inherently dissipative. More specifically, we observe in our problem that the energy dissipation is due to not only viscosity but also friction. As we discussed in the previous Section 5, the coefficient of friction was considered as a measure. However, in order to

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7089

investigate an energy balance, we impose a stronger condition that l is Lipschitz continuous in terms of the sliding velocity at the right end x = l. Thus, the Coulomb law of dry friction is formulated as follows:

_ lÞjÞjrN ðtÞj on ð0; T; jrT ðtÞj 6 lðjuðt; 

ð6:1Þ

_ lÞjÞjrN ðtÞj ) uðt; _ lÞ ¼ krT ðtÞ for some k P 0; jrT ðtÞj ¼ lðjuðt; _ lÞjÞjrN ðtÞj ) uðt; _ lÞ ¼ 0; jrT ðtÞj < lðjuðt;

ð6:2Þ

_ lÞ jÞ P 0 is called the coefficient of friction dependent on the slip rate juðt; _ lÞj. Assuming that there are where lðtÞ :¼ lðj uðt; Ml, Ll > 0 such that klkL1 ðRþ Þ 6 M l and klðr2 Þ  lðr 1 ÞkL1 ðRþ Þ 6 Ll jr2  r1 j for r 1 ; r2 2 Rþ , Kuttler and Shillor [19] showed the existence of solutions; they considered l as the set-valued function, and used the Lipschitz continuous approximation of the set-valued function, and extended the slip rate dependent friction into the d-dimensional case that a rigid obstacle moves. We consider the restriction or trace operator c such that c : H1L ð0; lÞ ! R with cðbwÞ ¼ wðlÞ 2 R. Then we set up the following identity

Z

T

hc n; bwidt ¼

0

Z

T

hn; cðbwÞidt;

ð6:3Þ

0

where c : R ! H1 ð0; lÞ with c⁄(n) = n(t)d(x  l) is an adjoint operator. Although the maximum dissipation principle can be defined in d-dimensional problems (its main idea is based on Theorem 2.1 in the paper [19]), we want to consider it in only one-dimensional case. Let n ¼ rT 2 Mð½0; TÞ. Instead of dealing with the averaging operator, we can define the maximum dissipation principle by

_ ¼ lðjuðt; _ lÞjÞjrN ðtÞkuðt; _ lÞj for almost all t 2 ½0; T: hn; cðbuÞi Recalling (3.5), it is easy to see from monotone convergence theorem that for any z – 0 2 V

Z

T

hc n; bwidt ¼

0

Z 0

T

  lim J 0n ðbzÞ; bu_ dt ¼ lim

n!1

Z

n!1

0

T

 0  J n ðbzÞ; bu_ dt:

ð6:4Þ

In the viscoelastic frictional (Coulomb law) contact problems, the rate of energy loss is attributed to viscous damping, external forces, and frictional heat generation. The energy balance is considered, based on the following variational inequality that is equivalent to the variational form (4.4): find u(t) 2 V with almost all t 2 [0, T]

€ ðtÞ þ Le uðtÞ þ Lv uðtÞ _ hu þ c nðtÞ  fðtÞ; wðtÞ  uðtÞiV  V P 0 for all wðtÞ 2 V: _ 2 K, where K is a closed and convex set in a separable reflexive Banach We also need an additional assumption that cðbuÞ _ in the strong sense, it is straightforward to get space. Provided that there exists cðbuÞ

€ ðtÞ þ Le uðtÞ þ Lv uðtÞ _ _ hu ¼ 0; þ c nðtÞ  fðtÞ; uðtÞi which plays a significant role in showing the energy balance. However, in dynamic contact problems the solutions are not differentiable in the strong sense. Therefore more sophisticated approach needs to be considered in the next Lemma 9. _ 2 K #V Lemma 9. Assume that c⁄n(t) 2 L2(0, T; V⁄) and the solutions u : [0, V] ? V satisfy the variational inequality: for cðbuÞ

Z 0

T

€ ðtÞ þ Le uðtÞ þ Lv uðtÞ _ hu þ c nðtÞ  fðtÞ; wðtÞ  uðtÞidt P 0 for all cðbwÞ 2 K:

Provided that there exists u_ 2 L2 ð0; T; VÞ in the distributional sense, then for almost all t 2 [0, T] we have the following form of energy balance;

1 1 _ _ _ _ ðhuðtÞ; uðtÞi þ hLuðtÞ; uðtÞiÞ  ðhuð0Þ; uð0Þi þ hLuð0Þ; uð0ÞiÞ 2 2 Z t _ sÞi  hLv uð _ sÞi  hc nðsÞ; uð _ sÞids: _ sÞ; uð hfðsÞ; uð ¼

ð6:5Þ

0

Proof. By applying Lemma 2.1 in the paper [27], we can obtain

€ ðtÞ þ Le uðtÞ þ Lv uðtÞ _ _ hu ¼ 0 for almost all t 2 ½0; T: þ c nðtÞ  fðtÞ; uðtÞi Thus we can rewrite it as follows;

1 d _ _ _ _ ðhuðtÞ; uðtÞi þ hLe uðtÞ; uðtÞiÞ ¼ hfðtÞ  Lv uðtÞ  c nðtÞ; uðtÞi: 2 dt

ð6:6Þ

Taking an integral on both sides of (6.6) over [0, t] with 0 < t 6 T, we use the fundamental theorem to arrive at Eq. (6.5). h

7090

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

The energy balance discussed above will be supported by numerical results in Section 8. The author conjectures that it may be possible to extend this energy balance for Coulomb friction into d-dimensions. 7. Fully discrete numerical scheme In order to compute fully approximations, the Galerkin method will be used in the spacial domain [0, l] and the time discretization on the time space [0, T], as we discussed in Section 5. When we employ the Galerkin method with piecewise linear functions on the interval [0, l], our finite dimensional space V hx  H1L ð0; lÞ becomes

V hx ¼ fwhx 2 C½0; ljwhx ð0Þ ¼ 0g; where the nodes xi are chosen with the uniform size of the subintervals, i.e., hx = l/m and xi = ihx for 1 6 i 6 m. Most basis functions wi(x) associated with nodes x1, x2, . . . , xm1 are the piecewise linear functions. Then the discrete approximate solu_ u; u_ in the space [0, l] are denoted by hh ; xh ; uh ; v h respectively. In order for them to satisfy tions for the actual solutions h; h; x x x x the natural boundary conditions, we construct the last basis function wm(x) for uhx ; v hx :

( wm ðxÞ ¼

1 hx

ðx  xm1 Þ if xm1 6 x 6 xm ;

0

otherwise;

and the last basis function /m(x) for hhx ; xhx :

( /m ðxÞ ¼

 h13 ðx  xm1 Þ3 þ h12 ðx  xm1 Þ2 þ h1x ðx  xm1 Þ if xm1 6 x 6 xm ; x

x

0

otherwise:

Thus setting wi = /i for 1 6 i 6 m  1, we consider two finite dimensional spaces W hx ¼ spanf/i j1 6 i 6 mg  H1L ð0; lÞ and V hx ¼ spanfwi j1 6 i 6 mg  H1L ð0; lÞ. Now for each time step tk we write the fully approximate solutions as

ukhx ðxÞ ¼

m X

ukj wj ðxÞ;

v kh ðxÞ ¼ x

j¼1

hkhx ðxÞ ¼

m X

m X

v kj wj ðxÞ;

j¼1

hkj /j ðxÞ;

xkhx ðxÞ ¼

j¼1

m X

xkj /j ðxÞ;

j¼1

where the coefficients of the basis functions are ukj ¼ uj ðt k Þ; v kj ¼ v j ðtk Þ; hkj ¼ hj ðt k Þ, and xkj ¼ xj ðtk Þ. Recalling the semidiscrete numerical formulations (5.4) and (5.5) and assuming that f = 0, we establish the fully discrete numerical formulations:

" kþ1 #  k 1 xhx  xhx 1 A1 ; ¼ ht v kþ1 2 A2  v kh h

ð7:1Þ

" kþ1 # " kþ1 # k k 1 h hx  h hx 1 xhx þ xhx ¼ ; k ht uhkþ1  ukh 2 v kþ1 hx þ v hx x x

ð7:2Þ

x

x

where

              k A1 ¼ hkþ1 þ hkhx þ ukþ1 þ ukhx x  hkþ1 xkþ1 þ xkhx xx þ v kþ1 þ v khx x  hhkþ1  hkhx and hx hx  hhx þ a2 hx hx hx x xx xx x xx x 

k  kþ1 k k           l rN uhx ðlÞ  uhx ðlÞ     ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : þ ukhx xx  hhkþ1  hkhx þ a2 v hkþ1 þ v khx xx  xkþ1  xkhx x  r A2 ¼ uhkþ1 hx x x x 2 xx x x xx x kþ1 k 2 uhx ðlÞ  uhx ðlÞ þ e Note that in actual computation the smoothing parameter e > 0 which is small enough is used instead of ht > 0. By using Eq. (7.2), from (7.1) we obtain the following equations at one time step:

2

þ 2

ht

!

 1 a2 kþ1 1 a2  kþ1  hhx  hhx þ þ ¼ xx 2 ht 2 ht

!



 1 a2  k  1 a2  k  2 1 a2  kþ1  hhx þ h hx uhx þ  þ þ xkhx þ xx x 2 ht ht 2 ht ht 2 ht

  1 a2 ukhx x  þ 2 ht 2

 2

ð7:3Þ

7091

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

and







 1 a2  kþ1  2 k 1 a2  k  1 a2  kþ1  1 a2  k  2 kþ1 u u h hhx þ v khx þ  þ  u  ¼ u þ   h h h h h 2 2 x x t x x xx xx x x 2 2 2 2 h h h h h t t t t t ht ht 

 kþ1 k k

k l rN uhx ðlÞ  uhx ðlÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi :  r 2 k 2 ukþ1 ðlÞ  u ðlÞ þ e hx hx 2

ð7:4Þ

Multiplying both sides of (7.3) by the basis functions /i(x) and both sides of (7.4) by wi(x) for 1 6 i 6 m respectively, two next ~ kþ1 ÞT can be computed via the Galerkin approximation: step solutions ðhkþ1 ; u

! " !

 #



 # 1 a2 1 a2 2 1 a2 a2 1 2 1 a2 kþ1 ~ kþ1 M1 þ M1 þ K1 h L1 u þ K1 hk þ M1 xl þ þ þ þ ¼  þ  2 2 2 ht ht 2 ht ht 2 ht 2 ht ht 2 ht

 1 a2 ~k L1 u þ  2 ht 2

30

30

25

25

15 10

15 10 5

5 0

ð7:5Þ

20

20

Energy

Energy

"

0 0

0.5

1

1.5

2

−5

0

0.5

Time

1

Time

Fig. 8.1. Energy functions.

Fig. 8.2. Motion of the right end (contact point) of the beam a = 0.001.

1.5

2

7092

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

Fig. 8.3. Motion of the right end (contact point) of the beam a = 2.

Fig. 8.4. Motion of the beam from t = 0 to t = 0.6 with a = 0.001.

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7093

and

"

2

M2 þ 2

ht

"

 #

 #



 1 a2 2 a2 1 2 1 a2 T kþ1 1 a2 T k k k ~ kþ1 ¼ ~ ~ K K2 u L L h M þ  þ M v þ h þ þ þ  u 2 2 2 2 2 ht ht 2 ht 2 2 ht 2 ht 2 ht 0 1T

 k lk rkN ukþ1 B m  um C ffiA ;  @0;    ; 0; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  kþ1  uk 2 þ e2 um m

ð7:6Þ

where the mass matrices M1, M2 and the stiffness matrices K1, K2 and the matrices L1 are presented as follows:

Z

M1 :¼ ðM 1 Þji ¼ K1 :¼ ðK 1 Þji ¼ L1 :¼ ðL1 Þji ¼

Z

Z

0 l

0

l 0 l

/j ðxÞ/i ðxÞdx and M2 :¼ ðM2 Þji ¼ /0j ðxÞ/0i ðxÞdx and K2 :¼ ðK 2 Þji ¼

Z

Z

l

wj ðxÞwi ðxÞdx;

0 l

0

w0j ðxÞw0i ðxÞdx;

w0j ðxÞ/i ðxÞdx

 T  T k  k ~ ¼ u1 ; . . . ; ukm1 ; ukm T ; Fully discrete approximate solutions are denoted by hk ¼ hk1 ; . . . ; hkm1 ; hkm ; xk ¼ xk1 ; . . . ; xkm1 ; xkm ; u  T v~ k ¼ v k1 ; . . . ; v km1 ; v km . Letting

!



 1 a2 2 1 a2 a2 1 M K1 ; B ¼ þ þ K2 ; þ þ þ 2 2 2 ht ht 2 ht ht 2 ht !



 2 1 a2 a2 1 2 a2 1 M1 þ C¼ K1 ; and D ¼ 2 M2 þ K2 ; þ    2 ht 2 ht 2 ht 2 ht ht

A¼

2

M1 þ 2

Fig. 8.5. Motion of the beam between t = 0.6 and t = 1.5 with a = 0.001.

7094

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

kþ1 we can set up the following nonlinear equation in terms of um :

"

hkþ1 P kþ1 ~ u

#

"

hk ¼Q k ~ u

#

 2 M1 þ ht 0

0 M2

"

xk

v~ k

#

2

3 0 kþ1 k 6 l jr jðum um Þ 7  4q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5; 2 ðumkþ1 ukm Þ þe2 k

k N

ð7:7Þ

where

2

A

6  P¼4   12 þ ah2t LT1

  3  12 þ ah2t L1 7 5; B

2 6 Q ¼ 4

1 2

C   ah2t LT1



1 2

 3  ah2t L1 7 5: D

Here the banded matrices P; Q ; S 2 R2m2m and 0 2 Rð2m1Þ1 in the last term of Eq. (7.7). Note that P is a positive definite matrix.  T ~ kþ1 from (7.7). We introduce the block matrix notation to deal with the nonLet the next step solution be ykþ1 ¼ hkþ1 ; u linear equation (7.7) efficiently. For instance, if we extract i1th row to i2th row and j1th column to j2th column for ði i þ1Þðj2 j1 þ1Þ 1 6 i1 6 i2 6 2m and 1 6 j1 6 j2 6 2m, then will be the block matrix. Let the block matrix T  Pði1 : i2 ; j1 : j2 Þ 2 R 2 1 T k k k k k  be T = P(1 : 2m  1, 1 : 2m  1) and y ¼ h1 ; . . . ; hm1 ; hm ; u1 ; . . . ; ukm1 . Then through several algebraic manipulations from (7.7) and we can set up the nonlinear equation in terms of ukþ1 m : 1

0 ¼ Pð2m; 1 : 2m  1ÞT



k

r 

ukþ1 m Pð1

 k lk rkN ukþ1 m  um kþ1 k ffi; : 2m  1; 2mÞ þ um Pð2n; 2nÞ  r m  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  k 2 þ e2 ukþ1 m  um

Fig. 8.6. Motion of the beam between t = 1.5 and t = 2 with a = 0.001.

ð7:8Þ

7095

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

where rkn ¼ Q ð2m; 1 : 2m  1Þyk þ Q ð2m; 2mÞukm þ Sð2m; m þ 1 : 2m  1Þv k þ v kn Sð2m; 2mÞ and rk ¼ Q ð1 : 2m  1; 1 : 2m  1Þ k þ v kn Sð1 : 2m  1; 2mÞ. To compute a better approximation of Eq. (7.8), yk þ Q ð1 : 2m  1; 2mÞ þ Sð1 : 2m  1; 1 : 2m  1Þy we define the nonlinear function F to be

 kþ1 F um ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    k 2 þ e2 ukþ1 Pð2n; 2nÞ  r k þ Pð2m; 1 : 2m  1ÞT1 rk  ukþ1 Pð1 : 2m  1; 2mÞ ukþ1 m  um m m m

 kþ1  ukm :  lk rkN um

ð7:9Þ



kþ1 kþ1 When we solve the next step solution of the end of the beam, um for the nonlinear equation F um ¼ 0, Newton’s method kþ1 kþ1  with Gaussian elimination method is used. When um is computed, y is obtained from the linear system:

 kþ1 kþ1 ¼ T1 rk  um y Pð1 : 2m  1; 2mÞ :

ð7:10Þ

~ kþ1 ÞT are obtained and the next step velocity vk+1 and next step angular speed xk+1 can be Thus the next step solutions ðhkþ1 ; u computed from (7.2):

v~ kþ1 ¼

2 kþ1 ~ ~kÞ  v ~ k; ðu u ht

xkþ1 ¼

2 kþ1 ðh  hk Þ  xkþ1 : ht

ð7:11Þ

We summarize our numerical scheme proposed above. See the following Algorithm 10, assuming that the number of time steps, T/ht are positive integers. ~ 0; v ~ 0 Þ are given Algorithm 10. The initial data ðh0 ; x0 ; u for k = 0 : T/ht  kþ1 ¼ 0 which is set up in (7.9) Solve the nonlinear equation F um ~ kþ1 Þ from the linear system (7.10) Compute the next step solutions ðhkþ1 ; u ~ kþ1 Þ from (7.11) Compute the next step solutions ðxkþ1 ; v end for

Fig. 8.7. Motion of the beam from t = 0 to t = 0.6 with a = 2.

7096

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

In order to consider how stable the fully discrete approximations are, we present the following Lemma 11 that will be supported by the numerical results concerning the energy balances in Fig. 8.1. The approximations which are computed with our fully numerical scheme have to guarantee dissipation of energy. ~k; v ~ k  in the fully discrete case: We also define the energy functional Eðt k Þ :¼ E½hk ; xk ; u

~k; v ~ k ¼ E½hk ; xk ; u

1 1 ~ k ÞT K2 u ~ k ÞT L1 hk þ ðhk ÞT M1 hk Þ: ~ k  2ðu ~ k Þ þ ððhk ÞT K1 hk þ ðu ~ k ÞT M2 v ððxk ÞT M1 xk þ ðv 2 2

ð7:12Þ

~k; v ~ k Þ with the viscosity quantity a2 P 0 satisfy Eqs. (7.1) Lemma 11. Suppose that fv = 0 and the fully approximations ðhk ; xk ; u ~ kþ1 ; v ~k; v ~ kþ1  6 E½hk ; xk ; u ~ k  for all integers and (7.2), where a smoothing parameter e > 0 is sufficiently small. Then E½hkþ1 ; xkþ1 ; u k P 0. Proof. Partitioning  the time space  [0, T], we first set up the numerical formulation of (7.1) and then use the fully discrete approximations ukhx ; v khx ; hkhx ; xkhx at each time step t = tk. Applying Eq. (7.2) and integration by parts, we can obtain

1 ~ kþ1  ðv ~kÞ ~ kþ1 ÞT M2 v ~ k Þ T M2 v ððxkþ1 ÞT M1 xkþ1  ðxk ÞT M1 xk þ ðv 2ht 1 1 ~ kþ1 ÞT K2 u ~ k ÞT K2 u ~ kþ1 ÞT L1 hkþ1  ðu ~ k ÞT L1 hk Þ ~ kþ1  ðu ~ k Þ  ððu ¼ ððhkþ1 ÞT K1 hkþ1  ðhk ÞT K1 hk þ ðu 2ht ht 1 kþ1 a a ðh þ hk ÞT M1 ðhkþ1  hk Þ  ðv kþ1 þ v k ÞT K2 ðv kþ1 þ v k Þ  ðxkþ1 þ xk ÞT K1 ðxkþ1 þ xk Þ  2ht 4 4 h kþ1  kþ1 i 1  kþ1 kþ1 k k k k kþ1 u  um um þ um  hm  hm þ a v m þ v m  xm  xkm : þ 2ht m

ð7:13Þ

Then the last factor in the last term of the previous Eq. (7.13) is replaced by the numerical tangential force by the natural boundary condition (2.14). For the sufficiently small e > 0, the numerical tangential force satisfying the Coulomb friction law is approximated by



 k lk rkN ukþ1 1 _ kþ1 _ k m  um ~ ~ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi: q ðu þ u Þ ¼  Je 2 k 2 þ e2 ukþ1  u m m 0

Fig. 8.8. Motion of the beam between t = 0.6 to t = 1.5 with a = 2.

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7097

Therefore we have

~ kþ1 ; v ~ kþ1 ÞT L1 hkþ1 þ ðhkþ1 ÞT M1 hkþ1 ~ kþ1 þ ðhkþ1 ÞT K1 hkþ1  2ðu ~ kþ1  ¼ ðxkþ1 ÞT M1 xkþ1 þ ðv ~ kþ1 ÞT M2 v 2E½hkþ1 ; xkþ1 ; u

 lk rkN umkþ1  ukm 2 ~ k ÞT L1 hk þ ðhk ÞT M1 hk  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ k þ ðhk ÞT K1 hk  2ðu ~ k Þ T M2 v ffi 6 ðxk ÞT M1 xk þ ðv  kþ1  uk 2 þ e2 um m ~k; v ~ k ; 6 2E½hk ; xk ; u as desired. h

8. Numerical results and discussion In this Section, we implement our numerical schemes that we discussed in the Section 7. Our numerical experiments are performed using the time interval [0, 2] and the spacial domain [0, 5] and the regularized parameter e = 102/5. Thus, the length of the Timoshenko beam becomes l = 5 and the end time T = 2. Since the beam is clamped at its left end and moves downward initially, we will use the following data: the initial displacement of the viscoelastic Timoshenko beam is u0(x) = x/200 and its initial velocity is v0(x) = x. Notice that we use dimensionless constants. When the numerical scheme is implemented, the size of subinterval hx = 0.005 and the time step size ht = 0.002 are used and all coefficients are used as follows: q I = qA = 1 and EI = KAG = 100. In addition, we let the coefficients of viscosity, a2 = 0.001 so the beam is nearly elastic and a2 = 2 for the viscoelastic beam with large damping. The contact force rN(t) may computed from (5.1)–(5.5). However, our main interest lies in the motion of the Timoshenko beam with the frictional contact at its right end. So the initial contact forces rN(t) are given by rN(0) = 100,000 and rN(t) = 100,000/100t. Those are chosen so that the time period that the end of the beam touches the _ lÞj=100 which satisfies the slip rate stationary obstacle is short. The coefficient of friction used in our simulation is lðtÞ ¼ juðt; dependent friction. According to our several tests, oscillating of the beam at the right end depends on the magnitude of contact forces; larger contact forces causes more oscillations. When contact is lost, the velocity of slip vanishes. In this numerical computation, Matlab software is used to compute the numerical solutions. Especially, sparse matrix utility is used to handle the banded matrices efficiently which are constructed from the linear systems, because it compresses the zero elements.

Fig. 8.9. Motion of the beam between t = 1.5 and t = 2 with a = 2.

7098

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

 kþ1 When we solve the nonlinear equation F um ¼ 0 (see (7.9)), we use the stopping criterion  = 1010. In our numerical experiments, the number of the average iteration counted at each time step in Newton’s method is 2.0250 for the almost elastic beam (a = 0.001) and 2.0125 for the viscoelastic beam (a = 2). Therefore we conclude that Newton’s method employed here works very well. We turn to the numerical results. Two graphs are drawn in Fig. 8.1, which supports dissipation of energy proved in Lemma 9. The outstanding difference between the two energy functions is that the energy function with viscoelastic case drops quite drastically, due to the loss of viscosity. The dashed line (in red) displays conservation of energy without viscosity and friction effects. In our numerical results, energy loss for friction seems to be negligible as compared to the viscosity damping. Numerical simulations of the frictional beam problem are presented through Figs. 8.2–8.9. We begin by describing the motion of the right end (the contact point x = 5). In Fig. 8.2 the nearly elastic beam (a = 0.001) is shown and in Fig. 8.3 the viscoelastic beam (a = 2) is shown. The motion of Timoshenko beams for both visco quantities a = 0.001 and a = 2 is simulated in Figs. 8.4, 8.5, 8.6, 8.7, 8.8, 8.9. Each subfigure consists of two pictures; the top one is zoomed out so that the behavior near the end can be seen clearly, and the bottom one is the whole picture of the beam. As we can see, the angular speed and velocity of the almost elastic beam at the pictures (B) and (D) in Fig. 8.4, waves are displayed. This is so since the equation for the almost elastic beam is close to the hyperbolic type. Fig. 8.4 shows that the (almost) elastic beam moves downward and Figs. 8.5 and 8.6 show that it moves upward and downward. We can observe the similar motion of the viscoelastic beam through Figs. 8.7, 8.8, 8.9. Comparing two kinds of the beam, the elastic beam moves more quickly in the vertical direction with more vibrations, as one might expect. Even though some indications for numerical stability of the scheme are discussed above, the error analysis for this dynamic contact problem remains an unresolved issue. In addition to the issue, proving the uniqueness of the solution for dynamic contact problems (with Signorini conditions) is an open question. It is interesting to note that numerical results get worse, as the beam gets close to being purely elastic, i.e., as a ; 0. When we attempt to compute the numerical solutions with very small a, numerical results show an unreasonable increase of the energy function, which does not make any sense in terms of our proposed numerical schemes, since the solutions are supposed to demonstrate the decreasing energy. Consequently, we conjecture that proving the existence of solutions with frictional contact for purely elastic bodies would be a most difficult task. 9. Conclusion In this work, we established the existence of weak solutions to the problem of a frictional viscoelastic Timoshenko beam and investigated an energy balance with Coulomb law of dry friction. The important issue in deriving the energy balance is that imposing the slip rate dependence on the coefficient of friction is necessary, which has been justified theoretically and numerically. Another possible interesting question is on whether our energy balance for long duration of friction works or not. A new form energy balance would have to be established taking into account thermal effects, which are important in many applications. In order to take a new energy balance into consideration, we would need an additional equation of parabolic type to deal with frictional heat generation. This matter will be studied in the sequel. Acknowledgments The author would like to thank the comments of anonymous referees which have improved the presentation of this paper. References [1] Jeongho Ahn, A vibrating string with dynamic frictionless impact, Appl. Numer. Math. 57 (8) (2007) 861–884. [2] Jeongho Ahn, Thick obstacle problems with dynamic adhesive contact, M2AN Math. Model. Numer. Anal. 42 (6) (2008) 1021–1045. [3] Jeongho Ahn, David E. Stewart, Euler–Bernoulli beam with dynamic contact: discretization, convergence, and numerical results, SIAM J. Numer. Anal. 43 (4) (2005) 1455–1480 (electronic). [4] Jeongho Ahn, David E. Stewart, Existence of solutions for a class of impact problems without viscosity, SIAM J. Math. Anal. 38 (2006) 37–63 (electronic). [5] Jeongho Ahn, David E. Stewart, Frictionless dynamic contact in linear viscoelasticity, IMA J. Numer. Anal. 29 (1) (2009) 43–71. [6] Jeongho Ahn, David E. Stewart, A viscoelastic Timoshenko beam with dynamic frictionless impact, Discrete Cont. Dyn. Syst. Ser. B 12 (1) (2009) 1–22. [7] Kevin T. Andrews, M. Shillor, S. Wright, On the dynamic vibrations of an elastic beam in frictional contact with a rigid obstacle, J. Elast. 42 (1996) 1–30. [8] Kevin T. Andrews, M. Shillor, S. Wright, On the dynamic vibrations of an elastic beam in frictional contact with a rigid obstacle, J. Elast. 42 (1) (1996) 1– 30. [9] Marius Cocou, Existence of solutions of a dynamic Signorini’s problem with nonlocal friction in viscoelasticity, Z. Angew. Math. Phys. 53 (2002) 1099– 1109. [10] Y. Dumont, K.L. Kuttler, M. Shillor, Analysis and simulations of vibrations of a beam with slider, J. Eng. Math. 47 (2003) 61–82. [11] M. Erdmann, On a representation of friction in configuration space, Int. J. Robotics Res. 13 (3) (1994) 240–271. [12] I.R. Ionescu, J.-C. Paumier, On the contact problem with slip rate dependent friction in elastodynamics, Eur. J. Mech. A/Solids 13 (4) (1994) 555–568. [13] J. Jarušek, C. Eck, Dynamic contact problems with small Coulomb friction for viscoelastic bodies. existence of solutions, Math. Methods Appl. Sci. 9 (1) (1990) 11–34 (electronic). [14] J.U. Kim, A boundary thin obstacle problem for a wave equation, Commun. Partial Differ. Equat. 14 (1989) 1011–1026. [15] J.U. Kim, A one-dimensional dynamic contact problem in linear viscoelasticity, Math. Methods Appl. Sci. 13 (1990) 55–79. [16] K.L. Kuttler, Y. Renard, M. Shillor, Models and simulations of dynamic frictional contact of a beam, Comput. Methods Appl. Mech. Engrg. 177 (1999) 259–272 (electronic).

J. Ahn / Applied Mathematics and Computation 218 (2012) 7078–7099

7099

[17] Kenneth Kuttler, Meir Shillor, A dynamic contact problems in one-dimensional thermoviscoelasticity, Nonlinear World 2 (3) (1995) 355–385. [18] Kenneth Kuttler, Meir Shillor, Dynamic bilateral contact with discontinuous friction coefficient, Nonliear Anal. Theory Math. Appl. 45 (2001) 309–327. [19] Kenneth Kuttler, Meir Shillor, Dynamic contact with Signorini’s condition and slip rate dependent friction, Electron. J. Differ. Equat. 2004 (83) (2004) 1– 21. [20] Serge Lang, Real and Functional Analysis, Graduate Texts in Applied Mathematics, vol. 142, Springer-Verlag, New York, Berlin, Heidelberg, 1993. [21] J.L. Lions, Quelques Methods de Resolution des Problemes aux Limites Non Lineaires, Dunod, Paris, 1969. [22] J.-C. Paumier, T. Renard, Surface perturbation of an elastodynamic contact problem with friction, Eur. J. Appl. Math 14 (2003) 465–483. [23] Michael Renardy, Robert C. Rogers, An Introduction to Partial Differential Equations, Texts in Applied Mathematics, vol. 13, Springer-Verlag, New York, Berlin, Heidelberg, 1993. [24] M. Rous, P. Chabrand, F. Lebon, Numerical methods for frictional contact problems and applications, J. Theoret. Appl. Mech. 7 (suppl) (1988) 111–128. [25] H.L. Royden, Real Analysis, Prentice Hall, Inc, New Jersey, USA, 1987. [26] M. Schatzman, M. Bercovier, Numerical approximation of a wave equation with unilateral constraints, Math. Comput 53 (1989) 55–79. [27] David E. Stewart, Energy balance for viscoelastic bodies in frictionless contact, Quarterly Appl. Math. 67 (4) (2009) 735–743. [28] Michael E. Taylor, Partial Differential Equations 1, Applied Mathematical Sciences, vol. 115, Springer-Verlag, New York, 1996. [29] J. Wloka, Partial Differential Equations, Cambridge University Press, 1987.