Finite element treatment of two-dimensional thermoelastic wear problems

The present paper concerns the numerical treatment of thermoelastic wear problems. The governing equations of thermoelasticity coupled to Signorini contact, Coulomb's friction and Archard's wear are formulated as a system of discrete equations. This equation system is solved, using a Bouligand differentiable Newton method, for five problems of didactic nature. © 1999 Elsevier Science B.V. All rights reserved.


Introduction
In the present paper, a Newton method for thermoelastic wear problems is suggested. The continuum thermodynamical wear model presented in [20] is considered for a two-dimensional thermoelastic body which is unilaterally constrained by a rigid flat tangential moving support. In addition, it is assumed that no heat can flow through the support. For this system, starting from the governing equations, a system of finite element equations, which can be associated with an augmented Lagrangian formulation (see e.g. [1,11]), is derived. This system of equations is numerically solved for five problems of didactic nature using a Bouligand differentiable Newton method given by Pang [16].
Augmented Lagrangian formulations of two-and three-dimensional elastic wear problems were solved in [21,22] using Pang's Newton method. Furthermore, an extensive comparison of Pang's algorithm with an interior point method and a method commercially available in ABAQUS was performed for two-dimensional friction problems in [4]; numerical tests showed that the Newton method is both effective and robust. In the present paper, thermal properties are added to the augmented Lagrangian formulations, discussed in the papers above, in order to achieve a robust and efficient algorithm for thermoelastic wear problems. The method achieved is a direct method without any splitting between mechanical and thermal subproblems (see e.g. [18]). Numerical tests have proven that the method works viably. However, a crucial detail for making the method effective is to eliminate the temperature by algebraical manipulations, such that the temperature is not explicitly involved in the Newton step, but rather appears as an inner variable in the algorithm. Otherwise, the temperature field will not be solved accurately, which in turn, of course, would affect the quality of the mechanical part of the solution.
During the past years, a certain degree of progress has been made concerning the finite element analysis of thermomechanical contact problems. ~ Zavarise et al. [29] formulated a finite element model for thermal contact conduction; frictional heat generation was briefly discussed using a frictional tangential stiffness term and Concerning the analysis of thermomechanical contact problems of simple geometries such as half-planes, see [2]. 0045-7825/99l$see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S0045-7825 (98)00392-2 1 additional empirical coefficients. In [26], frictionless thermal contact conduction problems were considered using a penalty formulation, where the penalty terms depended on micro-phenomenological assumptions. Thermoelastic contact problems including both thermal contact conduction and frictional heat generation were formulated and solved in [9] by splitting the problem in one mechanical part and one thermal part; the mechanical problem was solved using a frictional contact algorithm suggested by Klarbring and Bjrrkman [10]. Wriggers and Miehe [25] treated frictional contact problems including thermal contact conduction and frictional heat generation using a normal compliance law and a radial return approach for the friction. In [28], augmented Lagrangian techniques were used for thermomechanical frictional contact problems described with microscopical laws; unfortunately, as it seems (see also [27]), the theoretical contact formulation is for bilateral contact and not for the proper unilateral formulation; in addition, a simplification of Coulomb's law of friction is used and the frictional heat generation is excluded. Recently, finite element formulations of thermomechanical frictional sliding in large deformations were given by Laursen and Oancea [12,15].
In this paper, we focus on the finite element treatment of frictional heat generation. Therefore, thermal contact conduction is not considered; of course, it is straightforward to include heat transfer through the contact interface (see [20]). Actually, wear heat generation is also included in the thermomechanical model presented herein. However, for physically realistic wear coefficients, the wear heat generation can be neglected in comparison to the frictional heat generation. Still, the internal state variable measuring the wear gap has a major influence on the frictional heat generation and vice versa. This fact is illustrated by numerical examples. For instance, it is numerically shown that two frictional hot spots appearing in a frictional analysis are removed if wear also is included in the constitutive settings.
The contents of the present study is as follows: in Section 2, the governing equations for a unilaterally constrained thermoelastic body are given. In Section 3, the governing equations are put together to an initial boundary value problem, which in turn is discretized both in space and time in order to reach a system of finite element equations. The system of equations obtained is considered in Section 4, where the adoption of Pang's Newton algorithm also is discussed. In Section 5, five numerical problems are studied. The problems differ only in the constitutive settings. That is, five different combinations of elasticity, thermoelasticity, Signorini contact, Coulomb's friction and Archard's wear are investigated for a particular initial boundary value problem. Finally, in Section 6, concluding remarks are given.

