NUMERICAL TREATMENT OF CONTACT PROBLEMS WITH THERMAL EFFECT

. The paper deals with the formulation and the ﬁnite element approximation of a quasi-static thermoviscoelastic problem which describes fric- tional contact between a deformable body and a rigid foundation. The contact is modeled by normal damped response condition whereas the friction is de- scribed by the Coulomb law of dry friction. The weak formulation of the model consists of a coupled system of the variational inequality for the displacement and the parabolic equation for the temperature. The main aim of this paper is to present a fully discrete scheme for numerical approximation together with an error estimation of a solution to this problem. Finally, computational simulations are performed to illustrate the mathematical model.


1.
Introduction. Modeling of contact between deformable bodies has great practical value. These problems occur in many different settings. Some examples from automotive industry are contact between brake pads and rotors or between pistons and cylinders.
The classical approach to solve presented problem requires the use of the theory of variational or hemivariational inequalities. It was developed and presented in [6]. The adaptation of this theory to solve mechanical contact problems can be found for example in [7]. A good example of monographs summarizing recent progress in this field is [5].
In many cases in the field of contact mechanics existence and uniqueness of solution to particular problem was proved, but numerical estimation of finite element method error was not conducted. There are some examples of scientific papers and publications that provide this estimation for particular cases, for instance [4], [2] and [1]. Still, a wide variety of problems is not covered.
This article presents numerical treatment of quasi-static problem with thermal effects. We provide numerical error estimation and present computational simulations for an example of stated problem.
The body is clamped on Γ D with a displacement equal to 0. The contact with the foundation occurs on Γ C . This may cause friction and damping response depending on introduced conditions. Moreover, the force f N acts on the part Γ N of the boundary and the force f 0 acts in Ω. The temperature of the body is equal to 0 on Γ D and Γ N , and heat exchange with its foundation Γ C may occur. The body movement and internal heat sources may also cause change of its temperature. Forces, factors modeling temperature and contact conditions may be time dependent. We are interested in the body displacement and the temperature in the time interval [0, T ], with T > 0.
Let us denote Q = Ω × (0, T ), Σ D = Γ D × (0, T ), Σ N = Γ N × (0, T ), Σ C = Γ C × (0, T ), S d = R d×d sym and ν the outside normal vector to Γ. We assume that the boundary Γ is Lipschtz continuous, and therefore ν exists a.e. on the boundary. Indexes i and j run from 1 to d and the index after a coma represents the partial derivative with respect to the corresponding component of the independent variable. Summation over repeated indices is implied. We denote the divergence operator by Div σ = (σ ij,j ). We also define the linearized (small) strain tensor with dependence on displacement Let u ν = u · ν and σ ν = σν · ν be the normal components of u and σ, and let u τ = u − u ν ν and σ τ = σν − σ ν ν be their tangential components. Now let us introduce the mathematical model of quasi-static contact problem with temperature.
Problem P . Find a displacement field u : Q → R d , a stress field σ : Q → S d and a temperature θ : Q → R such that Div Here, equation (1) represents a short memory viscoelastic constitutive law, where the operators A and B are viscosity and elasticity operators, respectively, and C M = (c ij ) is the expansion tensor describing the thermal influence. Equilibrium equation (2) corresponds to the assumption that acceleration of the body is negligibly small. Equation (3) represents the fact that body is clamped on part of its boundary and (4) represents outside forces. Equation (5) describes the damping response of the foundation, whereas the friction is modeled by equation (6). Equation (7) represents Fourier's law of heat conduction, where K = (k ij ) is the thermal conductivity tensor. In general, the Fourier law is of the form: c p ρ θ −div(K∇θ) = q and takes into account the thermal mass of the material c p ρ, but we normalize it and include the influence of body movement to obtain (7). The boundary condition on the temperature is given by (9), where θ f is the foundation temperature and k f is the heat transfer coefficient which describes the heat exchange between the considered domain (structure) and the environment.
3. Variational formulation. We recall that C([0, T ]; X) is the space of continuous functions from [0, T ] to X. We also use the standard notions of Lebesgue and Sobolev spaces (L 2 (Ω; R d ), L 2 (Ω; S d ) and H 1 (Ω), respectively) needed for the weak formulation. We consider the following Hilbert spaces respectively, and corresponding norms · X with X being H, H, H 1 , H 1 , V , E and L 2 (Ω). Moreover, we denote by ·, · Y * ×Y the duality paring between a dual space Y * and Y , with Y being V and E. Now we define some operators and functions in order to present variational formulation of Problem P . Let A, B : Using the standard procedure and Green's formula we obtain the variational formulation of Problem P in the following form We present the hypotheses on the data of Problem P .

