Mathematical and numerical approaches to a one-dimensional dynamic thermoviscoelastic contact problem

Nonlinear Analysis: Real World Applications 22 (2015) 437–451

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa

Mathematical and numerical approaches to a one-dimensional dynamic thermoviscoelastic contact problem Jeongho Ahn ∗ , Natanya Clark Department of Mathematics and Statistics, Arkansas State University, P.O. Box 70, State University, AR 72467, United States

article

info

Article history: Received 28 April 2014 Accepted 29 September 2014 Available online 3 November 2014 Keywords: Dynamic contact Thermoviscoelasticity Normal compliance Barber’s heat exchange condition

abstract In this paper, we propose numerical schemes for solving a nonlinear system which consists of a coupled partial differential equations and two conditions, called normal compliance contact condition and Barber’s heat exchange condition. The convergence of numerical trajectories is shown by using a time discretization and passing the limit of the time step size. The uniqueness of the weak solution is proved as well. We derive the extensive form of an energy balance which will be a criterion to examine numerical stability. An example is provided to present and discuss numerical results. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction This one-dimensional dynamic contact problem is basically formulated by a coupled partial differential equations (PDEs) of linear thermoelasticity: an equation of motion and an equation of heat conduction, as was derived in Carlson’s paper [1]. In general, this complicated nonlinear system may be considered as a simple example of thermodynamics. In everyday life, its applications can be found easily: for example, heating and cooling houses, cooking, regulating and mainlining our body temperature and so on. The recent research on thermodynamics can provide a wide range of applications in applied sciences such as biology, biochemistry and automobile engineering. However, all laws of thermodynamics are not applicable to atom or molecular levels and thus it does not explain the rate of chemical reaction. The book [2, section 4.1] presents the general formalism of thermodynamical systems, and derives an equation which is thermodynamically consistent, by applying the second law of thermodynamics. The formalism can be applied to the PDEs to see their thermodynamic consistency. However, there are still some mathematical issues (see [2, section 4.5]), when we consider thermodynamics with side effects such as friction or wear. Concerning the thermodynamic consistency, readers may find more references (e.g., [3]). Indeed, William Day (see his book [4]) introduces more advanced mathematical treatment for the one-dimensional PDEs with linear thermoelasticity. Unlike the standard boundary conditions considered in Day’s book, we impose more applicable boundary conditions at two ends of the rod to complete our contact model. As for the usual contact problems of beams, the homogeneous essential boundary conditions are used at the left end. In addition to them, there are two natural boundary conditions at the right end which arise in our physical situation: one condition is normal compliance used in many papers (e.g., [5–8] and references therein) and another is Barber’s heat exchange condition (see the papers [9,10]). Similar types of quasistatic contact problems were previously studied in many papers (e.g., [11–15]), applying the Signorini’s contact

∗

Corresponding author. Tel.: +1 870 972 3090; fax: +1 870 972 3950. E-mail address: [email protected] (J. Ahn).

http://dx.doi.org/10.1016/j.nonrwa.2014.09.017 1468-1218/© 2014 Elsevier Ltd. All rights reserved.

438

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

Fig. 2.1. Dynamic contact of a thermoviscoelastic rod with a deformable obstacle.

conditions instead of normal compliance. Those papers except [12] use different but unrealistic heat exchange conditions which are simpler than Barber’s condition. Especially, in [15] a numerical scheme is proposed based on the Crank–Nicolson method and its convergence is proved as well. In this work, we attach the viscous term to reformulate the equation of motion, due to the fact that proving the existence of solutions for the elastic case does not seem to be successful. From an engineering point of view, the viscosity is neither realistic nor important for most of engineering materials. Therefore, we will perform numerical experiments with a small quantity to approximate the purely thermoelastic case. Another mathematical difficulty to show the existence of solutions is to handle the heat exchange coefficient function in the radiation condition. Barber’s three thermal contact conditions can be derived, by approximating a set-valued mapping. See the paper [12] for this issue. Even though the coefficient function is an unusual function, the existence of solutions for the quasistatic case can be proved in the paper [12], considering it as a measurable function. However, in the dynamic case, nicer assumptions on the coefficient function are required to show the existence results. If we apply Signorini contact conditions instead of normal compliance in this thermodynamic case, we have to impose auxiliary conditions into the input data of the coefficient function. Although the two PDEs system and physical situation seem to be more complicated than other contact problems, proving the convergence of numerical trajectories is less technically complicated. Indeed, in order to show the boundedness of contact forces in dynamic contact problems with Signorini contact conditions, we need to verify an additional condition, called the strong pointedness. This contact problem will play a fundamental role in being able to extend into more realistic contact models such as switches or actuators in Micro-electro-mechanical systems (MEMS). The V-shape electro-thermal actuator which is motivated by MEMS system is carried out in the paper [16]. In order to compute numerical approximations, we use the combination of a time discretization on the time interval and the Galerkin approximation over the spatial domain. This fully discrete numerical scheme leads us to establish the recursive formula where we can obtain the next time step numerical solutions per each time step. This paper is organized as follows. Section 2 introduces continuous formulations of our contact model. In Section 3, mathematical notations and background are illustrated. In Section 4, we set up a weak formulation which is equivalent to the original PDEs system. The main results of the existence are presented, as well. In Section 5, we present a detailed description of numerical schemes and the convergence theory, applying a time discretization into the weak formulation. In Section 6, the new form of energy balance is derived and the uniqueness of the weak solution is proved. In Sections 7 and 8, the fully discrete numerical schemes are proposed, numerical experiments are performed, and numerical results are presented. Section 9 concludes this with several important remarks. 2. Contact model Fig. 2.1 illustrates our physical situation in which a thermoviscoelastic rod fixed at the left wall stretches out, its right end touches the right wall, and shrinks back. The displacement and velocity of the rod are denoted by u = u(t , x) and v = ut (t , x), respectively, and its temperature distribution is denoted by θ = θ (t , x), where x ∈ [0, l] with its initial length l and the time variable t ∈ [0, T ] with the final time T > 0. The two walls are assumed to be kept at different temperatures. The right wall is assumed to be a stationary deformable obstacle located at x = ϕ > 0. The rod is allowed to penetrate the obstacle such that u(t , l) − ϕ ≤ ε for all t ∈ [0, T ], where a sufficiently small ε > 0 is related to wear and hardness of its surface. If ε is too big, the obstacles will disintegrate from a physical point of view. We assume that there are two given quantities, F = F (t , x), called a body force and Q = Q (t , x), called an internal heat source. As we mentioned above, we can expect four boundary conditions: two homogeneous Dirichlet boundary conditions on the left of the rod, x = 0, and two nonhomogeneous Neumann boundary conditions on the right of the rod, x = l. Our physical situation makes Neumann boundary conditions at x = l to be two contact conditions; one is the normal compliance and another is Barber’s heat exchange. The normal compliance condition is defined with the positive part function, denoted by (r )+ = max(r , 0) for r ∈ R. Since the reactive obstacle is penetrated as mentioned above, contact forces σ = σ (t ) can be characterized by the power

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

439

of the penetration

