Contact constraints within coupled thermomechanical analysis-A finite element model

Based on a sophisticated interface model a thermomechanical contact formulation is derived. The development addresses theoretical and numerical aspects. It starts from a continuous formulation of the local contact geometry and includes constitutive equations for the contact stresses, the contact heat flux and the frictional dissipation. Based on these considerations, a finite element model for large deformation processes is developed. An operator split technique concerning the solution of the global coupled equations is applied for the algorithmic treatment of the thermomechanical coupling. By means of three numerical examples, the theoretical and numerical formulations are validated.


Introduction
In engineering problems such as shrink fitting, deep drawing or hot rolling, the effects of heat generation and heat transfer have to be considered during the process.Since many of these simulations involve contact problems with dry friction, it is of interest to develop contact laws that are also able to predict, besides normal and tangential stresses, the amount of heat generated due to friction and the heat transfer across the contact surface.
During recent years, a series of publication was devoted to the derivation of thermomechanical interface models that take into account the microgeometrical shape of macroscopically smooth surfaces (see e.g.[l]).The associated theoretical studies consider a statistical characterization of the surface profile and of the related thermomechanical phenomena and usually pre-define the shape and the mechanical behaviour of the asperities.The derived constitutive relations for the thermomechanical contact interface also include effects due to the finishing treatments of the surface layer which involve different characteristics from the bulk material.Hence, the effective yield pressure (determined by hardness tests) depends on the depth of the plastic zone [2].
Connected to these sophisticated interface models for the heat transfer within the contact area are constitutive relations for the mechanical contact which are also based on statistical models (see e.g.[3] for an overview).Here we have to distinguish between constitutive equations for the normal contact and the tangential contact.For the normal contact, statistical models have been stated in [4].For the tangential contact, which splits into the stick and slip part, several constitutive relations have been developed (see e.g.141).
Numerical formulations of contact problems within the finite element method have been given in many papers.The standard approach is to model the contact surface as an idealized smooth interface which then leads to a condition describing the non-penetration by a pure geometrical constraint.Within this assumption, several approaches are possible for the variational formulation of contact problems in case of finite deformations.
The so-called Lagrangian multiplier technique enforces the contact constraints exactly, (see e.g.[S, 61).A consistent formulation for implicit computations using penalty methods is given in e.g.[7] based on a discrete model discussed in (81.Recently also constitutive equations within the contact interface have been considered in finite element approximations (see e.g.]91)* Finite element formulations for the frictional contact have been discussed e.g. in 13, g-121.Finite element models which include thermomechanical coupling for contact problems have been presented in e.g.[13] for lubricated contact, in [14] for large deformation frictionless contact and for small deformations in 1151 including frictional heating.
In this paper we focus on a finite element formulation of a sophisticated interface model that includes frictional contact, heat exchange in the contact area and frictional heating.For this purpose, we give (i) a continuous formulation of the local contact geometry based on a parametrization of a moving master surface, where we restrict our consideration to two-dimensional problems.This formulation extends the formulation proposed in fl2] for the rigid obstacle problem.Then we postulate (ii) constitutive equations for the contact stress, the contact heat flux and the frictional dissipation on the contact surface, These equations have model character but exhibit essential features of known experimental and microgeometrical investigations of contact interfaces.After having formulated (iii) the structure of the initial boundary value problem of finite strain thermoelasticity combined with thermomechanical contact, we discuss (iv) two main algorithmic aspects.Firstly we discuss the solution algorithm for the coupled global equations.Global solution algorithms for coupled thermomechanical analysis have been formulated in [16-B].Here we use an algorithm that is based on an operator split of the global thermomechanical evolution operator.This strategy was recently introduced in the context of finite strain thermoplasticity without contact constraints in [19,20] and leads to an algorithmic decoupling of the coupled problem within the time step by preserving the symmetry properties of the subproblems.Secondly we consider an extension of the return map algorithm proposed in [21] for the update of the contact stress and frictional dissipation.Finally we discuss (v) the spatial discretization of the contact surfaces and show as an interesting aspect the transition of the continuous contact geometry to the discrete contact geometry of the three node contact element used e.g. in [9] for frictional contact.
Several numerical examples show the performance and the accuracy of the derived finite numerical formulation.