Governing equations
Let us consider a two-dimensional body in a quasi-static small displacement setting, occupying a domain /2 @ ~t2 with a piecewise smooth boundary a~O, see the heat flux vector by q, the Cauchy stress by 0", the infinitesimal strain tensor by ~ and the location of a point by x. Vectors and tensors are represented by components given in a Cartesian basis (e l, e2) which is defined in Fig. 1. Moreover, the summation convention is applied to repeated indices and Kronecker's delta 6,j is used. The boundary 0/2 is dived in three disjoint parts/7, F, and ~. The body is subjected to prescribed tractions t on /7, and zero displacements on F . It is further assumed that the body is unilaterally constrained on /7,. by a rigid support in the direction n,.. No heat transfer can take place through the support. Furthermore, the support is flat ended with a normal n s = e 2 and subjected to a tangential displacement ~ = sCe~, where s c = ((t) is a function of time t. An important consequence of the small displacement assumption is that n, = -n , . This fact is utilized in the work by setting n, = -e 2.

Equilibrium equations and principles of thermodynamics
Neglecting inertial forces and body forces, the equilibrium equations for the system considered above read where n i is the i:th component of the outward unit normal on /7, and p =p,e~ -p , e 2 represents the contact traction vector acting on the support. Two work conjugate measures belong to p. These are u, = -u 2 and w, = u~ -s ~. That is, u n is the normal displacement and w, is the relative tangential displacement between the contact surface and the support. The balance of energy of the body and the contact surface are given by /~ = p~/ + p ,~b , -qzon ~, where p is the mass density of the bulk material, e is the specific energy of the body, E is the density of the internal energy of the contact surface and a superimposed dot represents a time derivative. Defining the absolute temperature T, the specific entropy of the body s and the entropy density of the contact surface S, the second principle of thermodynamics is given by the following Clausius-Duhem inequalities: ~t <~p, ti, + p,ffe -SJ" on ~, where gt = e -sT and q* = E -ST are the Helmholtz free energies of the body and the contact interface, respectively.

Constitutive assumptions--body
Established constitutive assumptions of isotropic thermoelasticity are utilized for the body. 2 That is, defining Lame's elasticity coefficients .~ and /2, the specific heat capacity c, the thermal dilatation coefficient a and a reference temperature T n, the following free energy is considered for the body: from which the following state laws are derived per definition: z The material given in this section is classical and can be found in textbooks (see e.g. [13,14]).
In agreement with this latter inequality, representing the second law of thermodynamics, Fourier's law of heat diffusion is adopted. That is, where k, is the thermal conductivity. Inserting (2.13) in the balance of energy represented by (11) yields This equation is rewritten, by neglecting ~kk and putting T = T o, to read k, Tx, =ocT, (15) which is the standard heat conduction equation of uncoupled thermoelasticity.