σ (t ) = λnc (u(t , l) − ϕ)p+ , where λnc > 0 and p ≥ 1 are called the stiffness of a reactive obstacle and the normal compliance exponent, respectively. Note that λnc is sufficiently large and p > 1 is helpful to regularize solutions. Unlike typical contact problems, in thermoviscoelasticity we require an extra condition, the radiation condition at the right free end. The boundary condition used in several papers [17,18,14] does not seem to be very realistic. Duvaut [19] studied a static problem with Signorini contact conditions where the thermal interaction is perfect when two surfaces are in contact and the two surfaces are perfectly insulated when they are separated. However, Duvaut’s thermal condition caused mathematical instability. Barber [9,10] improves the thermal condition in the limiting case; he postulates three boundary conditions which are perfect contact, imperfect contact, and non-contact. In Barber’s heat exchange condition, the heat exchange function κ = κ(r ) is ideally a step function which jumps down at r = 0. Then, the natural variable r = r (t ) depends on time t and is defined by r (t ) = ϕ − u(t , l) − σ (t ). The natural variable can be understood in this way: r > 0 when contact forces take place, and r < 0 when the rod is in separation from the reactive obstacle. Indeed, Barber’s three conditions are generalized to maximal monotone graph and the Neumann boundary condition (2.4) is set up in the paper [12]. In our contact model, the exchange function is assumed to be a continuous and decreasing function and particularly, drop down drastically around r = 0 to approximate the set valued graph. In addition to those conditions, the function y = κ(r ) has two horizontal asymptotes y = 0 and y = 1. A specific κ function will be sketched in Section 8 and used to implement numerical algorithms. We assume that rods are homogeneous and isotropic. The equations of linear thermoelasticity can consist of an equation of motion which is a hyperbolic type and an equation of energy which is parabolic type, as derived in Carlson’s paper [1]. All coefficients which appear in the two equations are the thermal conductivity k, the specific heat c, thermal expansion β , the elastic moduli λ, µ, and the density ρ . By the proper scale (see the Day’s book [4]), two equations are simplified with dimensionless constants a, b such that a2 =

θ0 β 2 (3λ + 2µ)2 , c (λ + 2µ)

b2 =

k2 ρ c (λ + 2µ)l2

,

where θ0 is the reference temperature. Adding the viscous term into the equation of motion, we are led to the following PDEs system which describes the physical setting explained above: b vt = uxx + α vxx − a θx + F

in (0, T ] × (0, l),

(2.1)

u(t , 0) = 0,

θt = θxx − a vx + Q in (0, T ] × (0, l), θ (t , 0) = 0 on (0, T ],

(2.3)

−λnc (u(t , l) − ϕ)p+ = ux (t , l) + α vx (t , l) − a θ (t , l) on (0, T ],

(2.4)

θx (t , l) = −κ(r (t ))θ (t , l) on (0, T ],

(2.5)

u(0, x) = u (x),

(2.6)

0

(2.2)

ut (0, x) = v (x), 0

θ (0, x) = θ (x) in (0, l), 0

(2.7)

where α > 0 is the viscous coefficient and all initial data are denoted by u (x), v (x) and θ (x). Note that subscripts denote partial derivatives with respect to the subscripted variables. 0

0

0

3. Mathematical preliminaries In this section, we introduce mathematical background and some notations. The spaces of solutions over the spatial domain (0, l) that we mostly deal with are based on parts of the Gelfand triples (see the book [20]): V ⊂ H = H ′ ⊂ V ′ , where all spaces which are separable Hilbert spaces are compactly and densely embedded. Let X be a Banach space. Then its dual space is denoted by X ′ and the duality pairing between X and X ′ is denoted by ⟨·, ·⟩X ′ ×X . We note that ⟨·, ·⟩ may be used if dual spaces are clear. It is well known from the extension of the inner product (·, ·)H ×H on V ′ × V over the Gelfand triples that for ς ∈ H and ξ ∈ V

⟨ς , ξ ⟩V ′ ×V = (ς , ξ )H ×V = (ς , ξ )H ×H .  l ∥ς ∥  0 ς ξ 1dx and the associated  norm 1 H = (ς , ς )H . For our convenience (·, ·)H is used instead of (·, ·)H ×H . Also let V = ς ∈ H (0, l) | ς (0) = 0 , where H (0, l) is a Sobolev  space. Similarly, its inner product and associated norm are defined to be (ς , ξ )V = (ςx , ξx )H and ∥ς ∥V = (ς , ς )V , In this problem, let H = L2 (0, l) with the inner product (ς , ξ )H =

respectively. Since solutions that we want to seek are vector-valued functions, we will extend those spaces with the Cartesian product. Throughout this paper, vectors and matrices are denoted by boldface letters. Let z(x) = (u(x), θ (x)). Then we define a couple of operators, L1 and L2 to be

L1 z = (uxx , θxx )T

L2 z = (α uxx − a θx , −a ux )T .

440

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

In order to show the results of existence, we list a couple of linear operators. Let V = V × V and V ′ = V ′ × V ′ . For z = (u, θ ) , w = (w, ϑ)T ∈ V , linear operators A1 , A2 : V → V ′ with appropriate boundary conditions are defined by

⟨A1 z, w⟩V ′ ×V := (L1 z, w)H = (ux , wx )H + (θx , ϑx )H , ⟨A2 z, w⟩V ′ ×V := (L2 z, w)H = α (ux , wx )H + a ((ux , ϑ)H − (θ , wx )H ) , where H = H × H. The elliptic self-adjoint operator A1 will play a significant role in proving the compactness and regularity of numerical trajectories for this dynamic contact problem. The operator  A2 will assist to improve the regularity of solutions. Since the vector-valued functions z : [0, l] → R2 and z ∈ H 1 0, l; R2 , its norm can be defined to be

 1/2 ∥z∥H 1 (0, l;R2 ) = ∥u∥2V + ∥θ ∥2V , where H 1 0, l; R2 is a Sobolev space for vector valued functions. It can also be seen easily that if z ∈ L2 0, l; R2 , then we









can have ∥z∥2L2 0, l;R2 = ∥u∥2H + ∥θ∥2H . We note that the duality pairing ⟨A1 z, z⟩1/2 is equivalent to the norm ∥z∥H 1 (0, l;R2 ) . ( ) In addition, the duality pairing ⟨A2 z, z⟩1/2 is equivalent to the norm ∥u∥H 1 (0, l) . In order to handle the normal compliance condition and Barber’s heat exchange condition properly in the abstract forms, we define the nonlinear operator P : V → V ′ which is decomposed by two operators P1 , P2 : V → V ′ ;

⟨P z, w⟩V ′ ×V = ⟨P1 u, w⟩V ′ ×V + ⟨P2 (u, θ ), ϑ⟩V ′ ×V , where

⟨P1 u, w⟩ = λnc (u(l) − ϕ)p+ w(l), ⟨P2 (u, θ ), ϑ⟩ = κ (u (l)) θ (l)ϑ(l) for fixed κ ∈ C (R) . Since the natural variable r is dependent of u(t , l) for all t ∈ [0, T ], the heat exchange function is defined by κ(r (u(t , l))) p with r (u(t , l)) = ϕ − u(t , l) − λnc (u(t , l) − ϕ)+ . Contact forces in Signorini contact conditions which are understood as the Dirac delta function δ can be approximated by the normal compliance, passing the limit λnc ↑ ∞. This implies that an obstacle is assumed to become perfectly rigid. This method combined with variational inequalities is a standard way to approach dynamic contact problems with Signorini conditions in the infinite dimension. On the other hand, another useful finite dimensional approach for Signorini type dynamic contact problems is to use complementarity conditions (CCs). In our future work, we will perform mathematical and numerical explorations on the comparison of the normal compliance condition to CCs. In this continuous dynamic thermoviscoelastic problem, for each t ∈ [0, T ] the energy function EC (t ) is defined to be EC (t ) := E [u(t ), v(t ), θ (t )]