Contact geometry
Let $8 (I, a = 1, 2, denote two bodies of interest and (0: : % a -+ R" the associated deformation maps at time tfR+.4~: maps points X" E W" of the reference configuration onto points _ra = y)y(Xa) of the current configuration.
Assume that both bodies come into contact.Motivated by micromechanical investigations of contact problems, we view mechanical contact as a (microscopical) penetration of the current mathematical boundaries qp(rz) where rz C ZBO are possible contact surfaces of the bodies B3" ; see Fig. 1 for an illustration of this idea.In what follows we denote q_!(rr) as the current slave surface which penetrates in the case of contact into the current master surface qf(fz).
is the orthogonal projection of a given slave point x1 onto the n := (a; x&/II -2 a, x tiill is the outward unit normal on the current master surface at the master point where Zz are tangent vectors at $(g', 5").Here and in the following we denote by a bar over a quantity its evaluation at the minimal distance point (c', k").
The penetration function gives two pieces of informations: Firstly it serves as a local contact check, i.e. we set contact @ggN > 0. Secondly it enters for g, > 0 as a local kinematical variable the constitutive function for the contact pressure (see Section 3.1).
By taking the time derivative of Eq. (2.2) at the minimal distance point (g', i'), one obtains, in the case of contact, the rate of penetration -& = [U: -$(<', C")] -ii* for given spatial velocities V: and Lj:(cl, 5") at the slave and master points.

Tangential relative velocity and tangential relative slip
The tangential relative slip between two bodies is related to the change of the solution point (z', i') of the minimal distance problem.Thus, we can compute the time derivative of 5" from Eq. (2.3).This yields the following result: with with cw llj7; = $*(5', 57 Jr $&*, I%",.we obtain ip from the following system of equations: Note that Ix: -if{<', E")] has been written with Eqs.(2.1) and (2.2) as --gNi*.A-well known result from differential geometry introduces for $,,( c', $") * ii2 the curvature tensor b,,.Thus, we can rewrite I?,, = [ct,, -g,G*J.We define as the second important kinematica function (with a view to problems with friction) a tangential relative velocity function on the current slave surface e:(ri) by setting &g, := i$"a, . (2.8) Eq. (2.8) determines per definition the evolution of the tangential slip g, which enters as a local kinematical variable the constitutive function for the contact tangential stress (see Section 3.2).The rate i" in Eq. (2.8) at the solution point (r', 8') has already been computed in Eq. (2.6), REMARK 1.Note that the second term on the right-hand side of Eq. (2.7) depends on the penetration g,.Thus, in the case of a strong enforcement of the non-penetration condition (g, = 0) with Lagrangian multipliers, this term vanishes.Then the evolution 2" in Eq. (2.8) is given by the projection of the spatial velocities vi and 6:(i) evaluated at the slave and master points onto the tangential direction of the master surface at the master point.For a penetration g, > 0, we have to take into account the second term in Eq. (2.8) and the scaling factors Hap, both consequences of the time dependence of ii.

Frame invariance
Consider the vectors g, : = -g,ti2 and g, := ia& normal and tangential to the master surface at the master point.Clearly these vectors transform objectively under rigid body motions superposed onto the spatial configuration (xf)' = Q(t)x: + c(t) VQ(t) E SO(3) w h Example: (gN)+ = QgN is a consequence of (6')' = Qi'.
ere SO(3) is the proper orthogonal group.
Observe furthermore that the (a priori objective) Lie derivative of the tangential vector g, has the representation based on the deformation gradient Ff of the master surface defined above.Thus, Eq. (2.3) represents an evolution equation for the objective rate JZUgT of the tangential vector introduced above.

Constitutive modelling on the contact surface
When the micromechanical behaviour of the contact area is studied, different phenomena have to be considered for the mechanical as well as for the thermal interface description.In this paper, we use constitutive models which have been derived based on micromechanical observations of physical contact surfaces.These models can be related to formulations relative to mathematical contact surfaces by an averaging process as symbolically indicated in Fig. 3.The goal of this section is to formulate local constitutive equations for the pressure, tangential stress and normal (outward) heat flux on the slave surface at the slave point x1 relative to the bases {CT, &, ri'} acting on body 9? ' ; see general structure in Eq. (4.7).

Contact normal stress
It is well known that the contact pressure is related to the approach of the physical surfaces which come into contact, i.e. the penetration of the mathematical surfaces results from the deformation of the micro-asperities (see Fig. 3).Based on experimental investigations Kragelsky et al. [4] formulated the following non-linear elastic constitutive equation for the contact pressure: in terms of the penetration g, defined in Eq. (2.1).Here cN and m are material parameters which have to be determined by experiment.

Contact tangential stress and plastic tangential slip
In this paper we restrict our consideration to classical Coulomb-type friction problems.Coulomb's constitutive equation can formally be formulated in the framework of elastoplasticity.This has been considered by several authors (e.g.[22,23]) and within a finite element formulation by Wriggers [21].
The key idea of this approach is a split of the tangential slip component g, into an elastic part gt (microdisplacement describing the stick behaviour) and a plastic part g! (frictional slip), (see Eq. (3.2b)).We assume a linear elastic constitutive equation for the tangential contact stress component To obtain the relation representing the heat transfer associated with the spots, the behaviour of one spot is analyzed in detail {24].Here, a heat flux tube having a narrowing in the contact zone is considered around each spot.Then the mean value of the heat transfer through the spots is obtained from a statistical model.The additional effect of taking into account the hardness variation with the mean planes approach has been suggested in [2] leading to the following interface conductivity for the spots: (3.9) Here c 1, c2 describe the harness variation and k is the mean thermal conductivity depending on the conductivities of the two bodies contact.E is the mean absolute asperity slope and (r denotes the surface roughness.All values have to be deduced from experimental data.
The heat transfer due to gas or liquid within the microcavities is mainly governed by conduction.This fact results from the small height of the cavities which do not allow for convective flow.Based on this assumption, Yovanovich [25] derived a relationship, which also takes into account the change of cavity height by the contact pressure pN.With the current gas temperature 8,, a constitutive constant crc for the gas and the Vickers hardness H,, the gas conductivity k, and the surface roughness rr, the heat transfer through the gas yields

Initial boundary value problem
In the previous sections we have discussed the contact geometry and the constitutive equations associated with contact interface.Let us now formulate the initial boundary value problem for non-linear thermoelasticity combined with thermomechanical (frictional) contact.In order to obtain a compact structure, we introduce locally at the material point X" E sA" and time t E 118, a vector of primary variables F(X", t) := {p*(Xa, t), V"(X", t), @(Xa, t)} (which contains the configuration Q~, the material velocity V* and the material temperature Oa) and locally at points X" E rz of the contact surfaces a vector of history variables sl(Xff, t) := (g$(Xa, t)} (which contains the plastic tangential slip g!).Internal variables in the domains W" do not appear since we restrict our presentation to thermoelastic constitutive response.The thermomechanical initial boundary value problem is governed by the local field equations in the domains BB",

Solution algorithms
The main solution strategies for coupled problems are monolithic schemes where the differential equations for the different variables are all solved together, or staggered schemes where the different variables are solved separately (see [26] for an overview).Global solution algorithms for coupled thermomechanical analysis have been formulated by Argyris and Doltsinis [16], Miehe [17], Doltsinis [18], Simo and Miehe [19], among others.In what follows, we adopt the staggered scheme proposed in [19] based on an operator split of the global thermomechanical evolution operator discussed in Section 4. This strategy yields an algorithmic decoupling of the thermomechani~al equations within a time step on the basis of (in the frictionless case) symmetric subproblems.
A further objective of this section is to formulate the local update algorithm for the contact tangential stress and the frictional dissipation.Here we extend the return mapping algorithm proposed e.g. in [21] for isothermal frictional contact to the non-isothermal case.
The central idea is a split of the evolution operator A" in Eq. (4.2) into its natural mechanical and thermal parts as indicated in Eq.Thus the operator split algorithm results in an algorithmic decoupling of the coupled thermomechanical equations within the time interval.The proposed split (5.1) yields an identical coupling algorithm to that suggested in [16], characterized by an isothermal deformation predictor followed by a heat conduction corrector.Within their work, this algo~thm has been interpreted as a one-pass Gauss-Seidel scheme in terms of the variables +Y and 0 (see also [18]).
In view of the finite element treatment, we construct the weak forms of the time discrete field equations (5.2a) and (5.3) by standard Galerkin procedures.The weak form of the mechanical subproblem ALGO, takes the form The algorithmic update of the tangential force is performed by the return algorithm as described in Section 5.2.Note that we have only one component in the two-dimensional case.Analogous to Eqs. tate (5.6)-(5.9),we compute the elastic trial s perform the return algorithm Eqs.(6.1) and (6.2) characterize the main kinematical relations of the contact element in Fig. 4.
In what follows, we formulate the constitutive relationships for the normal force PN, tangential force TT, the discrete heat flux QN and the discrete dissipation gap at the discrete contact point(s) of the contact element under consideration.Assume that &+ I and gNn+i are known from the exploitation of Eq. (6.1).Then analogous to Eq. (3.1), we compute the normal force Thus the virtual mechanical work, Eq. (6.8a), of the contact element can be written in the matrix formulation Gk = it .Rkn+, with the mechanical contact element residual, In Eq. (6.16), the contact normal force PNn+l follows from Eq. (6.3) whereas the tangential force TTntl is given by the return algorithm (6.5).Due to this approach, a pure displacement formulation of the contact problem is possible which is in contrast to the Lagrangian multiplier technique often used to enforce Eq. (2.1).For a global algorithmic treatment using Newton's method, we have to linearize Eq. (6.16).The associated formulation for this discretization can be found in [9].
The matrix formulation of the thermal part, Eq. (6.8b), is similar to the mechanical part.As a consequence of the global operator split algorithm discussed in Section 5 and the assumed simplified constitutive equation, Eq. (6.7) is the thermal part linear.
A formulation completely analogous to Eqs. (6.11)-(6.15)yields the matrix representation GG = & + RLI for the virtual thermal work Eq.(6.8b) with the thermal contact element residual Rk,cl:= QNn+lHn+I -gb+lHon+l* (6.17) Here we have introduced the matrices (6.18) analogous to Eq. (6.13) and (6.19) for the matrix formulation of the thermal variations (6.10).In Eq. (6.19) the discrete contact heat flux & Nnil follows from Eq. (6.7) whereas the frictional dissipation BL+l is given by the return algorithm (6.6).The nodal contact contribution Eqs.(6.16) and (6.19) have to be added to the global system of equations. 13

Numerical examples
The proposed contact-model has been implemented in the finite element code FEAP and was tested by means of several examples.For the following examples, we use finite element formulations for large strain coupled thermo-elastic and thermo-plastic response, which are documented in detail in [17] and [ 191, respectively.The first example is a test for the developed and implemented pressure dependent thermomechanical constitutive equation.The system and its finite element discretization is depicted in Fig. 5.A temperature of 100 K and a pressure p is applied to the upper surface of 24 2 whereas the bottom of L% ' is fixed and has a temperature of 0 K, the vertical boundaries are thermally isolated.Body 3' is assumed to be rigid but can conduct heat.Body 33' is the~oelastic.
The thermoelastic ~onstitutive parameters for aluminium are given in Table 1.These results clearly show the pressure dependency of the constitutive law for the heat transfer in the contact area.They depict also that the finite element approximation is in very good accordance with the analytical solution.

Frictional heating of a block on a rigid surface
The mechanism of heating due to friction is studied in this example.For this purpose, we consider an elastic block 9 * which is slid over a rigid block 93 '.Both bodies are heat conductors.The upper body moves within 3.75 X 10d3 s from the left to the right end of the lower block, the displacements are prescribed at the top of ~4 2. During this process a pressure of p = 10 N/mm* is applied on block 93 *.The system and its finite element discretization are shown in Fig. 6.In order to check the numerical result, we consider idealized thermally isolated surfaces of both bodies.Again we assume that both bodies consist of aluminium with the material parameters given in Table 1.The frictional coefficient is assumed to be p = 0.2.
A total of 100 time steps have been used to compute the evolution of the temperature during the process.The temperature distributions at time t = 0.75 x 10m3 s and t = 2.25 x 1O-3 s are shown in Fig. 7.After that time, the movement of the upper body was stopped.
Another 50 time steps with increasing step length were applied to compute the temperature evolution up to the final homogeneous steady state at the process time of 10 s.After that time, the temperature is balanced and thus equal at all elements of the bodies %a.Since in this example, we wanted to evaluate

Upsetting of a billet
A thermoplastic upsetting process of an aluminium block is considered under plane strain conditions in this example.The block is pressed between two rigid plates which are able to conduct heat.Within this process, heating of the block occurs due to plastic dissipation within the body and frictional dissipation within the contact surface.The material parameters are listed in Table 3.We have used a temperature-dependent yield stress, given by the function j(eP, 0) = [l -w,(O -O,,)]yO + [l -~~(0 -OO)]heP.The finite element mesh is shown in Fig. 8. Due to symmetry, only one quarter of the system was discretized by 340 finite elements with a total number of 591 unknowns.
The billet is deformed within T = 0.0035 s.A total of 100 time steps were needed to arrive at the final reduction of 42%.Fig. 9 shows the temperature distribution and the deformed system at the final --r---nIsPLfxzmwr 3 j.9'35Ef,3@ to 8.028E+00 reduction.The maximum temperature increase in the middle of the specimen is 28.2 K whereas the average temperature increase is 16 K within the contact interface.
Fig. I. Basic contact geometry.

REMARK 2 .
If we have a flat contact surface the curvature tensor baP is zero.

(
Fig. 3. Averaging of micromechanical contact relations. t,=c,gt with gt:=g,--g;, (3.2a,b)where cT is a material parameter.The plastic tangential slip gk is governed by a constitutive evolution equation which can be formally derived by using standard concepts of elastoplasticity theory.Let BP:= tTxYY"g;~O (3.3) be the dissipation due to the plastic slip.Now consider an elastic domain IE, : = {tT E R* 1 fp(tl.)s 0} in the space of the contact tangential stress.Here is the plastic slip criterion function for given contact pressure pN with materig parameter /A.Holding pN fixed, one obtains from the so-called maximum dissipation principle, (tT -tT) * Tug! 2 0 Vt, E E,, the constitutive evolution equation for the plastic slip --T with nT -,,~~,, (3.5) (normality rule for fixed contact pressure) along with the loading-unloading conditions in Kuhn-Tucker form h20, .(p(tT)so,Ajp(tT) = O, (3.6) which determine the plastic parameter A.3.3.Contact heat fluxDifferent techniques for the computation of the thermal contact resistance, taking into account the dependence ,on various parameters have been proposed (see e.g.[l]).We assume the following structure for the constitutive equation for the heat flux qN (outflux through slave surface):qN = G,(O', o*, pN) = h(O', o*, PN)(O1 -0').(3.7)The contact heat transfer coefficient h = h(O', O*, pN) depends on the current temperatures 0' and O* at the slave and master points, x1 and i:(s), respectively, as well as on the contact pressure pN defined in Eq. (3.1).The resistance against heat transfer is mainly due to the low percentage of physical surface area which is really in contact.The presence of a reduced set of spots surrounded by microcavities characterizes the contact zone, hence heat exchange is possible by heat conduction through the spots (S), heat conduction through the gas (G) contained in the cavities and radiation (R) between microca~ty surfaces.These mechanisms act in parallel which induces the additive decomposition (3.8) of the contact heat transfer function.Based on statistical micromechanical models, different relations for the gas and spot conductance have been developed (see e.g.121).
(3.10)In this paper, we use a simplified relation which neglects radiation effects and gas conductance but takes into account pressure dependency of the contact conductance ' and h,=h,=O, (3.11) based on the material parameters h,,, He, e.This model is valid for high pressures.A more sophisticated constitutiv~ model is implemented in f14] for the frictionless case.
, C. Miehe / Comput.Methods Appl.Mech.Engrg.113 (1994) 301-319 $d =V"+O, de~nition of the material velocity, the balance of linear momentum and the balance of internal energy in form of the temperature evolution equation.Pa is the first Piola-Kirchhoff stress tensor, Q" denotes the material heat flux vector and S" is a heat source which describes the Gough-Joule coupling effect in the framework of thermoelasticity.The field equations (4.1) form a coupled first order evolution system for the primary variable vector Y* introduced above, i.e. (4.2) where the non-linear evolution operator d represents the right-hand side of Eq. (4.1) by taking into account the dependency of the thermoelastic constitutive equations in the domains B3" ".Next the initial conditions for the primary and history variables for a time interval [t, , t, + ,] E R + of interest are given by Y""(XP, t,) = Y","(Xll) and .Y"(X", t,) = 3:(X*).(4.4)We consider boundary conditions for the deformation and the temperature on rz C &!BiP" and fg C FIB* va(Xa!, r) = 6*(X", t) on rs , Oa(Xu, t) = 0*(X", t) on r; .(4.5) as well as for the traction vector and the heat flux on I'; C ZZ?Y and I'% C &?V tn(4bzYXa), t) = i*(V;(X"), t) on qpP(rta) Y 4%2X?, @la(X% t> = S"(&(X*>, Op(xa>, t) on icop .

Fig. 6 .Fig. 7 .
Fig. 6.System of frictional heating problem ] = 9.375 N mm.From this quantity we obtain the average temperature of the isolated system S? A* for an infinite process time Sz = Wdi,,/(2pacaVa) which yields the results shown in Table2.Again we observe good agreement of the global energy check with the finite element solution based on the proposed model.

Table 1
Material parameters for examples in Sections 7.1 and 7.2

Table 3
Material parameters for the example in Section 7.3