Finite element analysis of transient thermomechanical rolling contact using an efficient arbitrary Lagrangian–Eulerian description

A theoretical and computational framework for the analysis of thermomechanically coupled transient rolling contact, based on an arbitrary Lagrangian–Eulerian (ALE) kinematical description, is developed. A finite element formulation featuring 2D cylinder–plate rolling contact is implemented. The implementation features penalty-type contact formulations for mechanical and thermal contact. It is noted that the ALE formulation allows for a simplified time description, a compact computational domain and localized mesh refinement. Numerical simulations considering stationary and transient rolling conditions are presented. Highlighted aspects include the influence of variations in thermal contact conductivity, rolling speed and external mechanical load on the contact interface heat flow. The model is shown to give predictions in qualitative agreement with results in the literature. For the velocity range studied, numerical issues such as spurious numerical dissipation/oscillations in the temperature field are noted to have a prominent influence. These phenomena are addressed using a Streamline-Upwind Petrov–Galerkin stabilization scheme together with a bubble function approach.

expansion, dissipative heat generation and frictional heat generation [1] are all modes of thermomechanical coupling which may significantly influence the mechanical and thermal behaviour of the contacting bodies, as well as their material properties. Examples of applications where a thermomechanical analysis might be crucial for an accurate description of the response include rolling mills, roller bearings and wheel-rail contacts [2]. It should here be noted that the flow of heat through the contact interface (in railway applications denoted "rail chill"), as well as the partitioning of the frictional heat generated at the contact interface, depend strongly on the rolling velocity [3,4]. This point will be further discussed in subsequent sections.
A common approach to numerical modelling of rolling contact is to employ a semi-analytical contact model based on Hertz theory [5][6][7]. Such an implementation is simple and fast, but limited by assumptions of elastic material response, smooth, continuous surfaces and a small relative dimension of the contact patch. In contrast, a contact formulation involving deformation dependent contact forces and an iterative contact search algorithm [7,8] adds a significant (often much-needed) versatility. However, this approach is relatively complicated to implement and increases computational demands.
In finite element modelling of solid mechanics problems, a Lagrangian kinematical description is usually employed. When modelling rolling contact in particular, a problem formulation based instead on an arbitrary Lagrangian-Eulerian (ALE) description provides important advantages [8,9]. These (elaborated in following sections) include the possibility to linearize the mechanical response, the compactness of the computational domain, the simplified time description, the possibility for highly localized mesh refinement and the elimination of the need for velocity-dependent contact conductivity modelling.
The mathematical foundations of mechanical rolling contact in the context of an ALE description were first presented in Nackenhorst [9]. Here, the relevant kinematical description, balance laws, weak forms and contact kinematics are discussed in elaborate detail. The paper highlights the advantages of the ALE approach, but also discusses complications stemming from the intrinsic difficulty of tracing material points in this case. These include the difficulty of handling inelastic material behaviour as well as keeping track of relative slip distances. Both of these issues were addressed in Ziefle and Nackenhorst [10], where a staggered solution scheme is suggested, i.a. involving solution of the advection equation to keep track of internal variable data. In a recent paper [11], the use of an ALE method for thermomechanically coupled stationary rolling contact is introduced. Here, a thermoviscoelastic constitutive model is employed, involving large deformations and temperature-dependent constitutive parameters.
The present paper aims to provide a comprehensive presentation of thermomechanically coupled transient rolling contact in the context of an ALE description. Specific attention is here given to the stress, the deformation-induced heat source, the heat flux and the mechanical and thermal boundary conditions as they appear in the ALE framework. Furthermore, issues of numerical inaccuracy related to the solution of the discretized energy balance equation are discussed in detail, and stabilization measures are suggested.
The paper is outlined as follows: In Sect. 2, the kinematical description of the rolling cylinder on the plate is presented. The thermomechanical boundary value problem, including its finite element formulation, is defined in Sect. 3. Section 4 contains a discussion about element choice and numerical stabilization measures necessary for an accurate numerical response. In Sect. 5, the mechanical and thermal contact formulations are described. Numerical examples are presented in Sect. 6. In particular, the model is validated and the influence of key operational parameters is evaluated. With Sect. 7, the paper is concluded with a summary and an outlook towards future work.