=

1 2

 p+1 (v, v)H + (u, u)V + (θ , θ )H + 2λnc (u(l) − ϕ)+ .

(3.1)

In (3.1) the first term is called the kinetic energy, the second term called the elastic energy, the third term called thermal energy, and the last term called the energy by the normal compliance. The energy function defined in (3.1) plays an important role in not only showing the existence of solutions but also justifying numerical stability. However, we note that the energy function does not represent the physical energy of the thermoviscoelastic rods. Let W 1,s (0, T ; X ) be Sobolev spaces with real numbers s ≥ 1. Then if u ∈ W 1,s (0, T ; V ), then we define its norm by

∥u∥W 1,s (0, T ;V ) :=

  

T

∥u(t )∥sV + ∥˙u(t )∥sV dt

1/s

if 1 ≤ s < ∞,

0

 

ess sup (∥u(t )∥V + ∥˙u(t )∥V )

if s = ∞.

0≤t ≤T

We note that (˙) is the time derivative. When we prove the convergence of numerical trajectories, considering Hölder spaces, denoted by C q (0, T ; V ) with the exponent 0 < q ≤ 1 will be essential to show the compactness of numerical trajectories. Its norm is defined to be

∥z∥C q (0, T ;V ) = ∥z∥C (0, T ;V ) + sup s̸=t

∥z(t ) − z(s)∥V . |t − s|q

In addition to the spaces, interpolation spaces Vη are required. By interpolation theories (see the books [21, Section 22.6 Exercises] or [22] for instance), Vη can be identified with Sobolev spaces with the fractional index; Vη := H η (0, l) and V−η = Vη′ for −1 ≤ η ≤ 1.

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

441

4. Existence results In this section, we present the main results for the existence of solutions. Section 5 will support to prove the convergence of numerical trajectories, using a time discretization. Throughout Sections 4–6, all coefficients will be one, i.e., a = b = α = 1 for our convenience. In Lemma 1, we derive a weak formulation which is equivalent to the original PDEs system (2.1)–(2.5). A pair of solutions, (u, θ ) is written by u = (u, θ ) and a body force F and heat source Q are written by F = (F , Q ). Now, let v = (v, θ ) and w = (w(t , x), ϑ(t , x))T be vector valued functions. Lemma 1. Assume that for almost all t ∈ [0, T ], the solutions (u, θ ) satisfy Eqs. (2.1)–(2.5) and F ∈ V ′ are given. Then we can derive the following equivalent variational formulation: find u ∈ V for almost all t ∈ [0, T ] such that

⟨˙v + A1 u + A2 v + P u − F, w⟩ = 0 for all w ∈ V .

(4.1)

Its proof shall be skipped. The similar approach to Lemma 4.1 in the paper [23] can be applied. Theorem 2. Assume that F ∈ L∞ (0, T ; H ) and the initial data u0 , θ 0 ∈ V and v 0 ∈ H. If κ ∈ C (R) and 0 < u(t , l) − ϕ ≤ ε with sufficiently small ε > 0, then there exists unique weak solution (u, θ ) satisfying the variational formulation (4.1) such that u ∈ C (0, T ; V ) ∩ W 1,∞ (0, T ; H ) ∩ C ([0, T ] × C [0, l])

θ ∈ C (0, T ; V ) ∩ W

1,2

and

(0, T ; H ) ∩ C ([0, T ] × C [0, l]) .

Theorem 2 will be proved through many steps in Sections 5 and 6. 5. Time discretization and its convergence We begin this section by partitioning the time interval [0, T ] uniformly so that 0 = t0 < t1 < t2 < · · · < tk−1 < tk < tk+1 < · · · < tn = T , where ht = tk+1 − tk is the time step size with positive integers k ≥ 0. The numerical trajectory of the displacement, denoted by uht (t , x) will be a piecewise continuous linear interpolant such that uht (tk , ·) = uk (·) and uht (tk+1 , ·) = uk+1 (·) in the subinterval [tk , tk+1 ]. The numerical trajectory of the temperature, denoted by θht will be also a piecewise continuous linear interpolant. The numerical trajectories vht (t , x) and θ˙ht (t , x) for velocity v and the rate of change of the temperature θ will be piecewise constant interpolants such that vht (t , ·) = v k+1 (·) and θ˙ht (t , ·) = θtk+1 (·) for t ∈ (tk , tk+1 ]. In addition to them, we assume that a given vector-valued function (F (t , x), Q (t , x)) can be approximated by piecewise constant interpolants in the same way as vht and θ˙ht do. Throughout this paper, constants M , C > 0 refer to a quantity that does not depend on the parameter ht > 0 of the approximations but may be different in each occurrence. Using the midpoint rule and the Euler methods, we establish the numerical formulations over the time interval [0, T ] in the distributional senses:

   k+1   k+1   1 1 uk u v + uk + vk v k+1 − v k     + Fk , = − A1 k+1 − A2 k+1 −P k+1 uk , θ k+1 + θ k /2 θ + θk θ + θk − θk ht θ 2 2     1 v k+1 + v k 1 uk+1 − uk = k+1 k+1 k , 2θ −θ 2 ht θ 1



(5.1)

(5.2)

u0 ,

θ0 ∈ V, v0 ∈ H , (5.3)   T where Fk = F k , Q k . We note that for our convenience each time step solutions are written in terms of column vectors. Now we can define the energy function which is generated by numerical approximations at each time step t = tk : ESk := ES (tk )

=

1  2

vk , vk

 H

     + uk , uk V + θ k , θ k H + 2λnc (uk (l) − ϕ)p++1 .

(5.4)

The energy function enables us to investigate a possibility of being the boundedness of energy function. Its successful procedure will provide an evidence of numerical stability. Lemma 3. Assume that F : [0, T ] → H and F ∈ L∞ (0, T ; H ). If one time step numerical solutions uk , v k , θ k ∈ V × H × V satisfy the numerical formulations (5.1)–(5.3), there are uht ∈ C (0, T ; V ), vht ∈ L∞ (0, T ; H ), and θht ∈ C (0, T ; H ), independent of ht > 0.





442

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

Proof. The first component equation in (5.1) is multiplied by the first component in the vector equation (5.2) and the second component equation in (5.1) is multiplied by θ k+1 + θ k /2. Then it is easy to see that

     T  k+1  T  k+1 u − uk /ht uk+1 + uk /2 v + vk v k +1 − v k   ,  k+1 = − A1  k+1 , θ k+1 + θ k θ k +1 − θ k 2ht θ + θ k /2 θ + θ k /2 V ′ ×V H        T  k+1 T  k+1 k k k + 1 k u − uk /ht 1 u v +v v +v   − − P  k k+1 ,  k+1 , A2 u ,θ + θ k /2 θ k+1 + θ k θ k+1 + θ k 4 θ + θ k /2 ′ V ′ ×V V ×V      k+1 k k T v + v /2 F  . ,  k+1 + Qk θ + θ k /2 1



H