H(K)
: the thermal conductivity operator K : Ω → S d is such that H(p ν , p τ ) : the normal compliance and friction bound functions p e : The body forces and surface traction, the heat sources density, the boundary thermic data and the initial data satisfy Finally, we need the following smallness assumption Under the above assumptions we have the following existence and uniqueness result for Problem P V .
Since the main result of the paper is the numerical analysis, we present only sketch of the proof of the existence of the unique solution to Problem P V .
Step 3. The operator Λ : , Theorem 2.9). Hence, u η * and θ η * being the unique solution to Problems P 1 η and P 2 η , respectively, are the solution to Problem P V . 4. Numerical scheme. In order to simplify notation we introduce the velocity variable v = u . Then for all t ∈ [0, T ] we have Let V h ⊂ V and E h ⊂ E be two families of finite dimensional subspaces with a discretization parameter h > 0. We introduce the time step k = T /N for N ∈ N, N > 0 and t n = nk, n = 0, 1, . . . , N . We also use notation g j = g(t j ) for any g ∈ C([0, T ]; X), where X is any introduced function space. We now present the following fully discrete scheme We remark that by a discrete analogue of Theorem 3.1 and a mathematical induction argument, the fully discrete scheme has a unique solution. Now we recall the Gronwall inequality in the following form (see [4], Lemma 7.26), which we use later. The main result of this section is the following theorem. and Proof. Throughout this proof we denote as c > 0 a generic constant (value of c may differ in different equations) and numerical solution errors for any j ∈ {1, . . . , N } by Let us fix j ∈ {1, . . . , N }. We start with the estimate of ε j . Since θ j ∈ L 2 (Ω), from equations (12) and (15), we obtain for all Writing then .

From the ellipticity of
.
From the fact that K and C 2 are Lipschitz continuous and from the Schwartz inequality, we get for all positive 1 , 2 , 3 .
Let us recall that from assumption (H 1 ), we have m A − m j > 0. Then, taking .
We sum this inequality from j = 1 to n for n ∈ {1, . . . , N } and take m ∈ {1, . . . , n} such that ε m L 2 (Ω) = max 1≤j≤n ε j L 2 (Ω) . Then, we multiply both sides by 2k and take into account the fact that T = N k to obtain Now, using the elementary inequality (a + b Then, from (18) we get Let us now fix j ∈ {1, . . . , N } and estimate e j . From (11) and (14) with w = v hk j and w h = w h j at t = t j , we obtain is given by (17). Because operators A, B and C 1 are Lipschitz continuous and j satisfies Summing this inequality from j = 1 to n for n ∈ {1, . . . , N }, multiplying it by 2k and adding to (19) we obtain We now use this result for n and n = m, and take into account the fact that Let us now estimate u j − u hk j V for j ∈ {1, . . . , N }. We observe Moreover, we get and from (21) we obtain Hence, Using (22) to estimate the first term on the right hand side of (20) and adding obtained inequality to (22) we get Since T = kN we obtain Then, applying the Gronwall inequality (cf. Lemma 4.1) to proves the theorem. We employ the Kelvin-Voigt short memory viscoelastic law for the isotropic body. The viscosity operator A and the elasticity operator B are defined by A(τ ) = 2φτ + ξtr(τ )I, τ ∈ S 2 , B(τ ) = 2µτ + λtr(τ )I, τ ∈ S 2 .
Here I denotes the identity matrix, tr denotes the trace of the matrix, λ and µ are the Lame coefficients whereas ξ and φ represent the viscosity coefficients, λ, µ, ξ, φ > 0. In our simulations we take We use spaces V h and E h of continuous piecewise affine functions as families of approximating subspaces. Numerical values denote temperature in a nearby vertex, the time step is equal to 0.01 and Figures 1, 2, 3 and 4 present outputs after 25, 50, 75 and 100 iterations, respectively. In our simulations we investigate the influence of heat sources represented by function q. We remark that heat sources cause the temperature of the body to change periodically. When we heat the body up, it increases its volume (Figure 1), whereas when we cool the body down, it returns to its original position ( Figure 2) and then it is squeezed (Figure 3). Finally, if we heat the body up once again, it returns to its original position (Figure 4). This behaviour corresponds with prescribed coefficients of operator C M .