ALE description of rolling contact
The employed ALE description of cylinder-plate rolling motion [12] can be described as follows: Two intermediate configurations are utilized in addition to the initial ("undeformed") configuration, Ω X, and the current ("deformed") configuration, ω x. The  urations feature a moving coordinate system that follows the cylinder centre. The map fromω to ω accounts for a pure translation of the system in going back to the original fixed coordinate system. The maps between the configurations are formulated aŝ X =φ(X, t),x =φ(X, t), x =φ(x, t) and x = ϕ(X, t), so that ϕ(X, t) =φ(φ(φ(X, t), t), t). The corresponding deformation gradients are, in respective order:f ,F,f and F. Figure 1 shows a schematic illustration of the employed configurations, with their intermediary maps and deformation gradients.
The rigid body mapsX =φ(X, t) and x =φ(x, t) can be expressed a priori from knowledge of the translational and rotational motion of the cylinder along the plate: where R is a rotation tensor, X 0 is the position of the cylinder centre in the undeformed configuration,X is the translation of the cylinder centre and Ω c and Ω p are the subsets of Ω representing the cylinder and plate domains, respectively. Note thatf = R andf = I in the domain of the cylinder and plate, respectively, and thatf = I. The Lagrangian problem of finding the map x = ϕ(X, t) (or the displacement u(X, t) = ϕ(X, t) − X) is thus narrowed down to the ALE problem of finding the mapx =φ(X, t) (or the displacementû(X, t) =φ(X, t) −X). For small strains, the ALE displacementsû will be small, which is generally not the case for the standard Lagrangian displacements u. Consequently, in the former case (but not the latter) it is possible to linearize the mechanical response. In particular, it is noted that the total deformation gradient can be expressed as sincef = I. The linearization of the material behaviour pertains to the assumption thatF ≈ I whilef is arbitrary. This point is elaborated in Sect. 3. Another advantage of the presented convective kinematical description is the fact that it allows for a compact computational model: only a relatively short section of the plate domain needs to be modelled, regardless of rolling distance. Further, the position (in the intermediate domains) of the contact region in both cylinder and plate is largely stationary throughout the rolling motion, allowing for localized mesh refinement. Figure 10 in Sect. 6 illustrates these points.
A potential difficulty in convective formulations is the tracking of boundaries. In the present case, since a round cylinder and a flat plate are considered, the boundaries are stationary and this is not a problem. However, numerical problems due to convective effects will require attention (see Sect. 4).

Thermomechanical problem
In the following, the strong and weak forms as well as the finite element formulation of the thermomechanically coupled problem are derived in terms of the ALE description. Isotropic, homogeneous and thermoelastic materials are initially assumed. Plane strain and linear elasticity are later assumed in the derivation of the FE formulation.

Strong form
The local balance of momentum equation with respect to a volume element in the initial configuration Ω is where P = P(F, θ) is the first Piola-Kirchhoff stress tensor, ρ is the density in the initial configuration, B is the external body force per unit volume (in the initial configuration) and a is the acceleration: . ∇ X is the vector differential operator with respect to Ω. Boundary conditions are imposed as (Γ := ∂Ω) where T := P · N (N is the outward normal), T P is the prescribed external traction per unit area with respect to Ω and u P represents prescribed displacements.
(·) denotes quantities related toΩ. In particular,P = P ·f T is the push-forward of the first Piola-Kirchhoff stress tensor toΩ. For a thermoelastic material, the constitutive relation P(F, θ) for the first Piola-Kirchhoff stress can be defined from a free energy function Ψ (F, θ) such that P(F, θ) = ∂Ψ (F, θ)/∂ F. For an isotropic material, Ψ should be independent of any rotation prior to the deformation, i.e.
for an arbitrary deformation gradient F, temperature θ and rotation tensor R. Consequently, the derivatives of Ψ with respect to F (i.e. the stress) must satisfy the conditions for anyF and any R. As a consequence, for anyf being either a rotation (as in the cylinder) or the identity (as in the plate). Hence, the original constitutive whereĤ =F − I and K is the bulk modulus. Hence, for an unconstrained specimen (P = 0), the temperature driven deformation isF = (1 + αθ)I.