By doing several algebraic manipulations, we use the definitions of the operators to obtain the following equation with any ht > 0:

        1  k+1 2 v  − v k 2 + 1 θ k+1 2 − θ k 2 H H H H 2ht 2ht  k 2  1  k+1  2 1 1  k+1 2 2 u  − u  − v + v k V − θ k+1 + θ k V =− V V 2ht 4 4   k+1  1  k k+1  1  p k k − λnc (u (l) − ϕ)+ u (l) − u (l) + F ,v + vk H ht 2  2 1  k k+1  1  − κ uk (l) θ k+1 (l) + θ k (l) + Q ,θ + θk H . 4

2

(5.5)

Now, we change the previous equation into the following inequality:

               1   k +1  2 v  + uk+1 2 + θ k+1 2 − 1 v k 2 + uk 2 + θ k 2 + ht κ uk (l) θ k+1 (l) + θ k (l) 2 H H V V H H 2 2 4    k+1  ht  k k+1 ht  k k+1 p k k k ≤ F ,v +v H + Q ,θ + θ H − λnc (u (l) − ϕ)+ u (l) − uk (l) . 2 2

(5.6)

In order to get the boundedness of the energy, the last term of (5.6) receives a special treatment:

  − λnc (uk (l) − ϕ)p+ uk+1 (l) − uk (l)      = −λnc (uk (l) − ϕ)p+ + (uk+1 (l) − ϕ)p+ − (uk+1 (l) − ϕ)p+ × uk+1 (l) − ϕ − uk (l) − ϕ   ≤ −λnc −(uk (l) − ϕ)p++1 + (uk+1 (l) − ϕ)p++1 − λnc (uk (l) − ϕ)p+ (uk+1 (l) − ϕ) + λnc (uk+1 (l) − ϕ)p++1   ≤ −λnc −(uk (l) − ϕ)p++1 + (uk+1 (l) − ϕ)p++1 + C λnc ht ε p .

(5.7)

Note that the previous estimates are obtained when the right end of a rod penetrates an obstacle. We apply the Cauchy–Schwarz inequality into the first and the second term of the right side of the inequality (5.6) and use the trapezoidal rule into them. Then the telescoping sum allows us to get tk



ESk ≤ ES0 + ∥F ∥L∞ (0,T ;H )

0

  vh  dt + ∥Q ∥L∞ (0,T ;H ) t H

tk

 0

  θh  dt + C T εp . t H

Recalling the energy function defined in (5.4) and the inequality (5.8) and using Hölder’s inequality, we can have

 k 2  k 2 v  + θ  ≤ C1 H H

tk

 0

 2  2 vh  + θh  dt + C2 , t H t H

where C1 , C2 > 0 are appropriate constants. It follows from Gronwall’s inequality that

 2

max v k H ≤ M1 < ∞ k≥0

and

 2

max θ k H ≤ M2 < ∞. k≥0

It also follows from the inequality (5.8) and the energy function (5.4) that

 2

max uk V ≤ M3 < ∞, k≥0

as desired.



(5.8)

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

443

We note that vht ∈ L2 (0, T ; V ) can be easily shown through (5.5)–(5.8). In order to consider the compactness of numerical trajectory θht , we require that θht be in the space C (0, T ; V ), which will be shown in Lemma 4. Lemma 4. Under the same assumptions of Lemma 3 with sufficiently small ht > 0, there are θ˙ht ∈ L2 (0, T ; H ) and θht ∈ C (0, T ; V ) satisfying the numerical formulations (5.1)–(5.3). Proof. Using the second component in (5.2), we multiply both sides of (5.1) by 0, θtk+1 we have



T

T

or 0, θ k+1 − θ k /2 . Then,







  k+1 2       θ  ≤ − 1 θ k+1 2 − θ k 2 + 1 v k+1 + v k , θ k+1 t t H V V H +

2ht 1  2ht

2

P2 uk , θ k+1 + θ



 k

   , θ k+1 − θ k + Q k , θtk+1 H .

Multiplying by 2 ht > 0, the Cauchy–Schwarz inequality and Cauchy’s inequality with a sufficiently small fixed value ϵ > 0 can provide the following: 2 2 2ht θtk+1 H + θ k+1 V − κ uk





 2 ≤ θ k V  2 ≤ θ k 

V

   k+1 2   2 θ (l) + κ uk θ k (l)          + ht θtk+1 H v k+1 H + v k H + 2 Q k H         2 1  k+1 2 v  + v k 2 + 2 Q k 2 + ht 2ϵ θtk+1 H + . H H H 2ϵ 



Thus, we can obtain

  k+1 2   2   2 θ (l) − κ uk θ k+1 (l) + κ uk θ k (l) 2  2  2    2  2 ht  ≤ −κ uk+1 θ k+1 (l) + θ k V + 2 v k+1 H + v k H + 2 Q k H . 2ϵ Since the numerical trajectory uht and the heat exchange function κ are continuous, the telescoping sum provides the following: for any k ≥ 1  k         i 2 C2 T   tk  θ˙h (t )2 dt + θ k 2 ≤ κ r k θ k (l) 2 + C1 ht θ (l) + 2 , 1 − 2ϵ 2 t V H 2ϵ 0 i=1 2 2 ht 1 − 2ϵ 2 θtk+1 H + θ k+1 V − κ uk+1













where C1 , C2 > 0 are independent of ht > 0. Since θht is continuous over the domain [0, T ]×[0, l], we can get the following estimate: each k ≥ 1



1 − 2ϵ

2



tk

 0

       θ˙h (t )2 dt + θ k 2 ≤ κ r k θ k (l) 2 + C3 T + C2 T ≤ M < ∞, t V H 2ϵ 2

as C3 > 0 is independent of ht > 0. Now, the proof is complete.



In order to confirm that there are numerical trajectories uht , , vht , θht for any ht > 0, it is essential to show the uniqueness of the one time step solutions satisfying all numerical formulations (5.1)–(5.3) per time step t = tk . This can be executed, using the extension of Gelfand triple and Lax–Milgram Lemma (e.g., see the book [24, Section 6.2]). The proof is straightforward and thus shall be skipped. Based on results from the previous lemmas, we begin to consider the compactness of the solutions. Although there are many useful theorems to prove it, using Hölder continuity is one of excellent ideas to show the convergence successfully. The self-adjoint operator A1 can use eigenfunctions to define the fractional index of Sobolev spaces H γ (0, l). Let Vγ = H γ (0, l) be interpolation spaces with 0 ≤ γ ≤ 1. Then the following useful inequality (see the book [21, Section 22.6]) will be applied:



1−γ

∥u∥Vγ ≤ cγ ∥u∥H



γ

∥u∥V .

Then, it follows from Lemmas 3 and 4 that the boundedness of energy function can allow us to get to the following with the exponent p of Hölder space with q = 1 − γ

  uh (t ) − uh (s) ≤ C1,γ |s − t |q t t Vγ

and

  θh (t ) − θh (s) ≤ C2,γ |s − t |q , t t Vγ