Constitutive assumptions--contact surface
The constitutive assumptions of the contact surface are given by Signorini's contact conditions, Coulomb's law of friction and Archard's law of wear. These laws were derived from a generalized standard material, by defining a certain internal state variable o~ measuring the wear gap, in [20]. Utilizing this internal state variable approach, the free energy is assumed to be given by where I c is the indicator function of the closed convex set C (i.e. I c = 0 if (u,,w) E C, 1 c = +oc if (u,,, o~) ~ C) and g is the initial gap between the body and the support. In addition, the following state laws are defined and derived: where °/4/" is the force associated to o~, called the wear driving force. O~ stands for the subdifferential of ~ and defines the normal cone of C. The contact law given in (17) represents an extension of Signorini's contact conditions that includes wear.
To evaluate the implication of the choice of free energy and state laws given in (16)- (18) in conjunction with the balance of energy in (5) and the Clausius-Duhem inequality in (7), we consider evolutions such that the left-hand time derivative and the right-hand time derivative of the state variables and the state functions exist with possibly different values, meaning that the time derivative of these functions and variables may not necessarily exist. That is, the state functions and the state variables are continuous in time, but not necessarily smooth. Thus, a superimposed dot stands for either the left-hand time derivative or the right-hand time derivative.
These two relationships represent the balance of energy and the second principle of thermodynamics for the contact interface. Furthermore, for later use, we remark that the contact conditions in (19) are optimality conditions to the following variational inequality: Next, Coulomb's law of friction and Archard's law of wear are introduced. Coulomb's law of friction is given by the following principle of maximal dissipation: where/z is a constant friction coefficient, I1 stands for the absolute value and (p.)+ = (p. + Ip.[)/2. Assuming that k is a constant wear coefficient, Archard's law of wear is represented by ;, = k(p.)+ Iw,I, (24) One might remark that the use of (p.)+ in (23) and (24), instead of just p., implies that these two representations of Coulomb's and Archard's laws are defined also for p. < 0 and not only for p,,/> 0.
Let us evaluate the consequences of Coulomb's law and Archard's law with respect to the principles of thermodynamics given in (20) and (21). Taking p~ = 0 in (23) implies that ptvb, >/0. Using this and (24), it is shown that the second law of thermodynamics represented by (21) is satisfied. Furthermore, (23) implies that p, =/z(p,)+sgn(vb,) when w, ¢ 0. This together with (24) in (20) result in the following energy balance: q2 :-(k(P,,)2+ +/x(P,,)+)lwtl, (25) which represents the friction and wear heat generation at the contact interface. Thus, all heat generated by friction and wear flows into the body.

The system of discrete equations
Given proper initial conditions, a load history t(t) and a displacement history ~(t) on a time interval [0, r], the following initial boundary value problem is stated: f. ov, oo; fr and 3-is a function space of temperatures. In addition, weak forms of the tribological laws given in (22)- (24) are assumed to be valid on F c (see [21]). Eq. (26) is the principle of virtual work, equivalent to the equilibrium equations given in (1)-(3), with the thermoelastic properties in (9) inserted. Finally, Eq. (27) is the balances of energy for the body and the contact surface given in (15) and (25) put together using (13) and q.n = -q2.