Remark 1
In the cylinder domain, it is only the relationP(F) -and not P(F) -that can be linearized, since F will here contain a finite rotation (see Eqs. (3) and (1)).
The boundary conditions can in the ALE framework be phrased as whereT TOT in the Robin-type boundary condition (12a) is the natural boundary traction obtained via integration by parts in the weak form (see Sect. 3.2).
The strong form of the energy balance equation with respect to the initial configuration Ω is [13,14] (θ ref +θ)β : where is the deformation-induced heat source, c is the mass specific heat capacity, q is the heat flux and r is the external heat power per unit volume (in the initial configuration). Boundary conditions are imposed as where q N := q · N. In terms of the ALE description, the energy balance equation takes the form where( ·) denotes quantities related toΩ. In particular,β = β ·f T andq =f · q were introduced. Furthermore, for isotropic materials, where Eq. (10) was used and a known constitutive relation β(F, θ) in the initial configuration was assumed. Hence, where the last equality is valid for the linearized case (see Eq. (11)). Furthermore, in analogy with the result for the first Piola-Kirchhoff stress, it can be shown that for an isotropic material (and forf being a rotation or the identity tensor), where it was tacitly assumed that the heat flux only depends on the gradient of the temperature with respect to the initial configuration. In particular, the linear Fourier's law is henceforth adopted, whereby it is obtained that where k is the constant heat conductivity. The boundary conditions can in the ALE framework be phrased as For isotropic materials, (18) is valid, and we obtain where it was used thatv is divergence free (see Eq. (7)).
In summary: For homogeneous, isotropic materials, the residuals of the linearized strong form of the transient thermomechanically coupled problem are where the linearized form ofP is given in (11). It is clear from (11) thatP is temperature dependent, due to the influence of thermal expansion. This constitutes the influence of the temperature field on the momentum balance equation. Furthermore, the deformation-dependent terms in the second equation above represents the Gough-Joule effect: reversible heating/cooling of the material resulting from a nonzero strain rate [15]. In the ALE context, this term is split into a referential derivative and a convective term (as seen above). In a stationary analysis, the former vanishes. It should be noted that the Gough-Joule effect is negligible for thermoelastic materials [13]. Consequently, the thermomechanical coupling is one-sided in this case. It can be seen that when the translational and rotational velocity of the system is constant in time,v andẌ are constant in time (specifically,Ẍ = 0) and the time dependence in the above equations is confined to the solution fields (and the external loads).
If stationary rolling conditions are assumed, all referential time derivatives (d t , d tt ) as well asẌ are zero, resulting in a time-independent problem involving the residuals

Weak form
The weak form is obtained by weighting the local expressions in Eq. (24) with arbitrary (time-independent) test functions (δû, δθ) ∈ V 0 u × V 0 θ , integrating over the whole domainΩ and performing integration by parts. The weak residuals are thus defined as wherê wherê and it was used that q = −k∇θ.
Remark 2 As shown in Nackenhorst [9], it is possible to obtain a higher degree of symmetry in the weak form (26) by partial integration of the second term on the right-hand side: where the divergence theorem and the fact thatv is divergence free were used. It is noted that the first two terms on the right-hand side constitute an antisymmetric contribution to the weak form.
The termT TOT emerges from the ALE formulation of the momentum balance equation. Thus, prescribingT TOT on Γ N u constitutes a natural (Neumann) boundary condition: If instead the intrinsic (physical) tractionT :=P ·N is prescribed, a Robin-type boundary condition is obtained: where the boundaryΓ N u was simply renamedΓ R u in order to reflect the type of boundary condition in effect. Similarly, prescribing the quantityq TOT N onΓ N θ constitutes a natural (Neumann) boundary condition: If instead the intrinsic (physical) heat fluxq is prescribed, a Robin-type boundary condition is obtained: where the boundaryΓ N θ was renamedΓ R θ .
In order to state the final version of the weak form, trial and test spaces are introduced for the respective solution fieldŝ u(X, t) andθ(X, t): The exact meaning of sufficiently regular is not elaborated here (see eg. Brenner and Scott [16]). The weak form of the ALE boundary value problem derived in the previous section can now be stated as: where the residuals are obtained by inserting the aforementioned Robin-type boundary conditions into the integral expressions in Eqs. (26), (27):