where 0 ≤ γ < 1. Thus, we can say that uht , θht ∈ C q (0, T ; Vγ ). By the Sobolev embedding theorem, Vγ with 1/2 < γ < 1 is compactly embedded in C [0, l]. We can apply the Arzela–Ascoli theorem to show that Hölder space C q [0, T ] is compactly embedded in C [0, T ]. Therefore, there are subsequences, denoted by uht and θht , such that uht → u and θht → θ in C [0, T ] × C [0, l], as ht ↓ 0. Since vht ∈ L2 (0, T ; V ) ∩ L∞ (0, T ; H ) and θ˙ht ∈ L2 (0, T ; V ) ∩ L∞ (0, T ; H ), we can apply the Alaoglu theorem to show that there are subsequences vht ⇀∗ v and θ˙ht ⇀∗ θ in L2 (0, T ; V ) ∩ L∞ (0, T ; H ), as ht ↓ 0.

444

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

6. Energy balance Unlike the Signorini contact conditions, we use a modified energy function which includes the energy associated with the normal compliance. We choose a sufficiently smooth function w = ⟨v, θ ⟩T from the weak formulation (4.1) and integrate it over [0, t ] ⊂ [0, T ]. Then, deriving an energy balance form is a simple task. We can see that for p ≥ 1 t



⟨˙v + A1 u + A2 v + P u − F, w⟩ dτ  t t = (ux , vx )H + (θx , θx )H dτ (vt , v)H + (θt , θ )H dτ + 0 0  t  t p + (u(τ , l) − ϕ)+ v(τ , l) dτ (vx , vx )H + (θ , vx )H − (vx , θ )H dτ + 0 0  t  t κ(r (τ )) (θ (τ , l))2 dτ − + (f , v)H + (Q , θ )H dτ 0 0  t  t  1 d  d 1 p+1 2 2 2 ∥v∥H + ∥θ∥H + ∥u∥V dτ + = (u(τ , l) − ϕ)+ dτ 2 0 dτ p + 1 0 dτ  t  t  t ∥v∥2V + ∥θx ∥2H dτ + + κ(r (τ )) (θ (τ , l))2 dτ − (f , v)H + (Q , θ )H dτ .

0 =

0



0

0

0

Therefore, defining the total energy function E to be

 λnc p+1 ∥v(t )∥2H + ∥θ (t )∥2H + ∥u(t )∥2V + (u(t , l) − ϕ)+ , p+1

1

E (t ) =

2

it follows from the fundamental theorem that E (t ) = E (0) −

t



∥v(τ )∥2V + ∥θx (τ )∥2H dτ −

t



0

κ(r (τ )) (θ (τ , l))2 dτ + 0

t



(f , v)H + (Q , θ )H dτ .

(6.1)

0

The first and second integrals on the right side of (6.1) produce the energy loss, while the last integral describes the potential energy. In fact, deriving the energy balance may play a fundamental role in showing the boundedness of solutions in suitable spaces for the continuous case. Furthermore, we can prove the uniqueness of the solution (u, θ ), based on this energy balance. However, we need to be more careful than the derivation of the energy balance. Lemma 5. Under the same assumptions of Theorem 2, there is a unique solution (u, θ ) in the spaces presented in Theorem 2 satisfying the weak formulation (4.1). Proof. Suppose that there are two solutions (u1 , θ1 ) and (u2 , θ2 ) in the space S1 × S2 , where S1 and S2 denote the spaces of solutions used in Theorem 2. Let u = u1 − u2 , v = v1 − v2 , and θ = θ1 − θ2 . Since v ∈ L2 (0, T ; V ) and θ ∈ C (0, T ; V ), it follows that on the interval t ∈ (0, T ]

 0=

t

⟨˙v + A1 u + A2 v, v⟩ + ⟨P1 u1 − P1 u2 , v⟩ dτ +

0

t



⟨P2 (u1 , θ1 ) − P2 (u2 , θ2 ) , θ ⟩ dτ . 0

Then, we define a new energy function E to be E (t ) =

1 2

 ∥v(t )∥2H + ∥θ (t )∥2H + ∥u(t )∥2V .

By the similar way as we derived the energy balance, we can have E (t ) ≤ E (0) −

t



∥v∥2V + ∥θx ∥2H dτ + λnc ε p 0

t



|v(τ , l)| dτ + 0

t



|κ (u1 ) − κ (u2 )| (θ (τ , l))2 dτ . 0

Using the trace theorem [24, section 5.5] we can obtain E (t ) ≤ −

t



∥v∥2V + ∥θx ∥2H dτ + C1 0

t



∥v∥2V dt + 0

t



|κ (u1 ) − κ (u2 )| ∥θ∥2V dτ , 0

where C1 > 0 is constant. Then, since κ and u are continuous, we can choose t1 ∈ (0, T ] such that t1 < 1/C2 . Thus

∥θ (t1 )∥2H + ∥u(t1 )∥2V ≤ C2

t1

 0

∥θ∥2H dt + C1

t2

 0

∥v∥2V dt .

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

445

By the net change theorem, we have

∥θ (t1 )∥2H + ∥u(t1 )∥2V ≤ C

t1



∥θ ∥2H dt + ∥u∥2V dt ,

0

where a suitable constant C > 0 can be determined. Therefore it follows from Gronwall’s inequality that u = 0 and θ = 0 on the interval [0, t1 ]. We apply the same argument on the next intervals [t1 , t2 ] , [t2 , t3 ], etc. and thus the proof is complete.  7. Fully discrete numerical schemes The fully discrete numerical schemes are proposed, based on time discretization and the Finite Element Method (FEM) over the spatial domain [0, l]. As we did on the time interval [0, T ], the domain [0, l] will be uniformly partitioned; 0 = x0 < x1 < x2 < · · · < xj < xj+1 < · · · < xm = l, where hx = xj+1 − xj is the size of subintervals and m = l/hx is the number of subintervals. Now, we choose the finite dimensional space over [0, l]; Vhx = whx ∈ C [0, l] | whz (0) = 0 ,





where whx is apiecewise linear function. So we use the piecewise linear basis functions ψj (x) with  1 ≤ j ≤ m to approximate  the solutions uk , v k , θ k at each time step t = tk . Thus, the fully numerical approximations uht ,hx , vht ,hx , , θht ,hx will be written by the linear combination of basis functions ψj : uhx (tk , x) =

m 

ukj ψj (x) , vhx (tk , x) =

j =1

m 

vjk ψj (x) , θhx (tk , x) =

j=1

m 

θjk ψj (x) .

j=1

Here simpler notations uhx (tk , x) = uht ,hx (tk , x), vht ,hx (tk , x) = vhx (tk , x), and θhx (tk , x) = θht ,hx (tk , x) are used. Recalling the semi-discrete numerical formulations (5.1)–(5.2) and assuming that F = 0 for our simplicity, we establish the fully discrete numerical formulations: b

vhkx+1 − v,khx

=

ht

θhkx+1 − θhkx ht

vhkx+1 + vhkx 2

1 ukh+ x



= =

′′

 ′′ + ukhx 2

 k+1 ′′  k ′′ θhx + θhx 2 1 − ukhx ukh+ x

ht

+α

 k+1 ′′  k ′′ vhx + vhx 2

1 ukh+ x

 −a

′

 ′ − ukhx

ht

 −a

′  ′ θhkx+1 + θhkx 2

,

,

(7.1)

(7.2)

,

(7.3)

where (′ ) means the derivative with respect to x ∈ [0, l]. Using the boundary conditions (2.3) and (2.4) and the linear combination of the fully approximations, we multiply the basis function ψi with 0 ≤ i ≤ m and apply integration by parts to get to the following linear system in a matrix from:



2b

M+ 2

K+

2 a T L ht

 h  t 



1

2b

M− 2

 h t = 

1

α ht



2



2 a T L ht

α ht



− L b ht

K+



a

K

 K

M+

2b ht

M

0

    k+1  + I uk+1   θ 1

2

K

a



 k  uk    v   k b 1 M− K θ 2

L

ht



+ 0, . . . , 0, −λnc (ukm − ϕ)p+ , 0, . . . , 0, −

2

  1

2

κ(r k )θmk

T

,

where I := Iij ∈ R2m×2m contain almost all zeros except that I2m2m =

1 2

 T κ(r k ) and uk = uk1 , uk2 , . . . , ukm ∈ Rm , vk =

 k k T  T v1 , v2 , . . . , vmk ∈ Rm , θ k = θ1k , θ2k , . . . θmk ∈ Rm . The symmetric matrices M, K, L ∈ Rm×m used above are defined by  l  l  l M= ψj (x) ψi (x) dx, K= ψj′ (x) ψi′ (x) dx, L= ψj (x) ψi′ (x) dx. 0

0

0

446

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

In the fully discrete case, energy function for each time step t = tk is defined to be EFk := EF (tk )

 T   p+1  T 1  k  k T + λnc ukm − ϕ + , + uk K uk + θ k M θ k bv M v 2

=

(7.4)

where λnc is a fixed large number. Lemma 6. Assume that the fully discrete solutions uht ,hx , vht ,hx , θht ,hx satisfy Eqs. (7.1)–(7.3) with the boundary conditions. Then we obtain the following estimate: for any k > 0



k k T  ht    T   ι   T  T  ht   ι + v + vι−1 K (vι ) + vι−1 θ + θ ι−1 K θ ι + θ ι−1 4 ι=1 4 ι=1

EFk + α

k ht 

+



2 ι=1

 2  κ r ι−1 θmι + θmι−1 ≤ EF0 + C T .

(7.5)

Proof. From (7.1), using the linear combination of the fully discrete approximations and integration by parts, we can obtain b



m 

2ht

v

k+1 j

−

j=1

=−

−

v

k j





2ht

ukj +1

+

j =1



4

m 

m 

v



l

ψj ψi dx

m 

+

j=1

m 

v

v

k+1 i

+

ukj





l

ψj ψi dx ′

′

0

k j



 v

k i

m 

uki +1

−

i =1



l

ψj ψi dx ′

′

0

j =1

m  i=1

i=1

j=1 k+1 j

m 

0

j =1

1

α

m 

m 



m 

uki

i=1

v

k+1 i

+

i=1

m 

 v

k i

+ I,

i =1

where I =



a 2ht

m  j =1

θ

k+1 j

+

m  j =1

θ

k j





l

ψj ψi dx

0

′

m 

uki +1

i=1

−

m 

 uki

+

i =1

 ′   1  k+1 um + ukm ψm ψm ukm+1 − ukm

2ht

      α  k +1 a  k+1 + vm + vmk ψm′ ψm vmk+1 + vmk − θm + θmk ψm ψm ukm+1 − ukm . 4

2ht

We can substitute thematrices and vectors to simplify this. We substitute Eq. (7.3) again into the next to last term, we can  factor out ukm+1 − ukm and have

  T  T    T  T  1  k+1 b  k+1 v − vk M vk+1 + vk u + uk K uk+1 − uk =− 2ht 2ht   k+1 T  k T    T  T  α  k+1 a  k+1 k − v +v K v + v + θ + θ k L uk+1 − uk 4 2ht  ′  k+1  ′     1  k+1 k k + um + um ψm ψm + α vm + vm ψm ψm − a θmk+1 + θmk ψm ψm ukm+1 − ukm . 2ht Now, we multiply by ht > 0 and apply the natural boundary conditions:

  T  T    T  T  b  k +1 1 v − vk M vk+1 + vk = − uk+1 + uk K uk+1 − uk 2 2   a       T  T   p   α ht  k+1 T T − v + vk K vk+1 + vk θ k+1 + θ k L uk+1 − uk + − λnc ukm − ϕ + ukm+1 − ukm . 4 2 We also multiply out the vector sums and rearrange the terms and have the following simplified equation: b  k+1   k+1 T 1  k+1   k+1 T b  k   k T 1  k   k T v M v + u K u = v M v − u K u 2 2 2      T  T   T  T  α ht  k+1 a  k+1 − v + vk K vk+1 + vk + θ + θ k L uk+1 − uk 4 2

2

p    − λnc ukm − ϕ + ukm+1 − ukm .

(7.6)

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

447

Similarly, it follows from (7.2) that 1



m 

2ht

θ

k+1 j

−

m 

j =1

=−

−

1



4 a

θ

k j



θ

k+1 j

+

j =1



m 

2ht

ψj ψi dx

m 

θjk

θ

j=1

+

m 

θ

k+1 i

+

i=1





l

θ

k j

ψj′ ψi′ dx

m 

θ

k+1 i



l

ψj ψi dx ′

0

 θ

k i

+

i=1



j =1

m  i =1

0

j=1 k+1 j

m 

0

j =1 m 



l

m  i=1

m 

 θik

i =1

uki +1

−

m 

 uki

i=1

+

   1  k+1 θm + θmk ψm′ ψm θmk+1 + θmk . 4

We substitute the matrices, vectors and use the heat exchange condition to get

    T 1  k+1 T 1  k+1 θ − θ k M θ k+1 + θ k + θ + θ k K θ k+1 + θ k 2ht 4 2   T 1    a  k+1 θ + θ k L uk+1 − uk − κ r k θmk+1 + θmk − 2ht 2 and thus we multiply by ht to obtain

  T  T  1  k+1   k+1 T 1  k   k T ht  k+1 − θ M θ θ M θ =− θ + θ k K θ k+1 + θ k 2 2 4   k+1 T  k T  ht  k   k+1 2 a  k+1 k +θ L u − u − κ r θm + θmk . − θ 2 2 Now, we add Eqs. (7.6) and (7.7) together to obtain

(7.7)

1  k+1   k+1 T 1  k+1   k+1 T b  k+1   k+1 T v M v + u K u + θ M θ 2 2 b  k   k T 1  k   k T 1  k   k T = v M v − u K u + θ M θ 2 2 2 T  T  ht  k+1 T  T      α ht  k+1 − v + vk K vk+1 + vk − θ + θ k K θ k+1 + θ k 4 4   k+1 T  k T  a  k+1   T  T  a  k+1 k + θ +θ L u − u θ + θ k L uk+1 − uk − 2 2  k p  k+1 2  ht  k   k+1 k − λnc um − ϕ + um − um − κ r θm + θmk . 2

2

Since the terms involving the L matrix is canceled out and there are the three perfect terms, we can have 1  k+1   k+1 T 1  k+1   k+1 T b  k+1   k+1 T v M v + u K u + θ M θ 2 2 2 b  k   k T 1  k   k T 1  k   k T = v M v − u K u + θ M θ 2 2 2     T  T   T  T  ht  k+1 α ht  k+1 − v + vk K vk+1 + vk − θ + θ k K θ k+1 + θ k 4 4   k p  k+1  ht  k   k+1 k 2 − κ r θm + θm − λnc um − ϕ + um − ukm . (7.8) 2 Modifying the last term of (7.8), the energy function (7.4) defined in the fully discrete case provides the following inequality:

α ht 

 T  + vk + C ht ε p 4   T  T  ht  k   k+1 2 ht  k +1 − θ + θ k K θ k+1 + θ k − κ r θm + θmk .

EFk+1 = EFk −

 

vk+1 + vk K

vk+1

T

4 2 Thus, the telescoping sum can allow us to get to the estimate (7.5). The proof is complete.



We recall the form of energy balance (6.1) in the continuous case. The second plus through the last term in (7.5) is called an energy loss in the fully discrete case which is perfectly due to the viscosity and the heat conduction and energy associated with heat exchange condition. Even if for sufficiently small time step size ht > 0 each of summation terms in (7.5) are almost close to each of integrals, there is a difference, i.e., C T between those energy balance forms. The gap seems to be a little big but it does not practically, which will be supported by numerical results in Section 8. The result of Lemma 6 will play a fairly reasonable role in verifying stability in the numerical solutions of the PDEs system.

448

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451 Table 1 Data used in numerical experiment. ht

hx

l

T

a

b

α

ϵ

ϕ

p

c

10−6

10−3

1

0.3

0.03

10−4

10−6

10−8

1.1001

4

100

(A) Energy.

(B) Energy + energy loss. Fig. 8.1. Energy.

(A) Initial temperature.

(B) Heat exchange function κ . Fig. 8.2. Data for the thermal effect.

8. Numerical results and discussion Table 1 summarizes the input parameters and their values used in our numerical experiments. We begin with the deformation of the rod as 0.1 × x where x is the unstretched position of the rod. This places the right end of the rod very near the obstacle, in order to see its immediate contact. In addition, the rod moves from rest initially and thus the initial velocity v 0 (x) = 0. Matlab TM utility is able to call the sparse function which are efficiently used to compute numerical solutions of the linear system each time step. In the linear system, the matrix size is 2000 × 2000 and numerical approximations are computed by the direct method at each time step, running through 300,000 time steps. The average time of computing one time step solution is approximately 0.07 s, which shows that the numerical approximations are reasonably stable. From our numerical performance, we obtain the maximum penetration of the deformable obstacle, ε = 0.007. According to the theoretical result of Lemma 6, we expect at least the energy function to be bounded. The two graphs of the energy function are shown in Fig. 8.1. The graph of the energy function in (A) shows that the energy defined in (7.4) decreases because of energy loss. The remainder C T on the right side of (7.5) does not seem to have any effect on getting monotonic decreasing. The graph including the energy loss terms in (B) shows the perfect boundedness over the whole period [0, 0.3]. Since the thermal effect is involved in the rod, we consider an initial temperature function for the rod. In order to obtain an initial temperature θ 0 , we solve the ODE which arises from the natural boundary conditions (2.5) and its solutions satisfy the essential boundary condition in (2.3). Here is one possible function θ 0 (x) = e−κ(r )(x−1) − eκ(r ) (x − 1)2 which is displayed in the picture (A) of Fig. 8.2. The picture (B) of Fig. 8.2 shows the heat exchange function κ(r ) = 1 + π2 arctan (−c · r ) which can approximate Barber’s idealized step function κ , as c ↑ ∞. We graph the velocity and temperature distribution every 1000 time steps to observe the behavior of the rod. For comparison, we subdivided the time interval into ten local time subintervals.

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451 25

1.2

15

20

1

10

449

1.2 1

15 0.8

10 5

0.6

0

0.4

5

0.6 0

–5

0.4 –5

0.2

–10 –15 –20

0

0.2

0.4

0.6

0.8

1

–10

0

–0.2

–15

–0.2

0

0.2

0.4

0.6

0.8

1

10

0.8

1

0.8

6

0.7

0.7

4

0.6

2

0.5

0

0.4

–2

0.3

–4

0.2

0.1

–6

0.1

0

–8

0.3

0

0.2 –5

0.4

0.6

0.8

1

–0.1

0

0.2

0.4

0.6

0.8

1

(C) [0.06, 0.09).

0

0.2

0.4

0.6

0.8

1

–10

0

0.2

0.4

0.6

0.8

1

0 0

0.2

0.4

0.6

0.8

1

–0.1

(D) [0.09, 0.12).

6 4

0.8

5

0.8

0.7

4

0.7

0.6

3

0.4 0

0.3 0.2

–2 –4

0.4

0.6

0.8

1

0.6

2

0.5

2

0.2

0.6

8

0.4

0

0.4

0.8

0.5

5

–6

0.2

0.9

0.6

0.2

0

(B) [0.03, 0.06).

15

0

0.2

0

(A) [0, 0.03).

–10

0.8

1

0.5

0

0.4

–1

0.3

–2

0.1

–3

0

–4

–0.1

–5

0

0.2

0.4

0.6

0.8

1

(E) [0.12, 0.15).

0.2 0.1 0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

(F) [0.15, 0.18).

4

0.7

3

0.6

2.5

0.7

2

0.6

1.5

2 1

0.5

1

0.5

0.4

0.5

0.4

0

0 0.3

–0.5

0.3

0.2

–1

0.2

–1 –2

–1.5

–3

0.1

–4

0

0

0.2

0.4

0.6

0.8

1

0.1

–2 0

0.2

0.4

0.6

0.8

1

(G) [0.18, 0.21).

–2.5

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

(H) [0.21, 0.24).

1.5

0.7

0.4

0.7

1

0.6

0.3

0.6

0.2

0.5

0.5

0.5

0.1

0.4

0.4 0

0 0.3 –0.5

–0.2

–1

0.1

–0.3

–1.5

0

0

0.2

0.4

(I) [0.24, 0.27).

0.6

0.8

1

0.3

–0.1

0.2

0

0.2

0.4

0.6

0.8

1

–0.4 0

0.2 0.1 0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

(J) [0.27, 0.3]. Fig. 8.3. Velocity and temperature.

As seen in the first pictures in (B) to (H) of Fig. 8.3, the velocity curves create more of a pattern over their intervals. The first graph in (A) is of the velocity over the first two largest contact points. One can see that the subintervals displayed on the velocity reflect the behavior of contact forces. The last two graphs (I) and (J) show the smallest changes of velocity and

450

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

(A) With thermal effect.

(B) Without thermal effect. Fig. 8.4. Contact forces.