Space discretization
The initial boundary value problem stated above is treated by introducing finite element approximations of 7/" and if (see e.g. [7]), and evaluating the integrals over I~ by an appropriate quadrature rule depending on the choice of finite elements. Let f be a function over F,, then the integrals are approximated by where I M are the weighting factors and r L is the set of integration points x M on 4-When applying this formula, the contact boundary is divided into contact elements. If these elements coincide with the edges of the displacement elements, then it is possible to choose an integration rule such that the integration points will coincide with the nodal displacement points of the contact surface. For two-dimensional problems, examples of such rules are the trapezoidal quadrature rule for linear elements and Simpson's rule for quadratic elements. A study of different integration rules can be found in [8].
After performing the space discretization outlined above, one obtains the following discrete equilibrium equation from (26): • M k to = ~ (P~)+l vb,~l.
Here, u~ = --NA(xM)d A, toM = to(X M) and gM = g(xM). In [21], the above relationships were derived from weak forms of (22)-(24) by using the integration rule in (28). Finally, we rewrite (31) in the following manner for any r > 0 (see [liD: find P,~ t>0 such that

The algorithm
The initial boundary value problem defined in the beginning of Section 3 can now be treated by solving a sequence of the discrete equations given in (29), (34), (38), (40) and (41) using Pang's algorithm presented below. However, numerically tests have proven that the search direction involved in the algorithm can not be solved accurately for this representation of the equations. Comparisons by analytical solutions given in [3] have shown that the Newton step smoothens out the temperature gradients in a wrong way; in addition, the algorithm converges slowly. This depends probably on the great differences in the numerical values of /~, /~ and (M + (tk+ j -t~)O), leading to a system of ill-condition equations. This problem is overcome by putting together (29) and (38) to form
Thus, letting y = (d, P,,, P,), the following system of equations are solved for each time increment:

Pang's Newton algorithm
In [21], it was shown that, for/~ = 0, the equation system in (43) is Bouligand differentiable, i.e. Lipschitz continuous and directionally differentiable. Obviously, this holds also for the case with /~ # 0 as/~(d, P,,) is Lipschitz continuous and directionally differentiable. For solving systems of Bouligand differentiable equations, Pang [16] suggested the following globally convergent algorithm (concerning convergence result of the algorithm, see also [4]): ALGORITHM. Let fl, y and • be given scalars with fl E (0, 1), y E (0, 1/2) and e small. Repeat the following steps for each time increment k + l" (0) Let yO be given from the previous time step k and set l = 0.
(1) Find a search direction z = (za, z,,,zt) such that 8 t I where H (y ;z) is the directional derivative. (2) Let a ~= tim,, where m~ is the first non-negative integer m for which the following decrease criterion holds:
(4) If g(f+l) ~< e, then terminate with y/+ 1 as an approximate zero of H(y). Otherwise, replace l by 1 + 1 and return to step 1. After convergence in each time step is achieved, the variables T, w~ and w M are updated.

Approximation of directional derivatives
The time-consuming part of the algorithm is to find the search direction in (46). Generally, this system of equations is non-linear due to the non-linearity of H'(y; z) in the second argument, which in turn depends on the non-differentiability of H(y). Of course, H'(y;z)=VH(y)z at states where H(y) is differentiable. The nonlinearity of (46) might limit the effectiveness of the algorithm. This drawback is overcome by introducing an appropriate approximation of H'(y;z) such that (46) results in a system of linear equations (see [4,21,22]). The

Numerical examples
Five different problems are considered for a structure subjected to boundary conditions and a displacement history as given in Fig. 2. The problems differ only in the constitutive settings, i.e. the same boundary conditions and displacement history (except for the number of cycles) are used for all five problems. The dimensions of the structure are 3.0 X 1.0 cm 2. Furthermore, the structure is fixed at the top and unilaterally constrained at the bottom by a rigid support using a negative initial gap of g = -0.1 txm. In addition, the support is subjected to a cyclic tangential displacement history with an amplitude of 1.0 mm according to Fig. 2. In the context of wear, this type of boundary condition can be called a fretting condition. Although ~: is small, this condition will result in global sliding fretting. It is a major principle difference between the contact stresses obtained for global sliding fretting compared to stick-slip fretting. This is discussed in more detail further on. The thermoelastic properties of the structure are approximated by using 42 X 28 bilinear quadrilateral elements. 3 The potential contact surface at the bottom of the structure is approximated by 42 contact elements using the trapezoidal rule in such manner that the integration points coincide with the nodal displacement points of the contact surface. Concerning the time discretization, each cycle is treated using four time steps, i.e. tk+ j = t k + 0.01 s.
For the system given above, the following problems, defined by five different constitutive settings, are considered: the first problem is pure elastic, meaning that no thermal properties are used, with /.t = k = 0, the second problem is pure elastic with /~ = 0.3 and k = 0, the third problem is pure elastic with /x = 0.3 and k= 1.0x 10 -tj Pa -~ (this is a typical wear coefficient for adhesive wear in steel assemblies), the fourth problem is thermoelastic with # = 0.3 and k = 0, and, finally, the fifth problem is thermoelastic with /z = 0.3 and k = 1.0 X 10-~ Pa -j.
The contact pressures of problem one (/1, = k = 0, no heat) and problem two (/x = 0.3, k = 0, no heat) are compared in Fig. 3. The only difference between these two problems is the value of the coefficient of friction. Obviously, p,, is constant in time for/z = 0, but, for # = 0.3, p,, depends on the sliding direction of the support. In Fig. 3, p,, is plotted at states when t = n X 0.04 s, where n is a positive integer. The solution will be the same at times when ~ > 0 and global sliding is developed. On the other hand, the reflected image of p,, will appear when there is global sliding with ~ < 0. Furthermore, the dissipation is non-zero for/z = 0.3 in contrast to the case for /z = 0 when it is zero. Actually, the frictional dissipation is constant for each cycle. Thus, if one is interested in studying the dissipation for a large number of cycles such as the dissipation in frictional dampers, then the constitutive assumptions of problem two might not be sufficient due to the fact that real engineering materials eventually will be broken during a repeated frictional process and consequently the dissipation will eventually be equal to zero. This fact is included in the constitutive settings of problem three and five. For those  x I tml Fig. 3. Comparison of the contact pressures for problem one (tz =k=0, no heat) and problem two (/z =0.3, k=O, no heat).
problems the frictional dissipation turns toward zero after many cycles depending on the surface break down due to wear.
In Fig. 4, the contact pressure and the wear gap are plotted for problem three (/z = 0.3, k = 1.0 X 10-j~ Pa -~, no heat) after 30 cycles (t = 1.2 s). At this instant, the contact pressure is approximately zero since the initial gap is almost worn away, i.e. w ~ -g = 0.1 ;xm over the contact surface. However, during the first cycle of this fretting process, the wear gap is approximately zero and the contact pressure is similar to the pressure obtained in problem two, see Fig. 3. Thus, although a physically realistic value on the wear coefficient is used, a large reduction of the contact pressure has been obtained within a few number of cycles. This depends not only on the global sliding conditions but also on the fact that the wear model used in this work does not include the transportation of wear debris in the contact interface but rather annihilates the wear debris. Of course, the rate of reduction in contact pressure can be decreased by diminishing the wear coefficient. Perhaps a better approach to obtain such a reduction would be to extend the formulation used in this work by a third body of wear debris, following the ideas of Godet [6]. Another approach in that direction is accomplished in problem five, where the rate of reduction in contact pressure is decreased by adding thermal properties to the constitutive settings of this problem. One might notice that the result obtained above would be completely different if the fretting conditions would be chosen such that no global sliding would be developed. If the amplitude on ~: would be taken sufficiently small such that no slip would be developed over the contact surface (global stick), then, of course, no wear would be initiated. However, in the intermediate with amplitudes on ( implying stick-slip, the contact pressure would turn toward zero in the slipping regions and the pressure would increase in the sticking regions (see [21,22]). In addition, large values and large gradients in the contact stresses would be developed in the regions between stick and slip. Concerning the field of fretting fatigue, such stress distributions might initiate cracks near the contact surface leading to fatigue failures of the structure. Similar observations have been found experimentally by using so-called fretting maps in e.g. [17,23].

11
In Fig. 5, two frictional hot spots and two peaks in the contact pressure are plotted after 30 cycles for problem four (/z = 0.3, k = 0, heat). After the first cycle the contact pressure is similar in shape to the pressure of problem two, see Fig. 3, but due to the thermal expansion induced by frictional heat generation the contact pressure changes rapidly from this shape during the evolution of cycles. The shape obtained in the temperature field and the contact pressure are very sensitive to the geometry and the boundary conditions. For instance, if the dimensions of the structure are changed to 1.0 × 1.0 cm 2, then only one frictional hot spot will appear. Furthermore, if the velocity of the support is decreased sufficiently, then no hot spots at all will appear. An obvious drawback of this constitutive model, concerning this type of evolution problems, is that no damages are induced due to the high temperatures and high contact stresses, neither in the bulk material nor in the contact surface. In real engineering materials such high temperatures and high stresses would result in damages of the material in some way, which in turn would reduce the amount of frictional heat generation and consequently result in lower temperatures and pressures than predicted here. This fact is included in the constitutive settings of problem five by adding Archard's law of wear to the constitutive laws used here. In that problem, the surface damages in form of wear imply that no frictional hot spots are developed.
In Fig. 6, the temperature field, the contact pressure and the wear gap are plotted after 210 cycles (t = 8.4 s) for problem five (/x = 0.3, k = 1.0 × 10 ~ Pa -~, heat). At this instant, the solution has almost reached a steady state. The temperature field is approximately constant at 0.0987 K. Notice, in comparison to problem four, that no frictional hot spots are developed. Instead the increase in temperature is astonishingly low. This fact depends on the wear process, which decreases the frictional and wear heat generation. However, following the discussion on problem three, this reduction is probably too large due to the annihilation of wear debris. If the wear debris would be included in the model as a third body, then the frictional and wear heat generation would be higher, but not as high as in problem four. Although the temperature is very low, it has a major effect on the wear rate. After seven times more cycles compared to problem three, the contact pressure is still much higher here. This depends on the thermal expansion. Furthermore, the thermal expansion also increases the amount of wear; the wear gap is no longer constant and the maximum in wear gap obtained is more than 20% larger in comparison to problem three.   The temperature rise in fretting is a subject of considerable interest. Some researchers claim that the presence of et-Fe20 3 in the wear debris must be taken as evidence that the interface temperature exceeds 300-600 K in fretting (see e.g. [24]). Obviously, such high temperatures are not predicted in problem five. Nevertheless, by neglecting the removal of material, as in problem four, even higher temperatures are predicted, see Fig. 5. However, measurements by other researchers indicate that the average temperature in contacting asperities subjected to fretting conditions is much lower; magnitudes of 20-30 K have been measured (see e.g. [5,19]). Such temperature increases are more comparable to the increases obtained by the wear model used in the present paper.
For fretting problems with stick-slip conditions, the frictional heat generation has a minor influence on the contact pressure and the wear gap using this wear model. Still, other initial and boundary conditions on the temperature and the heat flux vector can have a major effect on the contact pressure and the wear gap for stick-slip fretting problems. For instance, one might think of problems where external thermal boundary conditions induce stick-slip fretting. This type of problem can of course be treated with the algorithm suggested in this work.
Finally, some execution statistics for the five problems considered are provided. The number of Newton iterations per time increment (New./inc.) and the number of line-searches, i.e. Armijo iterations, per Newton iterations (Arm./New.) are given in Table 1. In order to get a fair comparison, the statistics are determined for Table 1 Some execution statistics for the five problems considered Problem-constitutive setting New./inc. Arm./New. 3.07 the ten first cycles of each problem. Apparently, problem four is the most time-consuming to solve, and problem one and two are the fastest ones. The reason for problem two to have the same statistics as problem one is that there is only global sliding and no stick-slip developed in the solution of this problem as well as in the solutions of the other problems. If there would be stick-slip in the solutions, then the statistics would be higher for problem 2-5, at least the number of New./inc. would be higher.

Concluding remarks
In this work thermoelastic wear problems have been solved by a finite element method. The method is a direct method, meaning that no iterative strategy between a mechanical and a thermal part of the problem is used, but instead the finite element equations are solved using a Bouligand differentiable Newton method. The method obtained is effective and robust.
Numerically, it has been found that the internal state variable measuring the wear gap has a great influence on the temperature field generated by frictional heat. For a particular frictional heat generation problem, it is numerically shown that when this internal state variable is included in the analysis, then the increase in temperature is very low. On the contrary, if wear is excluded, then two frictional hot spots are developed for the same problem. However, although the temperature increase is very low for the wear problem, the frictional heat generation has a major effect on the wear rate.
The wear model utilized in this work does not include the transportation of wear debris, but instead the internal state variable measuring the wear gap rather annihilates the wear debris. From the numerical result discussed above, it is clear that the transportation of wear debris will affect the frictional heat generation which in turn will affect the wear rate. Therefore, for future plans, it would be of interest to extend this wear model by a third body of wear debris and to study its influence on the wear rate.