Finite element formulation
A finite element formulation of the problem based on plane strain and linear elasticity is now presented. Voigt matrix notation is employed. Displacement and temperature fields are approximated by piecewise linear or piecewise quadratic functions: Shape function matrices for displacement and temperature are denoted by N u and N θ , respectively. Furthermore, B u :=∇ u N u and B θ :=∇ θ N θ , wherê Inserting solution field approximations and employing Galerkin test functions yields the FE formulation where Here, (E TOT ) i jkl = E i jkl −ρδ ikv jvl (E TOT is the Voigt matrix representation of E TOT ), 1 = [1 1 0 0] T , andv l and v r are the Voigt matrix representations of the tensors δ i jvk and δ i j (D tvk − (v kvl ) ,l ), respectively. In the above FE formulation, the "uθ "-and "θ u"-terms represent the thermomechanical coupling effects previously described under Eq. (24).
Recall from Sect. 3.1 that the time dependence of the problem is confined to the solution fields and the loads when the translational and rotational motion of the system is constant in time. It is clear from the above that this property is manifested as time-independent matrices in the FE formulation.
In the stationary case, the FE equations reduce to the timeindependent system

Element choice and numerical stabilization
Previous studies [12] indicated that, at least for rolling speeds up to the order of a few hundred km/h, no stability problems arise related to the discretized momentum balance equation for the present implementation. By contrast, various numerical problems have been found to have a prominent influence for the discretized energy balance equation, even for very modest rolling speeds. The numerical problems are manifested as node-to-node oscillations in the plate domain and oscillations together with a damped temperature profile in the cylinder domain. The former case, involving a uniform convective velocity field that intersects the boundaries of the domain, is well understood and amenable to a standard application of the SUPG method. The implementation and performance of this method will not be elaborated further in this paper, see instead the references below. The latter caseinvolving circular, closed convective streamlines-has been found to pose more of a challenge to these techniques, however.
The numerical stabilization methods employed in this paper is the SUPG method (see Donéa and Huerta [17] and Codina et al. [18] (for quadratic elements) ) and an approach employing a variant of residual-free bubbles [19,20].
The employed bubble function approach bears similarities to that of Brezzi et al. [21], and is implemented as follows: In each triangular element, an additional node is inserted, resulting in a subdivision into three triangles (see Fig. 2). The extra node is positioned along the directed line segment represented in the figure by a dashed arrow. This line segment passes through the element centroid (the black dot) and is parallel with the direction of the convective velocity evaluated at the centroid. The position of the node along this line segment is chosen to correspond to the stationary point of the analytical solution along the line segment, of the pertinent 1D convection-diffusion problem. The shape functions related to this node (the bubble functions) thus serve as a rough approximation of the shape of the solution in the element. Further, their support coincides with the given element, so static condensation can be used to keep the global degree of freedom set unchanged. The bubbles will be chosen as piecewise polynomials of the same order as that of the global approximation (even though it is possible to choose these functions independently).
The local subgrid is treated as a standard FE mesh-no two shape functions are nonzero at any node (not even the added node). This means that the stated subgrid enrichment scheme is equivalent to a standard Galerkin formulation featuring an enriched discrete function space (i.e. on a refined mesh) [22]. The choice of the positions of the extra nodes in this refined mesh have been informed by appropriate observations of the convective velocity field, as discussed above.
To illustrate the aforementioned numerical difficulties arising for sufficiently high convective velocities in the cylinder domain, a simple Eulerian formulation of stationary, pure heat transfer is considered (cf. Eq. (23)): where θ is the temperature, q is the heat flux and s is the external heat source. Further, where ρ is the density, c is the mass-specific heat capacity andv is the convective velocity. A two-dimensional annular where ω is the angular velocity, r is the distance from the center and e ϕ is the circumferential unit vector. The temperature is fixed to zero at both the inner and the outer boundary and the external heat source s is uniform. Figure 3 shows a schematic illustration of the considered problem. It is noted that the geometry, boundary conditions and loads result in a problem that is one-dimensional (radially symmetric) and has an analytical solution independent of ω [23]. Results presented below correspond to a finite element solution of Eq. (35). The mesh used (deliberately unstructured) is shown in Fig. 4. Gaussian quadrature of order 5 (7 integration points) is employed. Unless otherwise stated, the parameters used are as follows: Outer and inner radii of the cylinder r o = 50 cm and r i = 5 cm and external heat  Figure 5a, b show the temperature distribution along the radial line segment x = 0, y < 0 for varying ω and for linear and quadratic shape functions, respectively. Figures 6, 7 show the same thing, but implementing a SUPG stabilization method and the stated bubble function method, respectively.
The decay (due to spurious numerical dissipation) of the solution for an increasing rotation speed ω is clear from Fig.  5, as is the presence of numerical oscillations. As seen in Fig. 6, the SUPG method is able to smooth out the response, but unable to deal with the numerical damping effect: instead seemingly exacerbating it. The same can in general be said for the bubble function approach using linear elements (Fig.  7a), while the use of quadratic elements seems to work much better (see Fig. 7b). It should be noted that an integration order of at least four was found to be necessary in the latter case: below that, the solution started to decay.
In summary, the following choice of element and stabilization scheme has been found to be the most effective in reducing the influence of the aforementioned numerical issues: linear elements and a SUPG method in the plate domain (details . An integration order of four for the Gauss quadrature is used. Due to the similarity between the problem studied here and the rolling contact problem studied in Sect. 6, it is likely that the scheme described above will be suitable also for the latter, although a mesh convergence study is necessary to ascertain numerically robust results for the specific cases studied.

Remark 3 A comparison between
Figs. 5b and 7b shows that the implemented bubble function scheme serves to reduce both the spurious numerical oscillations and the amount of spurious numerical dissipation. It is plausible that an unbiased mesh refinement scheme-where the extra node is instead placed in the centroid of each element-would provide a stronger reduction of the numerical dissipation, but at the cost of a poorer ability of diminishing the numerical oscillations. The test problem studied here exhibits a high degree of symmetry, which means that the influence of numerical instability (oscillatory behaviour) is especially weak. A problem that is more true to life is likely to be less symmetrical and thus more sensitive to numerical instability. The need for the stabilizing effect of the bubble function scheme would be clearer in such a case.

Contact formulation
The employed contact formulation is presented below in the context of the ALE description. As the focus of this paper lies elsewhere than in realistic contact interface modelling, the simplest possible laws are chosen for this purpose. However, it is emphasized that the computational framework has been constructed with modularity and extensibility in mind. The implementation of more complex contact laws should therefore be straightforward. Examples of more advanced mechanical/thermal contact interface laws (for frictionless contact), based on microgeometrical and statistical considerations, can be found in [24][25][26][27][28]. Furthermore, [29] is noted, in which homogenization of thermal contact resistances is suggested.
The presentation given below of the mechanical and thermal components of the implemented contact formulation relies heavily on concepts described in Wriggers [8].

Mechanical contact formulation
A standard penalty method is employed for the formulation of normal mechanical contact. This involves applying a penalty traction at each pointX ∈Ω on the cylinder surface (here denoted the "slave" surface) that, in the deformed configurationω, penetrates the plate surface (here denoted the "master" surface). This traction is proportional to the gap function g (the negative of the penetration distance) and is directed normal to the master surface. An opposing traction is applied at the pointX m (X) =φ −1 (x m (x)) on the master surface, i.e. the point inΩ corresponding tox m (x): the normal projection ofx =φ(X) on the deformed master surface. Assigning superscripts "s" and "m" to terms related to the slave and master surfaces, respectively, enables formulation of the normal contact tractions as wheren is the normal of the deformed master surface at X m (X),Γ s cand is the candidate contact surface subset of the slave surface ∂Ω s and is the scalarn-component of the normal contact force acting on the slave surface. Here, is the penalty stiffness, is the gap function corresponding to the pair of pointsX and X m (X) and · are Macaulay brackets. Note that in theory, In a practical numerical implementation, an increased penalty stiffness leads to a decreased absolute value of the gap function, but an excessive increase leads to ill-conditioning of the discretized equation system. is typically taken as mesh-dependent, scaling inversely with some measure of the mesh size (thus having a higher value in more refined regions) [8]. For the purposes of the current implementation, the simple choice = n / h is deemed adequate, where n is a constant and h is a local measure of mesh size. Note that this choice results in the desired property lim h→0 = ∞ ⇒ lim h→0 g = 0. IfΓ s c ⊂Γ s cand is the contact surface-the subset of the slave surface corresponding to negative gap functions-the contribution to the weak form residual (30) due to the contact tractions is Introducing Galerkin test functions (δû → N u c u ) results in the following contribution to the left-hand side of the finite element formulation (33a): At this point, the discretization of the master and slave surfaces is introduced. Linear contact elements are used regardless of the order of the elements used for the discretization of the contacting bodies. This is for compability reasons in case the latter discretizations are of different orders. The above expression is in this context evaluated via a one-point quadrature scheme with the integration points coinciding with the nodes of the discretized slave surface (see Wriggers [8]). The result is the finite element load vector pertaining to the normal mechanical contact interaction, to be added to the right-hand side of Eq. (33)a: Here, g i = g(X a i ), n i = n(X a i ), W i are integration weights (related to edge element lengths) and {a i } n c i=1 is the active set: the set of nodes in contact.X m a i andX a i are the undeformed positions of the respective pointsx m a i andx a i in the deformed configuration, wherex m a i is the point on the master surface closest to the pointx a i . Figure 8 illustrates a contact element in the deformed discretized domain. Here,

Thermal contact formulation
The formulation of thermal conduction in the contact region is derived in a manner analogous to the case of normal mechanical contact. A "penalty" heat flux proportional to the temperature difference is imposed on surfaces in contact: where k c is the (velocity-independent) contact conductivity and is the temperature difference. Note that perfect thermal contact is represented by lim k c →∞ Δθ = 0. That is, the contact conductivity k c would play a role similar to that of the penalty stiffness in the mechanical normal contact formulation. However, k c should here in general be interpreted as a physical conductivity pertaining to the surface properties at the contact. In reality, this parameter exhibits a complex dependence on e.g. microgeometry, third body characteristics and contact pressure. However, following the stated ambition to keep the contact interface model as simple as possible, k c is taken as constant. This choice is adequate for the purposes of the present paper, but would obviously be an oversimplification in implementations striving for more realistic modelling of the contact interface. The contribution from the contact fluxes to the weak form residual (31) is Introducing Galerkin test functions (δθ → N θ c θ ) and employing one-point quadrature as above yields the finite element load vector pertaining to the interfacial thermal conduction, to be added to the right-hand side of Eq. (33)b: where, Δθ i = Δθ(X a i ).

Solution method
The addition of the (generally nonlinear) contact contributions (40) and (44) to the (otherwise linear) FE formulation of the thermomechanical boundary value problem (Eq. (33)) leads to a nonlinear equation system. This system is solved monolithically by the Newton method (which requires linearization of the contact contributions, although these expressions are not shown in this paper). The employed contact iteration scheme is identical to the one described in Wriggers [8]: In each iteration of the Newton solver, the residual vector and the tangent stiffness matrix are constructed, followed by an update of the solution guess. The construction of the residual and the tangent involves a contact search procedure, a central step of which being the identification of the set of active nodes (using the contact condition g i < 0).

Numerical investigations
6.1 Numerical model The following numerical examples are based on a 2D (plane strain) model of an annular cylinder rolling on a plate. Owing to the convective ALE description, the latter can be kept fairly short, regardless of the actual distance traversed by the cylinder during a simulation. As mentioned in Sect. 3, the model features an isotropic, homogeneous, linear elastic material. Pure rolling and constant rolling velocity (velocity of the cylinder centre relative to a fixed coordinate system) are assumed. A vertical distributed load is applied along the inner boundary of the cylinder. In addition, a constant normal heat flux is applied to given sections of the cylinder perimeter (each having an angular extension of 45 • ). Gravitational loads on the bodies are not included. The base of the plate is fixed in all degrees of freedom and the cylinder inner boundary is fixed in the horizontal direction. All other boundaries are free. The temperature at the plate ends is fixed (to the reference temperature) while all other boundaries (not in contact) are thermally insulated. Figure 9 shows a schematic illustration of the geometry, boundary conditions and loads in the employed model. It is noted that these exhibit vertical symmetry.
A standard parameter setup (used in the following simulations unless otherwise indicated) is now defined. Material parameters for cylinder and plate are: Young's modulus E = 200 GPa, Poisson's ratio ν = 0.3, specific heat capacity c = 460 Jkg −1 K −1 , thermal conductivity k = 45 Wm −1 K −1 , thermal expansion coefficient α = 4.8 · 10 −6 K −1 , density ρ = 8 · 10 3 kgm −3 . Outer and inner radii of the cylinder are r o = 50 cm and r i = 5 cm, respectively. The height of the plate is h = 0.1 m and the width of the plate domain is The rolling velocity is chosen as constant with magnitudē v = 50 km/h. The rotational velocity of the cylinder is then v/r o , due to the assumption of pure rolling.
As mentioned in the discussion following Eq. (24), the Gough-Joule effect is negligible for thermoelastic materials. It will therefore not be modeled in the following numerical examples.
Much of the following presentation will study the weighted mean temperatureθ m , which for an arbitrary domain V 0 can be defined as where |V 0 | is the volume of V 0 . In the following, V 0 will be chosen to represent the cylinder domain and the plate domain, alternatively. It should here be emphasized that in the latter case, the magnitude of the resulting mean temperature is rather arbitrary, since it depends strongly on the volume of the arbitrary domain over which the mean is taken.

Discretization
The finite element formulation of the problem is implemented in MATLAB. The element type used is a triangular element with two or one degrees of freedom per node in the mechanical and thermal problems, respectively. The approximation for both displacements and temperatures is piecewise linear Fig. 10 The employed finite element mesh, with a zoomed-in view of the refined contact region in the plate and piecewise quadratic in the cylinder. Note that the discrepancy in element order between cylinder and plate has implications for the contact formulation. The approach taken here is to regard both domains as linear as far as the contact formulation is concerned (as was discussed in Sect.

5.1).
The employed mesh is shown in Fig. 10. The mesh of the cylinder is constructed from a coarse basic mesh which is refined according to the following scheme: successive refinement in a series of gradually smaller domains centered at the point of initial contact → global refinements (2 are here used) → refinement of the largest elements along the periphery. After that, all nodes are remapped radially so that the outer nodes describe a circle. At this point, even though measures have been taken to ensure that the outer boundary is as round as possible despite local refinement near the contact region (via peripheral refinement), the centroid is inevitably slightly displaced compared to that of the enclosing circle. This has the effect of making the resultant of the centrifugal force vector (the second component of the vector f uv , see Eqs. (33), (26)) nonzero, which is unphysical for a circular domain. Further, this causes considerable errors in the contact computation. For this reason, the inner boundary is rigidly moved so that the position of the centroid of the discretized domain is corrected (thus slightly modifying the geometry). The positions of the inner nodes of the mesh are then determined by linear elastic equilibrium.
The plate domain mesh is constructed by starting from a structured mesh and refining it locally in a rectangular area centered at the point of initial contact. This area extends a distance 2d c horizontally and d c vertically, where d c is an analytical prediction (using Hertz theory) of the contact patch size. The minimum allowed size for an element in the final mesh of the plate is 1.5 times the size of the largest element in the cylinder mesh in the most refined region.
All mesh refinements are performed according to Rivara's longest-edge refinement technique [30], and the final mesh contains 8017 elements and 12608 nodes.

Numerical stabilization
SUPG stabilization is employed in the plate (linear elements) and a bubble function scheme (according to the description in Sect. 4) is used in the cylinder (quadratic elements). An integration order of four is employed for the Gauss quadrature.

Convergence study
A mesh convergence study is performed in which the number of global refinements of the cylinder mesh, n r , is varied. Figure 11 shows a plot of the temperature along the radial line segment x = 0, y < 0 (see Fig. 3) for varying values of n r . Note that the horizontal axis represents the dis-

Fig. 11
Mesh convergence study tance r from the cylinder center. The range is here limited to r o − d c /5 < r < r o . These results indicate that the employed choice n r = 2 results in a mesh that is sufficiently fine for the purposes of the subsequent numerical investigations. Figure 12 shows the stationary temperature distribution in two small regions close to the contact patch (each having the dimensions 2 mm × 1 mm). The regions are displaced to the left: the middle of the contact patch is highlighted by a vertical gray mark. The combined influence of heat flux across the contact interface and convective effects is here evident. In particular, the skewing influence of the latter, despite the symmetrical boundary conditions and loads, is noted. Figure 13 shows a comparison between the computed normal contact stress distribution and the analytical Hertzian solution [6,31]. It is clear from the close correspondence between the two curves that the influence of thermomechanical effects on the mechanical solution is not enough to visibly affect the contact stress distribution. Figure 14 shows weighted mean temperatures of cylinder and plate (calculated using Eq. (45) for each subdomain) versus the rolling speed. The figure shows that the model is successful in capturing the effect of an increased cooling rate of the cylinder with increasing rolling speed (even with a velocity-independent contact conductivity k c ). Note that the free boundary of the cylinder is insulated, implying that the cause of the cylinder mean temperature decrease can only be a higher heat flux through the contact interface, in turn caused by a decreased local temperature in the plate at the contact due to an increased convection in the plate. Furthermore, the latter phenomenon leads to a higher rate of heat extraction out of the modeled plate domain (which is bounded by the artificial edges denoted by C in Fig. 9). This explains the significant decrease in mean temperature in the plate with increasing rolling speed, seen in the figure.  Figure 15 (bottom) shows the temperature difference between the cylinder and plate nodes of initial contact. The temperatures in cylinder and plate are seen to approach each other (although the increase in temperature in the plate is very slight), and a saturation effect can be observed from the graphs. It can be concluded that the interval k c > 10 7 Wm −2 K −1 may be considered as resulting in a state of perfect thermal contact. A similar behaviour was observed in Vernersson [3], where a numerical model of wheel-rail heat transfer was used to model the rail chill effect on tread braked railway wheels.  temperature in the plate is very small), and is due to the increased contact patch width (as seen in the lower graph). Note that this effect would have been even stronger had the contact conductivity been modeled as pressure-dependent.