imply that the motion of the rod will stop soon. With contact, the velocity is mostly negative, meaning that the rod is being pushed back by the obstacle, but there are times when the velocity of the right end of the rod is still in the positive direction because the contact forces have not yet overcome the forward momentum of the rod. Without contact, the velocity is mostly positive, meaning that the rod is expanding toward the obstacle. It is evident that the average velocity is decreasing over time, since the scale of the graphs is getting smaller. The second pictures in each subfigure of Fig. 8.3 show the temperature distribution of the rod over each subinterval, beginning with the initial contact. As we observe the second picture in (A), a large contact force causes the wider range of change of temperature: especially, the right end of the rod shows a more drastic change in temperature when there is contact, often having a much lower temperature than the surrounding area. As time passes, a steady temperature distribution is reached as observed in the picture (J). Although another numerical simulation with the different initial temperature is not displayed in the paper, an interesting phenomenon is observed. When we use the initial temperature θ 0 (x) which does not satisfy the homogeneous boundary condition at the left end of the rod, the effect of the homogeneous boundary condition gradually appears and eventually the temperature of the left end becomes zero, i.e., the left end is insulated. Unlike velocity around x = 0 shown in all pictures in Fig. 8.3, it is observed that there are considerably more vibrations near x = 0. The contact forces of the thermoviscoelastic rod and almost purely elastic rod are graphed in (A) and (B) of Fig. 8.4, respectively. There is a relatively very small contact force initially due to the zero initial velocity. So we note that the initial contact force is not visible in (A) of Fig. 8.4. The magnitude of velocity is very large on this interval, and individual curves show many irregular waves without any pattern right after the second contact. The velocity at the endpoint seems unstable because it is immediately following the first contact. From our physical configuration the contact forces are in the negative direction. After the thermoviscoelastic rod springs back from the largest contact force, it does not loose contact with the obstacle again until close to T = 0.26, but it oscillates back and forth while contacting the obstacle. Due to the initial position of the rod and the obstacle, the first instance of contact occurs relatively quickly, but it cannot be seen on the graph of the picture (A) because of the scale. The magnitude of contact forces is compared in the two cases. The two pictures show the numerical evidence on how a nearly perfectly elastic rod would behave in the absence of the thermal conditions. We tested the same data and found that the rod does not contact the obstacle at all during the time (0, 0.3]. By changing the position of the obstacle to ϕ = 1.09, the rod experiences an initial contact force at time t = 0 and continues to experience contact forces as time passes, as shown in the picture (B). The initial contact force is about four times larger than the largest contact force in the thermoelastic case. The other contact forces are comparable, but in the nearly elastic case, the contact forces do not decrease as quickly as in the thermoelastic case. Therefore, we can conclude from the numerical results that the thermal effect makes the contact forces more regular and the rod to stop touching the obstacle more quickly. We note that the viscosity produces more singularity of contact forces. 9. Conclusions This work demonstrates an efficient numerical scheme for solving the PDEs system of thermoviscoelastic dynamic contact. The contact conditions are considered by Barber’s heat exchange condition and normal compliance. The convergence of numerical trajectories cannot be shown without using the viscosity. In addition, the heat exchange function which becomes a set valued type of function from Barber’s idealized condition causes our convergence theory to be unsuccessful. Further investigation on the issues will be carried out in future work. Mathematical proofs on whether the thermal effect makes an influence on the behavior of contact forces or not would be an interesting and important task for a research community of contact problems. This contact model can be switched into Signorini contact conditions and then the existence of solutions can be shown, passing the limit of normal compliance coefficient (e.g., see papers [25,26]). If a deformable obstacle is assumed to be a Winker foundation, maybe the thermal effect causes more technical complexity to show

J. Ahn, N. Clark / Nonlinear Analysis: Real World Applications 22 (2015) 437–451

451

the relation between the impact period and the foundation’s elastic coefficient. The related work can be found in the papers [27,28]. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

D.E. Carlson, Linear Thermoelasiticity, in: Handbuch der Physik. Bd., vol. Vla/2, Springer-Verlag, Berlin, 1972. M. Sofonea, M. Shillor, J.J. Telega, Models and Analysis of Quasistatic Contact, in: Lecture Notes in Physics, Springer, Berlin, 2004. M. Frémond, Non-Smooth Thermomechanics, Springer, Berlin, 2002. William Alan Day, Heat Conduction within Linear Thermoelasiticity, vol. 30, Springer-Verlag, New York, 1885. A. Klarbring, A. Mikelić, M. Shillor, Frictional contact problems with normal compliance, Int. J. Engng. Sci. 26 (8) (1988) 811–832. Kenneth L. Kuttler, Meir Shillor, Dynamic contact with normal compliance wear and discontinuous friction coefficient, SIAM J. Math. Anal. 34 (1) (2002) 1–27. J.A.C Martins, J.T. Oden, Models and computational methods for dynamic frcition phenomena, Comput. Meth. Appl. Engng. 52 (1985) 527–634. J.A.C Martins, J.T. Oden, Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws, Nonlinear Anal. TMA 11 (3) (1987) 407–428. J.R. Barber, Contact problems invloving a cooled punch, J. Elasticity 8 (4) (1978) 409–423. J.R. Barber, J. Dundurs, M. Comninou, Stability consideration in thermoelastic contact, J. Appl. Mech. 47 (1980) 871–874. K.T. Andrews, Andero Mikelić, P. Shi, M. Shillor, S. Wright, One-dimensional thermoelastic contact with stress-dependent radiation condition, SIAM J. Math. Anal. 23 (6) (1992) 1396–1416. K.T. Andrews, P. Shi, M. Shillor, S. Wright, Themroelastic contact with Barber’s heat exchange conditions, Appl. Math. Opt 28 (1) (1993) 11–48. Kevin T. Andrews, Andro Mikelic, Peter Shi, Meir Shillor, Steve Wright, One-dimensional themroelastic contact with a stress-dependent radiation condition, SIAM J. Math. Anal. 23 (6) (1992) 1393–1416. P. Shi, M. Shillor, S. Wright, A quasistatic contact in thermoelasticity with radiation conditions for the temperature, J. Math. Anal. Appl. 172 (1993) 147–165. Peter Shi, Meir Shillor, Xiu-Lin Zou, Numerical solutions to one-dimensional problems of themroelastic contact, Comput. Math. Appl. 22 (10) (1991) 65–78. Jeongho Ahn, Kenneth Kuttler, Meir Shillor, Existence and simulations of a dynamic thermoviscoelastic rod and beam system, in progress. R.P. Gilbert, P. Shi, M. Shillor, A quasistatic contact problem in linear thermoelasticity, Rend. Mat. 7 (10) (1990) 785–808. P. Shi, M. Shillor, Uniqueness and stability of the solution to a thermoelastic contact problem, European J. Appl. Math. 1 (1990) 371–387. G. Duvaut, Free Boundary Problems Connected with Thermoeelasticity and Unilateral Contact, in: Free Boundary Problems: Proceedings of a Seminar held in Pavia, vol. II, Instituto nazionale, Roma, 1979. J. Wloka, Partial Differential Equations, Cambridge University Press, 1987. Kenneth L. Kuttler, Modern Analysis, CRC Press, Boca Raton, FL, 1998. H. Triebel, Interpolation Theory, Function Spaces, Differential Operators, North Holland, Amsterdam, 1978. Jeongho Ahn, David E. Stewart, A viscoelastic Timoshenko beam with dynamic frictionless impact, Discrete Contin. Dyn. Syst. Ser. B 12 (1) (2009) 1–22. Lawrence C. Evans, Partial Differential Equations, in: Graduate Studies in Mathematics, vol. 19, AMS, Providence, Rhode Island, 1998. Jeongho Ahn, Kenneth L. Kuttler, Meir Shillor, Dynamic contact of two gao beams, Electron. J. Differential Equations 2012 (194) (2002) 1–42. Kenneth L. Kuttler, Meir Shillor, Vibrations of a beam between two stops, Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms 8 (1) (2001) 93–110. R.J. Knops, Piero Villaggio, On oblique impact of a rigid rod against a winkler foundation, Continum Mech. Thermodyn. 24 (2012) 559–582. Peter Shi, The restitution coefficient for a linear elastic rod, Comput. Model. 28 (4–8) (1998) 427–435.