Transient analysis
For the transient simulations, a backward Euler time integration scheme is used. Three different scenarios involving transient processes are simulated: (i) the external heat flux W in is applied at time t = 0, (ii) the external heat flux is retracted at t = 0, (iii) a sharp hole of width 0.2 mm in the plate is traversed. In all cases, the appropriate stationary solution (i.e.  The plate profile used in case iii is implemented in the manner described in Draganis et al. [12], which in turn is based a methodology described in Thompson [32]. In this approach, a given offset function is imposed on the gap func- In cases i and ii, the time evolution of the weighted mean temperatures in both bodies is studied, while in case iii, the analysis is focused on the evolution of contact stresses and -forces as the hole in the plate is traversed. The former phenomenon occurs on a vastly larger time scale than the latter. The time steps used in the respective cases are i/ii: Δt = 500 s and iii: Δt = 10 −6 s. These time step sizes have been verified (details not given here) to give a numerically convergent response with respect to phenomena of interest in the present analysis. It is noted that the former time step size implies that the cylinder undergoes many revolutions per time step. In particular, this means that the simulation is unable to resolve mechanical phenomena in this case, which occur on vastly smaller time scales. However, their influence on the studied quantity in cases i and ii: the comparatively very slow evolution of the temperature distribution, is negligible. Taking into account also that the constitutive model does not include inelastic material parameters, it is concluded that the given choice of time step is admissible in this case. Figure 17 shows the time evolution of the weighted mean temperatures in cylinder and plate for the respective cases i, ii. Figure 17a also shows the stationary solution for reference. Due to the large difference in magnitude of the cylinder and plate temperatures (discussed in previous sections), a logarithmic axis is used for the vertical axis in these figures. The exponential decay of the temperature in 17b is noted. Figure 18 shows the contact stress distribution at four distinct points in time as the hole is traversed. Figure 19 shows the time evolution of the resultant contact force. Note that in order to resolve higher frequencies of the oscillations resulting from the contact interaction at the discontinuity, a smaller time step would be required.

Concluding remarks
A theoretical and computational framework governing thermomechanically coupled transient rolling contact based on an ALE kinematical description has been developed. The ALE formulation allows for linearization of the mechanical response, localized mesh refinement and a compact computational domain. Further, it was shown to simplify the time-description of the transient rolling contact problem and enable the formulation of the stationary rolling problem as time-independent.
Numerical simulations featuring both mechanical and thermal loads were performed. The results showed the thermomechanical contact model (featuring a velocity-independent contact conductivity) to be able to capture the effect of convective chilling of the cylinder due to the contact with the plate. A study of the influence of the contact conductivity was performed, and results were found to correspond qualitatively to results in the literature. Further, the relation between the magnitude of the heat flux through the contact interface and the applied mechanical load (owing to the influence of the latter on the contact patch width) was emphasized. Transient simulations showed the model to be able to capture phenomena occurring on disparate time scales, as well as simulations featuring very rough contact geometries.
The convective ALE formulation of the energy balance equation was found to be sensitive to stability problems and other numerical issues in its discretized form. Numerical stabilization techniques were implemented, satisfyingly addressing these problems. Since the numerical issues manifest themselves in essentially different ways in the two domains (due to the differences in the convective velocity fields), the numerical stabilization techniques implemented in these domains had to be designed thereafter.
Upcoming work will be focused on modelling frictional contact. Particular applications of interest include modelling of stick/slip phenomena and frictional heat generation.