# Engineering Software Solution for Thermal Elastohydrodynamic Lubrication Using Multiphysics Software

Authors: Thomas Lohner; Andreas Ziegltrum; Johann-Paul Stemplinger; Karsten  Stahl
Journal: Advances in Tribology
 (2016)

---

## Abstract

The complexity of thermal elastohydrodynamic lubrication (TEHL) problems has led to a variety of specialised numerical approaches ranging from finite difference based direct and inverse iterative methods such as Multilevel Multi-Integration solvers, via differential deflection methods, to finite element based full-system approaches. Hence, not only knowledge of the physical and technical relationships but also knowledge of the numerical procedures and solvers is necessary to perform TEHL simulations. Considering the state of the art of multiphysics software, the authors note the absence of a commercial software package for solving TEHL problems embedded in larger multiphysics software. By providing guidelines on how to implement a TEHL simulation model in commercial multiphysics software, the authors want to stimulate the research in computational tribology, so that, hopefully, the research focus can be shifted even more on physical modelling instead of numerical modelling. Validations, as well as result examples of the suggested TEHL model by means of simulated coefficients of friction, coated surfaces, and nonsmooth surfaces, highlight the flexibility and simplicity of the presented approach.


---

## Body

## 1. Introduction

Measurements, numerical simulations, and derived analytical solutions have, for decades, created detailed insights into the behaviour of thermal elastohydrodynamically lubricated (TEHL) contacts. Due to the large number of numerical investigations, it is astonishing that, so far, no commercial software package embedded in larger multiphysics software has been available to solve TEHL problems.The simulation of TEHL contacts is, mainly due to elastic deformations and the large pressure-induced increase of lubricant viscosity, highly nonlinear and has led to various specialised numerical simulation approaches. Great computational efforts, and the tendency to instabilities at high loads when using computational fluid dynamics (CFD) [1], make the use of simulations based on the Reynolds equation necessary. Hartinger et al. [2] thereby compare TEHL solutions obtained with CFD with solutions based on the Reynolds equation which show very good agreement. Lubrecht [3] and Venner [4] develop a direct iterative finite difference based Multilevel Multi-Integration (MLMI) solver for isothermal EHL problems. Habchi [1] introduces the finite element method (FEM) based full-system approach to solve the fully coupled system of the highly nonlinear elastohydrodynamic equations, which are coupled with the energy equations by an iterative procedure. Raisin et al. [5] extend the approach based on line contacts to transient operating conditions. Although both authors [5, 6] mention the use of the commercial software, COMSOL Multiphysics [7], to solve the TEHL problem, explicit implementation guidelines are not given. For the isothermal EHL line and point contact with Newtonian fluid behaviour under transient and steady-state operating conditions, Tan et al. [8] describe some information for implementation in COMSOL Multiphysics in more detail.Considering the large number of self-developed approaches and solvers, the modelling of TEHL contacts is hardly transparent. Hence, the promotion of using commercial software packages for solving TEHL problems is necessary to make simulation techniques available to a broader audience of the tribology society and to focus research more on physical relationships rather than on numerical procedures. Accordingly, the authors provide guidelines for the implementation and solving of TEHL problems in the commercial software, COMSOL Multiphysics, to stimulate and accelerate research in computational tribology. The possibility of developing a complex TEHL model with moderate effort is demonstrated in this paper.


## 2. Theory

In this section, the theory and governing equations for simulation of TEHL contacts are addressed. For simplification, a two-dimensional line contact is considered.


### 2.1. General

Contact conditions of, for example, gears are characterised by varying load, motion, and geometry along the path of contact. This transient contact can be interpreted as a transient model contact of two rolling elements with varying load, motion, and geometry for each meshing point. In TEHL simulations, the transient contact is usually further simplified to an equivalent contact between a single (inelastic) roller and an elastic flat body. This is adopted in this study and exemplarily shown for a line contact of a spur gear in Figure1.Figure 1
Line contact in a spur gear and its simplification to a model (disc) contact and to an equivalent contact.Characteristic quantities of the equivalent contact are described by the Hertzian contact parameters, for example,pH, bH, and Rx [9]. Poisson’s ratio νeq and Young’s Modulus Eeq of the equivalent elastic flat body are calculated according to Habchi [1]. The kinematics are characterised by the sum velocity vΣ and the siding velocity vg, whereby, for all analyses in this study, the velocity of the lower body is higher than the velocity of the upper body v1>v2:(1)vg=v1-v2,vΣ=v1+v2.The coefficient of friction μ of the TEHL contact is evaluated by integrating the shear stress in the middle of the lubricant film (see also Habchi [1]): (2)μ=∫xinxexτzxz=h/2dxFN/leff.


### 2.2. Generalised Reynolds Equation

Derived from the two-dimensional compressible Navier-Stokes equations by considering reasonable assumptions for TEHL contacts [1, 10] and no-slip boundary conditions, the shear rate γ˙x and velocity distribution vx in the lubricant film can be expressed by(3)γ˙xx,z,t=∂vx∂z=1η∂p∂xz-∫0hz/ηdz∫0h1/ηdz+v2-v1η∫0h1/ηdz,(4)vxx,z,t=∂p∂x∫0zz~ηdz~-∫0z1ηdz~∫0hz/ηdz∫0h1/ηdz+v2-v1∫0h1/ηdz∫0z1ηdz~+v1.Inserting the velocity equation in the integrated transient compressible continuity equation results in the generalised Reynolds equation according to Yang and Wen [11]:(5)∂∂x∂p∂x∫0hρ∫0zz~ηdz~-∫0z1ηdz~∫0hz/ηdz∫0h1/ηdzdz+∂∂xv2-v1∫0h1/ηdz∫0hρ∫0z1ηdz~dz+∫0hρdz·v1+∂∂t∫0hρdz=0.In TEHL line contacts, the generalised Reynolds equation determines the hydrodynamic pressure. Due to the integrations when deriving the generalised Reynolds equation, the pressure solution for line contacts is one-dimensional. The integral terms in (3) to (5) can be understood as integral viscosity and density terms. As a boundary condition, pxin=pxex=0 has to be prescribed. In addition, the cavitation model according to Wu [12] fulfils the Reynolds cavitation boundary condition [13]: p≥0 in ΩP and pxcav=∂p/∂xxcav=0.


### 2.3. Contact Mechanics

The calculation of the elastic deformation of the equivalent body is based on the finite element method as introduced by Habchi [1]. This has advantages, in terms of accuracy and resolution, over the half-space assumption [10]. The elastic deformation of the equivalent body in gap height direction δ is calculated by applying the linear elasticity equation and neglecting the dynamic response of the solids [5]:(6)∇·σ=0withσ=C·εU,δx,t=Uvx,twiththedisplacementvectorU=UuUv.Boundary conditions of the linear elasticity equation are zero displacement on the bottom of the calculation domain Ωδ and the hydrodynamic pressure p between xin and xex on the top of the calculation domain Ωδ, which is applied as a normal stress (Figure 3, top left). A free boundary condition assuming zero normal and tangential stress is applied to the remaining boundaries of the computational domain Ωδ.The film thickness equation describes the height of the separating lubricant film and consists of the constant parameterh0, the undeformed geometry, the elastic deformation of the equivalent body δ(x,t), and the deviation from the smooth profile R:(7)hx,t=h0+x22Rx+δ-R.Hydrodynamic pressure and the applied load have to be in balance. This is ensured by the load balance equation, which determines the constant parameter h0 in (7):(8)∫ΩPpx,tdΩP=FNlefft.


### 2.4. Energy Conservation

Similarly to the generalised Reynolds equation, the transient energy equations are simplified with reasonable assumptions typically used in TEHL calculations [1, 10]. The resulting energy equation of the lubricant can be written as(9)ρcp·∂T∂t+vx∂T∂x-∂∂xλ∂T∂x+∂∂zλ∂T∂z=-Tρ∂ρ∂T∂p∂t+vx∂p∂x+η∂vx∂z2.The transient energy equations of the solids (1,2) are shown in (10). The transitions between the solids and the lubricant are described by a temperature and conductive heat flux continuity condition:(10)ρicp,i·∂T∂t+vi∂T∂x-∂∂xλi∂T∂x+∂∂zλi∂T∂z=0withi=1,2.Thus, based on heat sources (shearing and compression of the lubricant), heat sinks (expansion of the lubricant), and heat transfer due to convection and conduction, the temperature distribution in lubricant and solids can be calculated. Boundary conditions ensure that lubricant entering the computational domain ΩT (vxxin,z,t≥0 and vxxex,z,t≤0) has bulk temperature TM and for lubricant leaving ΩT (vxxin,z,t≤0 and vxxex,z,t≥0) the conductive heat flux is zero (-λ·∂T/∂x=0). Similarly, solids entering the computational domains ΩT,1 and ΩT,2 (vixin,t≥0 and vixex,t≤0) have bulk temperature TM and for solids leaving ΩT,1 and ΩT,2 (vixin,t≤0 and vixex,t≥0) the conductive heat flux is zero (-λi·∂T/∂x=0). The height of solids ΩT,1 and ΩT,2 is with d=3.15·bH [14] sufficient so that the conductive heat flux into the solids becomes zero and bulk temperature TM can be assumed at d.


### 2.5. Lubricant Properties

Pressure, temperature, and shear rate distributions in TEHL contacts have a significant influence on lubricant properties. The accompanying changes in viscosity have the largest influence on the TEHL contact itself. The models for pressure and temperature dependency of the viscosity as suggested in [15] are used. The dynamic viscosity for a given temperature and ambient pressure ηT is described by the Vogel model [16]:(11)ηT=A·exp⁡BC+T-273.15K.The pressure dependency of the viscosity ηT,p is modelled by the Roelands [17] equation(12)ηT,p=ηT·exp⁡ln⁡ηT+9.67·-1+1+ppη0zTwith a temperature dependent pressure exponent z(T) according to(13)zT=αpT·pη0ln⁡ηT+9.67.The pressure-viscosity coefficient αp(T) is temperature dependent:(14)αpT=Eαp1·exp⁡Eαp2·T.The rheological behaviour of the lubricant is described by an Eyring model [18]:(15)ηT,p,γ˙x=τcγ˙xasinh⁡ηT,p·γ˙xτcwithτT,p,γ˙x=ηT,p,γ˙·γ˙x.The temperature and pressure dependency of lubricant density ρT,p is modelled by the Bode model [19]:(16)ρT,p=ρs·1-αs·T1-Dρ0·ln⁡Dρ1+Dρ2·T+p/Dρ1+Dρ2·T.Based on measurements, Larsson and Andersson [20] derived a model for the temperature and pressure dependency of the thermal conductivity λ(p) and of the specific heat capacity per volume (cp·ρ)(T,p). Temperature dependency of the thermal conductivity is not considered; as in [20], the influence is assessed as negligible. The equations for λ(p) and (cp·ρ)(T,p) are(17)λp=λ0·1+dλ1·p1+dλ2·p,(18)cp·ρT,p=cp·ρ295K·1+Ac1·p1+Ac2·p1+Ac3·1+Ac4·p+Ac5·p2·T-295K.Among others, parameters for paraffinic mineral oil are available for (17) and (18) in [20].


## 2.1. General

Contact conditions of, for example, gears are characterised by varying load, motion, and geometry along the path of contact. This transient contact can be interpreted as a transient model contact of two rolling elements with varying load, motion, and geometry for each meshing point. In TEHL simulations, the transient contact is usually further simplified to an equivalent contact between a single (inelastic) roller and an elastic flat body. This is adopted in this study and exemplarily shown for a line contact of a spur gear in Figure1.Figure 1
Line contact in a spur gear and its simplification to a model (disc) contact and to an equivalent contact.Characteristic quantities of the equivalent contact are described by the Hertzian contact parameters, for example,pH, bH, and Rx [9]. Poisson’s ratio νeq and Young’s Modulus Eeq of the equivalent elastic flat body are calculated according to Habchi [1]. The kinematics are characterised by the sum velocity vΣ and the siding velocity vg, whereby, for all analyses in this study, the velocity of the lower body is higher than the velocity of the upper body v1>v2:(1)vg=v1-v2,vΣ=v1+v2.The coefficient of friction μ of the TEHL contact is evaluated by integrating the shear stress in the middle of the lubricant film (see also Habchi [1]): (2)μ=∫xinxexτzxz=h/2dxFN/leff.


## 2.2. Generalised Reynolds Equation

Derived from the two-dimensional compressible Navier-Stokes equations by considering reasonable assumptions for TEHL contacts [1, 10] and no-slip boundary conditions, the shear rate γ˙x and velocity distribution vx in the lubricant film can be expressed by(3)γ˙xx,z,t=∂vx∂z=1η∂p∂xz-∫0hz/ηdz∫0h1/ηdz+v2-v1η∫0h1/ηdz,(4)vxx,z,t=∂p∂x∫0zz~ηdz~-∫0z1ηdz~∫0hz/ηdz∫0h1/ηdz+v2-v1∫0h1/ηdz∫0z1ηdz~+v1.Inserting the velocity equation in the integrated transient compressible continuity equation results in the generalised Reynolds equation according to Yang and Wen [11]:(5)∂∂x∂p∂x∫0hρ∫0zz~ηdz~-∫0z1ηdz~∫0hz/ηdz∫0h1/ηdzdz+∂∂xv2-v1∫0h1/ηdz∫0hρ∫0z1ηdz~dz+∫0hρdz·v1+∂∂t∫0hρdz=0.In TEHL line contacts, the generalised Reynolds equation determines the hydrodynamic pressure. Due to the integrations when deriving the generalised Reynolds equation, the pressure solution for line contacts is one-dimensional. The integral terms in (3) to (5) can be understood as integral viscosity and density terms. As a boundary condition, pxin=pxex=0 has to be prescribed. In addition, the cavitation model according to Wu [12] fulfils the Reynolds cavitation boundary condition [13]: p≥0 in ΩP and pxcav=∂p/∂xxcav=0.


## 2.3. Contact Mechanics

The calculation of the elastic deformation of the equivalent body is based on the finite element method as introduced by Habchi [1]. This has advantages, in terms of accuracy and resolution, over the half-space assumption [10]. The elastic deformation of the equivalent body in gap height direction δ is calculated by applying the linear elasticity equation and neglecting the dynamic response of the solids [5]:(6)∇·σ=0withσ=C·εU,δx,t=Uvx,twiththedisplacementvectorU=UuUv.Boundary conditions of the linear elasticity equation are zero displacement on the bottom of the calculation domain Ωδ and the hydrodynamic pressure p between xin and xex on the top of the calculation domain Ωδ, which is applied as a normal stress (Figure 3, top left). A free boundary condition assuming zero normal and tangential stress is applied to the remaining boundaries of the computational domain Ωδ.The film thickness equation describes the height of the separating lubricant film and consists of the constant parameterh0, the undeformed geometry, the elastic deformation of the equivalent body δ(x,t), and the deviation from the smooth profile R:(7)hx,t=h0+x22Rx+δ-R.Hydrodynamic pressure and the applied load have to be in balance. This is ensured by the load balance equation, which determines the constant parameter h0 in (7):(8)∫ΩPpx,tdΩP=FNlefft.


## 2.4. Energy Conservation

Similarly to the generalised Reynolds equation, the transient energy equations are simplified with reasonable assumptions typically used in TEHL calculations [1, 10]. The resulting energy equation of the lubricant can be written as(9)ρcp·∂T∂t+vx∂T∂x-∂∂xλ∂T∂x+∂∂zλ∂T∂z=-Tρ∂ρ∂T∂p∂t+vx∂p∂x+η∂vx∂z2.The transient energy equations of the solids (1,2) are shown in (10). The transitions between the solids and the lubricant are described by a temperature and conductive heat flux continuity condition:(10)ρicp,i·∂T∂t+vi∂T∂x-∂∂xλi∂T∂x+∂∂zλi∂T∂z=0withi=1,2.Thus, based on heat sources (shearing and compression of the lubricant), heat sinks (expansion of the lubricant), and heat transfer due to convection and conduction, the temperature distribution in lubricant and solids can be calculated. Boundary conditions ensure that lubricant entering the computational domain ΩT (vxxin,z,t≥0 and vxxex,z,t≤0) has bulk temperature TM and for lubricant leaving ΩT (vxxin,z,t≤0 and vxxex,z,t≥0) the conductive heat flux is zero (-λ·∂T/∂x=0). Similarly, solids entering the computational domains ΩT,1 and ΩT,2 (vixin,t≥0 and vixex,t≤0) have bulk temperature TM and for solids leaving ΩT,1 and ΩT,2 (vixin,t≤0 and vixex,t≥0) the conductive heat flux is zero (-λi·∂T/∂x=0). The height of solids ΩT,1 and ΩT,2 is with d=3.15·bH [14] sufficient so that the conductive heat flux into the solids becomes zero and bulk temperature TM can be assumed at d.


## 2.5. Lubricant Properties

Pressure, temperature, and shear rate distributions in TEHL contacts have a significant influence on lubricant properties. The accompanying changes in viscosity have the largest influence on the TEHL contact itself. The models for pressure and temperature dependency of the viscosity as suggested in [15] are used. The dynamic viscosity for a given temperature and ambient pressure ηT is described by the Vogel model [16]:(11)ηT=A·exp⁡BC+T-273.15K.The pressure dependency of the viscosity ηT,p is modelled by the Roelands [17] equation(12)ηT,p=ηT·exp⁡ln⁡ηT+9.67·-1+1+ppη0zTwith a temperature dependent pressure exponent z(T) according to(13)zT=αpT·pη0ln⁡ηT+9.67.The pressure-viscosity coefficient αp(T) is temperature dependent:(14)αpT=Eαp1·exp⁡Eαp2·T.The rheological behaviour of the lubricant is described by an Eyring model [18]:(15)ηT,p,γ˙x=τcγ˙xasinh⁡ηT,p·γ˙xτcwithτT,p,γ˙x=ηT,p,γ˙·γ˙x.The temperature and pressure dependency of lubricant density ρT,p is modelled by the Bode model [19]:(16)ρT,p=ρs·1-αs·T1-Dρ0·ln⁡Dρ1+Dρ2·T+p/Dρ1+Dρ2·T.Based on measurements, Larsson and Andersson [20] derived a model for the temperature and pressure dependency of the thermal conductivity λ(p) and of the specific heat capacity per volume (cp·ρ)(T,p). Temperature dependency of the thermal conductivity is not considered; as in [20], the influence is assessed as negligible. The equations for λ(p) and (cp·ρ)(T,p) are(17)λp=λ0·1+dλ1·p1+dλ2·p,(18)cp·ρT,p=cp·ρ295K·1+Ac1·p1+Ac2·p1+Ac3·1+Ac4·p+Ac5·p2·T-295K.Among others, parameters for paraffinic mineral oil are available for (17) and (18) in [20].


## 3. Model Description and Implementation

A guideline for the implementation of a TEHL simulation model in COMSOL Multiphysics (abbr. by COMSOL) and MathWorks MATLAB [21] (abbr. by MATLAB) is given in the following sections. All COMSOL-specific notations are written in italics.


### 3.1. Dimensionless Parameters

For a convenient numerical solution procedure of the nonlinear system of equations with good conditioning, the variables are transferred to a dimensionless form. Following Habchi [1], Tan et al. [8], Venner [4], and Lubrecht [3], the definitions are(19)
X
=
x
b
H
,
T
-
=
T
T
M
,
D
=
d
b
H
,
P
=
p
p
H
,
H
=
h
·
R
x
b
H
2
,
δ
-
=
δ
·
R
x
b
H
2
,
ρ
-
=
ρ
ρ
T
M
,
η
-
=
η
η
T
M
,
t
-
=
t
·
v
Σ
2
b
H
,
C
v
=
v
Σ
t
v
Σ
t
=
0
,
C
R
=
R
x
t
R
x
t
=
0
,
C
w
=
F
N
/
l
e
f
f
t
F
N
/
l
e
f
f
t
=
0
,
Z
=
z
b
H
f
o
r
s
o
l
i
d
s
z
h
f
o
r
l
u
b
r
i
c
a
n
t
.


### 3.2. Numerical Solution Procedure

The numerical solution scheme shown in Figure2 is similar to other schemes of, for example, Habchi [1] and Wang et al. [22]. In this work, it is based on two self-developed decoupled COMSOL FEM-models, denoted as FEM-model (P,H) and FEM-model (T-), and a MATLAB sequence controller, which also performs additional calculations.Figure 2
Numerical solution scheme.Figure 3
FEM-model (P,H) for pressure and film thickness calculation with corresponding COMSOL model tree (left) and FEM-model (T-) for temperature calculation with corresponding COMSOL model tree (right); element numbers are adjusted for presentation purposes.After reading all the required input parameters, an initial solution based on a simple steady-state isothermal Newtonian approach is calculated in the FEM-model (P,H) for the initial operating conditions. The procedure thereby is the same as described in the following paragraph. Based on the initial solution, the time loop is launched and repeated until the last time-step has been calculated. The coupling of two consecutive time-steps is realised by an implicit backward differentiation formula (BDE) scheme of the first order as used in [23]. It is implemented by passing the required solution variables of the previous time-step for calculating the time derivatives of the generalised Reynolds equation and energy equations from the sequence controller to the FEM-model (P,H) and FEM-model (T-).Within the time loop, a global loop is launched for each time-step. Firstly the pressure and film thickness distribution under consideration of non-Newtonian fluid behaviour is calculated for a given temperature distributionT- in the FEM-model (P,H) described in Section 3.3.1. Then, the lubricant properties acc. to Section 2.5, the integral terms in (5), and the velocity distribution vx (4) are updated in the sequence controller. Thereby, the shear stress τ is analytically derived (15). After calling up the FEM-model (P,H) again, the new calculated pressure distribution is compared with the previous calculated pressure distribution. The iterative procedure is repeated until the maximum absolute difference of two consecutive pressure distributions is smaller than 10-3.After the convergence ofP, the pressure and film thickness distributions P and H are kept constant and an iterative calculation of the temperature distribution T- in lubricant and solid bodies, based on the FEM-model (T-) described in Section 3.3.2, is launched. For each new temperature distribution T- calculated, the lubricant properties acc. to Section 2.5, the integral terms in (5), the heat sources in (9), and the velocity distribution vx (4) are updated in the sequence controller. Thereby, the shear stress τ is analytically derived (15). The iterative calculation is repeated until the maximum absolute difference of two consecutive temperature distributions in the lubricant is smaller than 10-3.Finally, the converged pressure distribution and temperature distribution in the lubricant of two consecutive global loops are compared. Again, convergence is assumed when the maximum absolute difference between two consecutive solutions of pressure and temperature distributions is smaller than10-3. After saving the solution of the considered time-step, the next global loop is calculated with updated operating conditions and time derivatives of the generalised Reynolds equation and energy equations.


### 3.3. Implementation in COMSOL

In the following, the characteristics of the FEM-model (P,H) and FEM-model (T-), as well as of the solver settings and sequence controller, are described in detail. To a great extent, the structure of the FEM-models follows the outstanding work of Habchi [1], who introduces the FEM based full-system approach to TEHL achieving very high convergence rates without underrelaxation required. Thereby, the FEM approach offers the possibility of using a mesh with nonregular unstructured elements and high order approximation functions. As shown in the following, it is possible to implement the approach of Habchi in commercial multiphysics software.


#### 3.3.1. FEM-Model (P,H)

Figure3 (top left) shows the computational domains Ωδ and ΩP on which the pressure and film thickness distribution are calculated fully coupled. Figure 3 (bottom left) illustrates in parts the corresponding model tree. For each new call-up of FEM-model (P,H), the solution of the previous call is used as an initial value.The elastic deformation of the equivalent body inZ-direction δ-(X,t-) is calculated on a quadratic two-dimensional computational domain Ωδ with the extensions X∈-30,30 and Z∈-60,0. Hence, an influence of the boundary conditions of Ωδ on the contact area can be excluded. Thephysics Solid Mechanics is used with the displacement vector U-(X,t-) as the dependent variable. Alinear elastic isotropic material specified by Young’s Modulus Eeq and Poisson’s ratio νeq is assumed. The boundary conditions described in Section 2.3 are implemented as afixed constraint condition on the bottom boundary of Ωδ, aboundary load condition with theload type pressureP on ΩP, and afree boundary condition on the remaining boundaries of Ωδ. Lagrange quadratic elements are used for the linear elasticity equations.The pressure and film thickness distributionsP and H, as well as the constant parameter H0 based on the load balance in (8), are calculated on the one-dimensional computational domain ΩP with an extension from Xin=-4.5 to Xex=3. Thereby, Xex of ΩP has already been enlarged for simulations with varying load over time, for example, [5]. The dependent variables are the pressure P, the constant parameter H0, and the displacement vector U-. The film thickness in (7) is implemented as avariable and the deviation from the smooth profile R is read as the input parameter by aninterpolation function. The load balance in (8) is implemented as aGlobal Equation with acomponent coupling integration on ΩP for the left-hand side of the load balance equation. Its result is the constant parameter H0. As no specific module for the generalised Reynolds equation in (5) is available in COMSOL, it is defined in thephysics Weak Form Boundary PDE by using Lagrange quintic elements. Based on the dimensionless parameters in (19), the weak form of the generalised Reynolds equation in dimensionless form is written as(20)-∫ΩPε-+ID∂P∂X∂WP∂XdX-∫ΩPρ-HH∂WP∂XdX+∫ΩP∂ρ-Ht-∂t-WPdΩP+∫ΩPQξWPdΩP+SUPGGLS=0withSUPGGLS=SUPGorGLS.WP is the test function for the pressure P and is denoted as test(P) in COMSOL. ID and SUPG/GLS represent numerical stabilisation terms, which are implemented asvariables and given here for completeness. Note that, for moderate cases, the TEHL simulation model can work without stabilisation terms. The numerical stabilisation terms have been applied to the generalised Reynolds equation by Habchi [1, 24] in order to avoid the oscillatory behaviour of the Reynolds equation at high loads. SUPG (Streamline-Upwind-Petrov-Galerkin) andGLS (Galerkin-Least-Squares) are consistent stabilisation terms whereasID (Isotropic Diffusion) is a nonconsistent stabilisation term. For implementation of ID and SUPG/GLS, the reader is referred to [1, 24]. The penalty term Qξ=Pf·P- as an additional source term with Pf in the order of 106 [1] ensures the Reynolds cavitation boundary conditions [13] and represents the cavitation model according to Wu [12]. Zero pressure on the boundaries of ΩP is ensured by aDirichlet boundary condition. The coefficients ε-, ρ-HH, and ρ-Ht- represent the first, second, and third term of (5) in dimensionless form:(21)ε-=H3λ-∫01ρ-∫0ZZ~η-dZ~-∫0Z1η-dZ~∫01Z/η-dZ∫011/η-dZdZwithλ-=Rx2·ηTM·vΣ2·pH·bH3,ρ-HH=2HvΣv2-v1∫011/η-dZ∫01ρ-∫0Z1η-dZ~dZ+∫01ρ-v1dZ,ρ-Ht-=H∫01ρ-dZ.Note thereby that all the parts of the coefficients in (21) and of the stabilisation terms in (20) can, despite the film thickness H, be prepared in the sequence controller and transferred to the FEM-model (P,H) byinterpolation functions.The mesh in the FEM-model (P,H) isfree triangular with a normal element size and a refinement at ΩP by adistribution with typically 1000 elements.


#### 3.3.2. FEM-Model (T-)

Figure3 (top right) shows the computational domains ΩT, ΩT,1, and ΩT,2 on which the temperature distribution T- in lubricant and solid bodies is calculated fully coupled. Figure 3 (bottom right) illustrates in parts the corresponding model tree. The lubricant domain ΩT has the dimensionless height of 1 with Z∈0,1, whereas the solid domains ΩT,1 and ΩT,2 have the dimensionless heights of D=3.15 with Z∈-D,0 and Z∈1,D+1, respectively. The length of the computational domains extends from Xin to Xex (see Figure 3, top left). For each new call-up of the FEM-model (T-), the solution of the previous call-up is used as initial value. The transient energy equations in dimensionless form are implemented and calculated in thephysics Heat Transfer in Fluids for the lubricant domain ΩT and inHeat Transfer in Solids for the solid domains ΩT,1 and ΩT,2. The dependent variable T- is solved based on Lagrange quadratic elements. The numerical stabilisation terms ID, SUPG, and GLS are already implemented in the consideredheat transfer physics in COMSOL and applied when required. The velocity distribution (4), viscosity distribution (15), density distribution (16), thermal conductivity distribution (17), and heat capacity distribution (18) in ΩT are calculated in the sequence controller and transferred to the FEM-model (T-) byinterpolation functions. Furthermore, theheat sources in thephysics Heat Transfer in Fluids, due to the shearing and compression of the lubricant, are calculated in the sequence controller and transferred to FEM-model (T-). Note thereby that, similar to the coefficients in (21), T- as a dependent variable needs to be applied in the FEM-model (T-). The implementation of the energy equations on ΩT,1 and ΩT,2 in theHeat Transfer in Solids modules with the translation motions of solids v1 and v2 as input parameters is straightforward.The boundary conditions in Section2.4 are implemented by defining bulk temperature T-M=1 at the left and upper boundary of ΩT,2 and at the left and lower boundary of ΩT,1, anOutflow boundary condition at the right-hand boundaries of ΩT,1 and ΩT,2 and anOpen boundary condition at the left- and right-hand boundary of the lubricant domain ΩT with lubricant entering ΩT having a bulk temperature T-M=1. As a common coordinate system in the FEM-model (T-) is used, manual adjustments are necessary to respect the different dimensionless definitions of Z in ΩT and ΩT,1 and ΩT,2, respectively (19). To be precise, the dimensionless scale of the Z axis of the solid domains ΩT,1 and ΩT,2 has to be adjusted according to the dimensionless scale of the Z axis of the lubricant domain ΩT, so that the continuity condition of the conductive heat flux and temperature across the lubricant-solid boundaries is respected. This is handled in COMSOL by multiplying the dimensionless energy equations of the solids by the conversion factor Zcf. This factor is the ratio of the dimensionless scales of the Z axes of the lubricant domain and the solid domains:(22)Zcf=ZlubricantZsolid=z/hz/bH=RxH·bH.The mesh in the FEM-model (T-) is divided into the lubricant domain ΩT and the solid domains ΩT,1 and ΩT,2. Whereas ΩT is meshed by amapped equidistant mesh, ΩT,1 and ΩT,2 are meshed withfree triangularelements with a refinement on the solid-lubricant boundaries and anextra-fine element size. The lubricant domain ΩT may consist of 40 elements in height and 1000 elements in the length direction.


#### 3.3.3. Sequence Controller

The communication between the MATLAB sequence controller and the FEM-model (P,H) and FEM-model (T-) is achieved by COMSOL’sLiveLink for MATLAB function. The MATLAB sequence controller basically represents the implementation of the global numerical scheme shown in Figure 2. There, it is very convenient, for example, to implement the models for the lubricant properties, calculate the integrals by adaptive Simpson quadrature rule with typically 40 elements in Z-direction, and perform other calculations with the current values for pressure, film thickness, and temperature. As no other possibility is available, the data exchange of vectors and matrices from MATLAB to COMSOL is handled in a somewhat cumbersome manner: The quantities, calculated in the sequence controller, for example, velocity distribution, are written into a text file and read by the COMSOLinterpolation functions. For data exchange from COMSOL to MATLAB, thempheval command is used.


#### 3.3.4. Solver

For solving the nonlinear system of equations, two separate COMSOL solvers for the FEM-model (P,H) and for the FEM-model (T-) are required. The solver for the FEM-model (P,H) includes the dependent solution variables P, H0, and U-, whereas the solver for the FEM-model (T-) includes the temperature T- as the dependent solution variable. For both FEM-models, a direct COMSOLMUMPS (Multifrontal Massively Parallel Direct Solver) [25] is used with a damped Newton-Raphson approach [26] and 10-3 as the tolerance factor as defined in [7].


## 3.1. Dimensionless Parameters

For a convenient numerical solution procedure of the nonlinear system of equations with good conditioning, the variables are transferred to a dimensionless form. Following Habchi [1], Tan et al. [8], Venner [4], and Lubrecht [3], the definitions are(19)
X
=
x
b
H
,
T
-
=
T
T
M
,
D
=
d
b
H
,
P
=
p
p
H
,
H
=
h
·
R
x
b
H
2
,
δ
-
=
δ
·
R
x
b
H
2
,
ρ
-
=
ρ
ρ
T
M
,
η
-
=
η
η
T
M
,
t
-
=
t
·
v
Σ
2
b
H
,
C
v
=
v
Σ
t
v
Σ
t
=
0
,
C
R
=
R
x
t
R
x
t
=
0
,
C
w
=
F
N
/
l
e
f
f
t
F
N
/
l
e
f
f
t
=
0
,
Z
=
z
b
H
f
o
r
s
o
l
i
d
s
z
h
f
o
r
l
u
b
r
i
c
a
n
t
.


## 3.2. Numerical Solution Procedure

The numerical solution scheme shown in Figure2 is similar to other schemes of, for example, Habchi [1] and Wang et al. [22]. In this work, it is based on two self-developed decoupled COMSOL FEM-models, denoted as FEM-model (P,H) and FEM-model (T-), and a MATLAB sequence controller, which also performs additional calculations.Figure 2
Numerical solution scheme.Figure 3
FEM-model (P,H) for pressure and film thickness calculation with corresponding COMSOL model tree (left) and FEM-model (T-) for temperature calculation with corresponding COMSOL model tree (right); element numbers are adjusted for presentation purposes.After reading all the required input parameters, an initial solution based on a simple steady-state isothermal Newtonian approach is calculated in the FEM-model (P,H) for the initial operating conditions. The procedure thereby is the same as described in the following paragraph. Based on the initial solution, the time loop is launched and repeated until the last time-step has been calculated. The coupling of two consecutive time-steps is realised by an implicit backward differentiation formula (BDE) scheme of the first order as used in [23]. It is implemented by passing the required solution variables of the previous time-step for calculating the time derivatives of the generalised Reynolds equation and energy equations from the sequence controller to the FEM-model (P,H) and FEM-model (T-).Within the time loop, a global loop is launched for each time-step. Firstly the pressure and film thickness distribution under consideration of non-Newtonian fluid behaviour is calculated for a given temperature distributionT- in the FEM-model (P,H) described in Section 3.3.1. Then, the lubricant properties acc. to Section 2.5, the integral terms in (5), and the velocity distribution vx (4) are updated in the sequence controller. Thereby, the shear stress τ is analytically derived (15). After calling up the FEM-model (P,H) again, the new calculated pressure distribution is compared with the previous calculated pressure distribution. The iterative procedure is repeated until the maximum absolute difference of two consecutive pressure distributions is smaller than 10-3.After the convergence ofP, the pressure and film thickness distributions P and H are kept constant and an iterative calculation of the temperature distribution T- in lubricant and solid bodies, based on the FEM-model (T-) described in Section 3.3.2, is launched. For each new temperature distribution T- calculated, the lubricant properties acc. to Section 2.5, the integral terms in (5), the heat sources in (9), and the velocity distribution vx (4) are updated in the sequence controller. Thereby, the shear stress τ is analytically derived (15). The iterative calculation is repeated until the maximum absolute difference of two consecutive temperature distributions in the lubricant is smaller than 10-3.Finally, the converged pressure distribution and temperature distribution in the lubricant of two consecutive global loops are compared. Again, convergence is assumed when the maximum absolute difference between two consecutive solutions of pressure and temperature distributions is smaller than10-3. After saving the solution of the considered time-step, the next global loop is calculated with updated operating conditions and time derivatives of the generalised Reynolds equation and energy equations.


## 3.3. Implementation in COMSOL

In the following, the characteristics of the FEM-model (P,H) and FEM-model (T-), as well as of the solver settings and sequence controller, are described in detail. To a great extent, the structure of the FEM-models follows the outstanding work of Habchi [1], who introduces the FEM based full-system approach to TEHL achieving very high convergence rates without underrelaxation required. Thereby, the FEM approach offers the possibility of using a mesh with nonregular unstructured elements and high order approximation functions. As shown in the following, it is possible to implement the approach of Habchi in commercial multiphysics software.


### 3.3.1. FEM-Model (P,H)

Figure3 (top left) shows the computational domains Ωδ and ΩP on which the pressure and film thickness distribution are calculated fully coupled. Figure 3 (bottom left) illustrates in parts the corresponding model tree. For each new call-up of FEM-model (P,H), the solution of the previous call is used as an initial value.The elastic deformation of the equivalent body inZ-direction δ-(X,t-) is calculated on a quadratic two-dimensional computational domain Ωδ with the extensions X∈-30,30 and Z∈-60,0. Hence, an influence of the boundary conditions of Ωδ on the contact area can be excluded. Thephysics Solid Mechanics is used with the displacement vector U-(X,t-) as the dependent variable. Alinear elastic isotropic material specified by Young’s Modulus Eeq and Poisson’s ratio νeq is assumed. The boundary conditions described in Section 2.3 are implemented as afixed constraint condition on the bottom boundary of Ωδ, aboundary load condition with theload type pressureP on ΩP, and afree boundary condition on the remaining boundaries of Ωδ. Lagrange quadratic elements are used for the linear elasticity equations.The pressure and film thickness distributionsP and H, as well as the constant parameter H0 based on the load balance in (8), are calculated on the one-dimensional computational domain ΩP with an extension from Xin=-4.5 to Xex=3. Thereby, Xex of ΩP has already been enlarged for simulations with varying load over time, for example, [5]. The dependent variables are the pressure P, the constant parameter H0, and the displacement vector U-. The film thickness in (7) is implemented as avariable and the deviation from the smooth profile R is read as the input parameter by aninterpolation function. The load balance in (8) is implemented as aGlobal Equation with acomponent coupling integration on ΩP for the left-hand side of the load balance equation. Its result is the constant parameter H0. As no specific module for the generalised Reynolds equation in (5) is available in COMSOL, it is defined in thephysics Weak Form Boundary PDE by using Lagrange quintic elements. Based on the dimensionless parameters in (19), the weak form of the generalised Reynolds equation in dimensionless form is written as(20)-∫ΩPε-+ID∂P∂X∂WP∂XdX-∫ΩPρ-HH∂WP∂XdX+∫ΩP∂ρ-Ht-∂t-WPdΩP+∫ΩPQξWPdΩP+SUPGGLS=0withSUPGGLS=SUPGorGLS.WP is the test function for the pressure P and is denoted as test(P) in COMSOL. ID and SUPG/GLS represent numerical stabilisation terms, which are implemented asvariables and given here for completeness. Note that, for moderate cases, the TEHL simulation model can work without stabilisation terms. The numerical stabilisation terms have been applied to the generalised Reynolds equation by Habchi [1, 24] in order to avoid the oscillatory behaviour of the Reynolds equation at high loads. SUPG (Streamline-Upwind-Petrov-Galerkin) andGLS (Galerkin-Least-Squares) are consistent stabilisation terms whereasID (Isotropic Diffusion) is a nonconsistent stabilisation term. For implementation of ID and SUPG/GLS, the reader is referred to [1, 24]. The penalty term Qξ=Pf·P- as an additional source term with Pf in the order of 106 [1] ensures the Reynolds cavitation boundary conditions [13] and represents the cavitation model according to Wu [12]. Zero pressure on the boundaries of ΩP is ensured by aDirichlet boundary condition. The coefficients ε-, ρ-HH, and ρ-Ht- represent the first, second, and third term of (5) in dimensionless form:(21)ε-=H3λ-∫01ρ-∫0ZZ~η-dZ~-∫0Z1η-dZ~∫01Z/η-dZ∫011/η-dZdZwithλ-=Rx2·ηTM·vΣ2·pH·bH3,ρ-HH=2HvΣv2-v1∫011/η-dZ∫01ρ-∫0Z1η-dZ~dZ+∫01ρ-v1dZ,ρ-Ht-=H∫01ρ-dZ.Note thereby that all the parts of the coefficients in (21) and of the stabilisation terms in (20) can, despite the film thickness H, be prepared in the sequence controller and transferred to the FEM-model (P,H) byinterpolation functions.The mesh in the FEM-model (P,H) isfree triangular with a normal element size and a refinement at ΩP by adistribution with typically 1000 elements.


### 3.3.2. FEM-Model (T-)

Figure3 (top right) shows the computational domains ΩT, ΩT,1, and ΩT,2 on which the temperature distribution T- in lubricant and solid bodies is calculated fully coupled. Figure 3 (bottom right) illustrates in parts the corresponding model tree. The lubricant domain ΩT has the dimensionless height of 1 with Z∈0,1, whereas the solid domains ΩT,1 and ΩT,2 have the dimensionless heights of D=3.15 with Z∈-D,0 and Z∈1,D+1, respectively. The length of the computational domains extends from Xin to Xex (see Figure 3, top left). For each new call-up of the FEM-model (T-), the solution of the previous call-up is used as initial value. The transient energy equations in dimensionless form are implemented and calculated in thephysics Heat Transfer in Fluids for the lubricant domain ΩT and inHeat Transfer in Solids for the solid domains ΩT,1 and ΩT,2. The dependent variable T- is solved based on Lagrange quadratic elements. The numerical stabilisation terms ID, SUPG, and GLS are already implemented in the consideredheat transfer physics in COMSOL and applied when required. The velocity distribution (4), viscosity distribution (15), density distribution (16), thermal conductivity distribution (17), and heat capacity distribution (18) in ΩT are calculated in the sequence controller and transferred to the FEM-model (T-) byinterpolation functions. Furthermore, theheat sources in thephysics Heat Transfer in Fluids, due to the shearing and compression of the lubricant, are calculated in the sequence controller and transferred to FEM-model (T-). Note thereby that, similar to the coefficients in (21), T- as a dependent variable needs to be applied in the FEM-model (T-). The implementation of the energy equations on ΩT,1 and ΩT,2 in theHeat Transfer in Solids modules with the translation motions of solids v1 and v2 as input parameters is straightforward.The boundary conditions in Section2.4 are implemented by defining bulk temperature T-M=1 at the left and upper boundary of ΩT,2 and at the left and lower boundary of ΩT,1, anOutflow boundary condition at the right-hand boundaries of ΩT,1 and ΩT,2 and anOpen boundary condition at the left- and right-hand boundary of the lubricant domain ΩT with lubricant entering ΩT having a bulk temperature T-M=1. As a common coordinate system in the FEM-model (T-) is used, manual adjustments are necessary to respect the different dimensionless definitions of Z in ΩT and ΩT,1 and ΩT,2, respectively (19). To be precise, the dimensionless scale of the Z axis of the solid domains ΩT,1 and ΩT,2 has to be adjusted according to the dimensionless scale of the Z axis of the lubricant domain ΩT, so that the continuity condition of the conductive heat flux and temperature across the lubricant-solid boundaries is respected. This is handled in COMSOL by multiplying the dimensionless energy equations of the solids by the conversion factor Zcf. This factor is the ratio of the dimensionless scales of the Z axes of the lubricant domain and the solid domains:(22)Zcf=ZlubricantZsolid=z/hz/bH=RxH·bH.The mesh in the FEM-model (T-) is divided into the lubricant domain ΩT and the solid domains ΩT,1 and ΩT,2. Whereas ΩT is meshed by amapped equidistant mesh, ΩT,1 and ΩT,2 are meshed withfree triangularelements with a refinement on the solid-lubricant boundaries and anextra-fine element size. The lubricant domain ΩT may consist of 40 elements in height and 1000 elements in the length direction.


### 3.3.3. Sequence Controller

The communication between the MATLAB sequence controller and the FEM-model (P,H) and FEM-model (T-) is achieved by COMSOL’sLiveLink for MATLAB function. The MATLAB sequence controller basically represents the implementation of the global numerical scheme shown in Figure 2. There, it is very convenient, for example, to implement the models for the lubricant properties, calculate the integrals by adaptive Simpson quadrature rule with typically 40 elements in Z-direction, and perform other calculations with the current values for pressure, film thickness, and temperature. As no other possibility is available, the data exchange of vectors and matrices from MATLAB to COMSOL is handled in a somewhat cumbersome manner: The quantities, calculated in the sequence controller, for example, velocity distribution, are written into a text file and read by the COMSOLinterpolation functions. For data exchange from COMSOL to MATLAB, thempheval command is used.


### 3.3.4. Solver

For solving the nonlinear system of equations, two separate COMSOL solvers for the FEM-model (P,H) and for the FEM-model (T-) are required. The solver for the FEM-model (P,H) includes the dependent solution variables P, H0, and U-, whereas the solver for the FEM-model (T-) includes the temperature T- as the dependent solution variable. For both FEM-models, a direct COMSOLMUMPS (Multifrontal Massively Parallel Direct Solver) [25] is used with a damped Newton-Raphson approach [26] and 10-3 as the tolerance factor as defined in [7].


## 3.3.1. FEM-Model (P,H)

Figure3 (top left) shows the computational domains Ωδ and ΩP on which the pressure and film thickness distribution are calculated fully coupled. Figure 3 (bottom left) illustrates in parts the corresponding model tree. For each new call-up of FEM-model (P,H), the solution of the previous call is used as an initial value.The elastic deformation of the equivalent body inZ-direction δ-(X,t-) is calculated on a quadratic two-dimensional computational domain Ωδ with the extensions X∈-30,30 and Z∈-60,0. Hence, an influence of the boundary conditions of Ωδ on the contact area can be excluded. Thephysics Solid Mechanics is used with the displacement vector U-(X,t-) as the dependent variable. Alinear elastic isotropic material specified by Young’s Modulus Eeq and Poisson’s ratio νeq is assumed. The boundary conditions described in Section 2.3 are implemented as afixed constraint condition on the bottom boundary of Ωδ, aboundary load condition with theload type pressureP on ΩP, and afree boundary condition on the remaining boundaries of Ωδ. Lagrange quadratic elements are used for the linear elasticity equations.The pressure and film thickness distributionsP and H, as well as the constant parameter H0 based on the load balance in (8), are calculated on the one-dimensional computational domain ΩP with an extension from Xin=-4.5 to Xex=3. Thereby, Xex of ΩP has already been enlarged for simulations with varying load over time, for example, [5]. The dependent variables are the pressure P, the constant parameter H0, and the displacement vector U-. The film thickness in (7) is implemented as avariable and the deviation from the smooth profile R is read as the input parameter by aninterpolation function. The load balance in (8) is implemented as aGlobal Equation with acomponent coupling integration on ΩP for the left-hand side of the load balance equation. Its result is the constant parameter H0. As no specific module for the generalised Reynolds equation in (5) is available in COMSOL, it is defined in thephysics Weak Form Boundary PDE by using Lagrange quintic elements. Based on the dimensionless parameters in (19), the weak form of the generalised Reynolds equation in dimensionless form is written as(20)-∫ΩPε-+ID∂P∂X∂WP∂XdX-∫ΩPρ-HH∂WP∂XdX+∫ΩP∂ρ-Ht-∂t-WPdΩP+∫ΩPQξWPdΩP+SUPGGLS=0withSUPGGLS=SUPGorGLS.WP is the test function for the pressure P and is denoted as test(P) in COMSOL. ID and SUPG/GLS represent numerical stabilisation terms, which are implemented asvariables and given here for completeness. Note that, for moderate cases, the TEHL simulation model can work without stabilisation terms. The numerical stabilisation terms have been applied to the generalised Reynolds equation by Habchi [1, 24] in order to avoid the oscillatory behaviour of the Reynolds equation at high loads. SUPG (Streamline-Upwind-Petrov-Galerkin) andGLS (Galerkin-Least-Squares) are consistent stabilisation terms whereasID (Isotropic Diffusion) is a nonconsistent stabilisation term. For implementation of ID and SUPG/GLS, the reader is referred to [1, 24]. The penalty term Qξ=Pf·P- as an additional source term with Pf in the order of 106 [1] ensures the Reynolds cavitation boundary conditions [13] and represents the cavitation model according to Wu [12]. Zero pressure on the boundaries of ΩP is ensured by aDirichlet boundary condition. The coefficients ε-, ρ-HH, and ρ-Ht- represent the first, second, and third term of (5) in dimensionless form:(21)ε-=H3λ-∫01ρ-∫0ZZ~η-dZ~-∫0Z1η-dZ~∫01Z/η-dZ∫011/η-dZdZwithλ-=Rx2·ηTM·vΣ2·pH·bH3,ρ-HH=2HvΣv2-v1∫011/η-dZ∫01ρ-∫0Z1η-dZ~dZ+∫01ρ-v1dZ,ρ-Ht-=H∫01ρ-dZ.Note thereby that all the parts of the coefficients in (21) and of the stabilisation terms in (20) can, despite the film thickness H, be prepared in the sequence controller and transferred to the FEM-model (P,H) byinterpolation functions.The mesh in the FEM-model (P,H) isfree triangular with a normal element size and a refinement at ΩP by adistribution with typically 1000 elements.


## 3.3.2. FEM-Model (T-)

Figure3 (top right) shows the computational domains ΩT, ΩT,1, and ΩT,2 on which the temperature distribution T- in lubricant and solid bodies is calculated fully coupled. Figure 3 (bottom right) illustrates in parts the corresponding model tree. The lubricant domain ΩT has the dimensionless height of 1 with Z∈0,1, whereas the solid domains ΩT,1 and ΩT,2 have the dimensionless heights of D=3.15 with Z∈-D,0 and Z∈1,D+1, respectively. The length of the computational domains extends from Xin to Xex (see Figure 3, top left). For each new call-up of the FEM-model (T-), the solution of the previous call-up is used as initial value. The transient energy equations in dimensionless form are implemented and calculated in thephysics Heat Transfer in Fluids for the lubricant domain ΩT and inHeat Transfer in Solids for the solid domains ΩT,1 and ΩT,2. The dependent variable T- is solved based on Lagrange quadratic elements. The numerical stabilisation terms ID, SUPG, and GLS are already implemented in the consideredheat transfer physics in COMSOL and applied when required. The velocity distribution (4), viscosity distribution (15), density distribution (16), thermal conductivity distribution (17), and heat capacity distribution (18) in ΩT are calculated in the sequence controller and transferred to the FEM-model (T-) byinterpolation functions. Furthermore, theheat sources in thephysics Heat Transfer in Fluids, due to the shearing and compression of the lubricant, are calculated in the sequence controller and transferred to FEM-model (T-). Note thereby that, similar to the coefficients in (21), T- as a dependent variable needs to be applied in the FEM-model (T-). The implementation of the energy equations on ΩT,1 and ΩT,2 in theHeat Transfer in Solids modules with the translation motions of solids v1 and v2 as input parameters is straightforward.The boundary conditions in Section2.4 are implemented by defining bulk temperature T-M=1 at the left and upper boundary of ΩT,2 and at the left and lower boundary of ΩT,1, anOutflow boundary condition at the right-hand boundaries of ΩT,1 and ΩT,2 and anOpen boundary condition at the left- and right-hand boundary of the lubricant domain ΩT with lubricant entering ΩT having a bulk temperature T-M=1. As a common coordinate system in the FEM-model (T-) is used, manual adjustments are necessary to respect the different dimensionless definitions of Z in ΩT and ΩT,1 and ΩT,2, respectively (19). To be precise, the dimensionless scale of the Z axis of the solid domains ΩT,1 and ΩT,2 has to be adjusted according to the dimensionless scale of the Z axis of the lubricant domain ΩT, so that the continuity condition of the conductive heat flux and temperature across the lubricant-solid boundaries is respected. This is handled in COMSOL by multiplying the dimensionless energy equations of the solids by the conversion factor Zcf. This factor is the ratio of the dimensionless scales of the Z axes of the lubricant domain and the solid domains:(22)Zcf=ZlubricantZsolid=z/hz/bH=RxH·bH.The mesh in the FEM-model (T-) is divided into the lubricant domain ΩT and the solid domains ΩT,1 and ΩT,2. Whereas ΩT is meshed by amapped equidistant mesh, ΩT,1 and ΩT,2 are meshed withfree triangularelements with a refinement on the solid-lubricant boundaries and anextra-fine element size. The lubricant domain ΩT may consist of 40 elements in height and 1000 elements in the length direction.


## 3.3.3. Sequence Controller

The communication between the MATLAB sequence controller and the FEM-model (P,H) and FEM-model (T-) is achieved by COMSOL’sLiveLink for MATLAB function. The MATLAB sequence controller basically represents the implementation of the global numerical scheme shown in Figure 2. There, it is very convenient, for example, to implement the models for the lubricant properties, calculate the integrals by adaptive Simpson quadrature rule with typically 40 elements in Z-direction, and perform other calculations with the current values for pressure, film thickness, and temperature. As no other possibility is available, the data exchange of vectors and matrices from MATLAB to COMSOL is handled in a somewhat cumbersome manner: The quantities, calculated in the sequence controller, for example, velocity distribution, are written into a text file and read by the COMSOLinterpolation functions. For data exchange from COMSOL to MATLAB, thempheval command is used.


## 3.3.4. Solver

For solving the nonlinear system of equations, two separate COMSOL solvers for the FEM-model (P,H) and for the FEM-model (T-) are required. The solver for the FEM-model (P,H) includes the dependent solution variables P, H0, and U-, whereas the solver for the FEM-model (T-) includes the temperature T- as the dependent solution variable. For both FEM-models, a direct COMSOLMUMPS (Multifrontal Massively Parallel Direct Solver) [25] is used with a damped Newton-Raphson approach [26] and 10-3 as the tolerance factor as defined in [7].


## 4. Model Validation

In order to show the functionality and the plausibility of the TEHL simulation model described in Section3, its solutions are compared with steady-state solutions from Hartinger et al. [2] and with transient solutions presented by Venner [4].


### 4.1. Steady-State Validation

Hartinger et al. [2] investigated the TEHL contact by CFD simulations for various conditions. They show that the assumptions of the generalised Reynolds equation are reasonable and provide representative results for thermal and isothermal contacts. For comparison of their results in [2] with the presented TEHL simulation model, the same lubricant property models and properties for lubricant and solids are used. Note that cavitation effects considered by the CFD approach cannot be depicted by the cavitation model of Wu [12]. Figure 4 shows the calculated pressure and film thickness distributions for a thermal (TEHL) and isothermal (EHL) contact with high sliding proportion (vΣ=5m/s, vg=2.5m/s) in comparison to the presented TEHL simulation model (a) and CFD model of Hartinger et al. [2] (b). The results show very good agreement. Furthermore, the local distributions of shear rate, temperature, and viscosity (not shown) correlate very well with the results of Hartinger et al. [2]. The calculation time for the TEHL solution in Figure 4 is approximately 6 to 8 minutes on a computer with a 3.5 GHz processor for the presented TEHL simulation model and in the order of 24 hours for the CFD model [2].Figure 4
Comparison of simulated pressure and film thickness distributions for a thermal (TEHL) and an isothermal (EHL) contact between the TEHL simulation model of this study (a) and the CFD model of Hartinger et al. [2] (b).
(a)
(b)


### 4.2. Transient Validation

Venner [4] investigated an EHL contact with a moving surface indentation under high sliding conditions by using a transient isothermal Newtonian approach. For the sake of comparison, the authors also assume isothermal and Newtonian simplifications in their presented TEHL simulation model. Input parameters and the investigated surface indentation are thereby adopted from [4]. Figure 5 shows the comparison of pressure and film thickness distribution for the example time, when the surface indentation is located at the dimensionless position Xd=-0.25. It is thereby shown that the moving surface indentation initiates a significant film thickness change, which hurries ahead the actual position of the indentation. This can occur under high sliding conditions. The comparison in Figure 5 shows very good agreement. Small differences in the results may be caused by different elastic deformation calculations (FEM in the TEHL simulation model in this study and elastic half-space assumption in [4]).Figure 5
Comparison of simulated pressure and film thickness distribution for a transient isothermal Newtonian EHL contact between the TEHL simulation model of this study (a) and the solution of Venner [4] (b).
(a)
(b)The validation of the TEHL simulation model shows its functionality and plausibility.


## 4.1. Steady-State Validation

Hartinger et al. [2] investigated the TEHL contact by CFD simulations for various conditions. They show that the assumptions of the generalised Reynolds equation are reasonable and provide representative results for thermal and isothermal contacts. For comparison of their results in [2] with the presented TEHL simulation model, the same lubricant property models and properties for lubricant and solids are used. Note that cavitation effects considered by the CFD approach cannot be depicted by the cavitation model of Wu [12]. Figure 4 shows the calculated pressure and film thickness distributions for a thermal (TEHL) and isothermal (EHL) contact with high sliding proportion (vΣ=5m/s, vg=2.5m/s) in comparison to the presented TEHL simulation model (a) and CFD model of Hartinger et al. [2] (b). The results show very good agreement. Furthermore, the local distributions of shear rate, temperature, and viscosity (not shown) correlate very well with the results of Hartinger et al. [2]. The calculation time for the TEHL solution in Figure 4 is approximately 6 to 8 minutes on a computer with a 3.5 GHz processor for the presented TEHL simulation model and in the order of 24 hours for the CFD model [2].Figure 4
Comparison of simulated pressure and film thickness distributions for a thermal (TEHL) and an isothermal (EHL) contact between the TEHL simulation model of this study (a) and the CFD model of Hartinger et al. [2] (b).
(a)
(b)


## 4.2. Transient Validation

Venner [4] investigated an EHL contact with a moving surface indentation under high sliding conditions by using a transient isothermal Newtonian approach. For the sake of comparison, the authors also assume isothermal and Newtonian simplifications in their presented TEHL simulation model. Input parameters and the investigated surface indentation are thereby adopted from [4]. Figure 5 shows the comparison of pressure and film thickness distribution for the example time, when the surface indentation is located at the dimensionless position Xd=-0.25. It is thereby shown that the moving surface indentation initiates a significant film thickness change, which hurries ahead the actual position of the indentation. This can occur under high sliding conditions. The comparison in Figure 5 shows very good agreement. Small differences in the results may be caused by different elastic deformation calculations (FEM in the TEHL simulation model in this study and elastic half-space assumption in [4]).Figure 5
Comparison of simulated pressure and film thickness distribution for a transient isothermal Newtonian EHL contact between the TEHL simulation model of this study (a) and the solution of Venner [4] (b).
(a)
(b)The validation of the TEHL simulation model shows its functionality and plausibility.


## 5. Example Results

In this section, examples of possible application opportunities, such as simulating coated surfaces, coefficients of friction, and nonsmooth surfaces, of the described TEHL simulation model are demonstrated.


### 5.1. Thin Surface Coatings

Firstly, the possibility of investigating the effect of thin surface coatings on the TEHL contact behaviour is shown. Results based on the TEHL model described in this study with slightly different lubricant properties have already been published by Lohner et al. [27]. There, the influence of thin surface coatings on the steady-state TEHL contact temperature was analysed by adding coating layers with separate thermophysical properties to the FEM-model (T-). The FEM-model (P,H) is unchanged. In the FEM-model (T-), additional geometries with thephysics Heat Transfer in Solids including adopted boundary and continuity conditions as described in Section 3.3.2, additional interpolation functions in COMSOL, and additional meshes for the thin surface coatings are integrated. The corresponding handling in the sequence controller is the same as for solids ΩT,1 and ΩT,2. This is very convenient to implement and shows the flexibility of the presented TEHL model approach. For simplification, it has been assumed that the deformation of the solid bodies is not affected by the very thin coatings (1.5 μm) [6].Figure6 shows an example of the pressure and film thickness (a) and the corresponding temperature distribution of the lubricant and solid bodies (b) for a TEHL contact including a coating with a thickness of 1.5 µm on the upper solid body. The bulk material is case-carburised steel 16MnCr5 and the coating material SiO2 with a significant different thermal conductivity (λ(16MnCr5) = 44 W/m/K, λ(SiO2) = 1.4 W/m/K). The SiO2 coating shows an insulating effect and leads to significantly higher temperatures compared to an uncoated TEHL contact. This also affects the coefficient of friction as shown by, for example, Habchi [6] for TEHL point contacts.Figure 6
Pressure and film thickness (a) and temperature distribution (b) of a TEHL contact with coated (SiO2) upper solid acc. to Lohner et al. [27].
(a)
(b)


### 5.2. Coefficient of Friction

For comparing measurement and simulation results, the coefficient of friction is a suitable quantity with integral character. Coefficient of friction measurements in fluid film lubrication regime can be performed at the FZG twin disk test rig with, for example, cylindrical discs (r1=r2=40mm, leff=5mm) made from case-carburised steel (16MnCr5) with a polished surface structure. Further information about the test rig is documented in [28], for example. Table 1 shows the main properties of the considered lubricant MIN100 and the physical properties of the bulk material. Model constants of Section 2.5 of the considered lubricant have been obtained from regression of measurement data. The Eyring shear stress of τc=6N/mm2 has been adopted from [29] as there a mineral oil of the same viscosity has been considered. Note that based on the Eyring shear model the increase of the shear stress with shear rate is not limited. To achieve a more realistic representation of the rheological behaviour of the lubricant, a shear stress model based on a limiting shear stress is currently being developed by the authors. However, this is not the focus of this study as it can be easy implemented once it is ready.Table 1
Lubricant and bulk material properties (i=1,2).
Bulk material properties
16MnCr5
E
i in N/mm2
206000
ν
i
0.30
ρ
i in kg/m3
7760
λ
i in W/m/K
44
c
p
,
i in J/kg/K
431
Lubricant properties
MIN100
υ (40°C) in mm2/s
95
υ (100°C) in mm2/s
10.5
ρ (15°C) in kg/m3
885
λ
0 in W/m/K
0.137
c
p (15°C) in J/kg/K
1921
τ
c in N/mm2
6A comparison of measured and simulated coefficients of the friction of lubricant MIN100 with very good correlation is shown in Figure7. The considered sliding velocities vg={0,0.08,0.21,0.42,0.64,0.89,1.41,2.00,2.67}m/s are thereby accompanied with corresponding bulk temperatures of the slower disc of ϑM=42,47,56,69,80,90,109,123,140°C.Figure 7
Simulated and measured coefficients of frictionμ over the sliding velocity vg.


### 5.3. Nonsmooth Surfaces

Surface features have a significant influence on TEHL contact behaviour [30, 31]. As an example, the contact of a deterministically structured surface (upper body) with repetitive protruding bars with a height of 0.2 µm and a length of 50 µm every 30 µm shown in Figure 8 with a smooth surface (lower body) is investigated. The deterministically structured surface is represented by the term R in (7).Figure 8
Schematic representation of the deterministically structured surface.The considered lubricant and material parameters are shown in Table1 and the operation conditions are pH=1200N/mm2, vg=0.89m/s, vΣ=8m/s, and ϑM=89.5°C. The model constants of the considered lubricant MIN100, as well as the macrogeometry of the discs, are the same as in Section 5.2. Figure 9 shows the simulated pressure, film thickness, and temperature distribution, which are significantly influenced by the deterministically structured surface. Very local changes of the pressure, film thickness, and temperature distributions are observed thereby. The maximum temperature rise in the lubricant is simulated to be 4 K higher for the deterministically structured surface than for a TEHL contact with smooth surfaces. The transient calculation time is approximately 5 to 10 hours on a computer with a 3.5 GHz processor. Postprocessing the simulated pressure and shear stress distribution by applying it to the upper and lower solid body with bulk material properties acc. to Table 1 easily delivers the stress distributions in the solid bodies based on linear elastic isotropic material behaviour. Figure 10 shows an example of the von Mises stress distribution in the lower (a) and upper (b) solid body corresponding to the pressure and shear stress distribution in Figure 9. Stresses due to temperature are omitted. Due to the deterministically structured surface, stress maxima close to the surface occur in addition to the well-known stress maxima in the deeper material depth.Figure 9
Temperature (a) and pressure and film thickness distribution (b) of a transient TEHL contact with a deterministically structured surface shown in Figure8.
(a)
(b)Figure 10
von Mises stress in the lower solid body (a) and in the upper solid body (b) correlating to Figure9.
(a)
(b)


## 5.1. Thin Surface Coatings

Firstly, the possibility of investigating the effect of thin surface coatings on the TEHL contact behaviour is shown. Results based on the TEHL model described in this study with slightly different lubricant properties have already been published by Lohner et al. [27]. There, the influence of thin surface coatings on the steady-state TEHL contact temperature was analysed by adding coating layers with separate thermophysical properties to the FEM-model (T-). The FEM-model (P,H) is unchanged. In the FEM-model (T-), additional geometries with thephysics Heat Transfer in Solids including adopted boundary and continuity conditions as described in Section 3.3.2, additional interpolation functions in COMSOL, and additional meshes for the thin surface coatings are integrated. The corresponding handling in the sequence controller is the same as for solids ΩT,1 and ΩT,2. This is very convenient to implement and shows the flexibility of the presented TEHL model approach. For simplification, it has been assumed that the deformation of the solid bodies is not affected by the very thin coatings (1.5 μm) [6].Figure6 shows an example of the pressure and film thickness (a) and the corresponding temperature distribution of the lubricant and solid bodies (b) for a TEHL contact including a coating with a thickness of 1.5 µm on the upper solid body. The bulk material is case-carburised steel 16MnCr5 and the coating material SiO2 with a significant different thermal conductivity (λ(16MnCr5) = 44 W/m/K, λ(SiO2) = 1.4 W/m/K). The SiO2 coating shows an insulating effect and leads to significantly higher temperatures compared to an uncoated TEHL contact. This also affects the coefficient of friction as shown by, for example, Habchi [6] for TEHL point contacts.Figure 6
Pressure and film thickness (a) and temperature distribution (b) of a TEHL contact with coated (SiO2) upper solid acc. to Lohner et al. [27].
(a)
(b)


## 5.2. Coefficient of Friction

For comparing measurement and simulation results, the coefficient of friction is a suitable quantity with integral character. Coefficient of friction measurements in fluid film lubrication regime can be performed at the FZG twin disk test rig with, for example, cylindrical discs (r1=r2=40mm, leff=5mm) made from case-carburised steel (16MnCr5) with a polished surface structure. Further information about the test rig is documented in [28], for example. Table 1 shows the main properties of the considered lubricant MIN100 and the physical properties of the bulk material. Model constants of Section 2.5 of the considered lubricant have been obtained from regression of measurement data. The Eyring shear stress of τc=6N/mm2 has been adopted from [29] as there a mineral oil of the same viscosity has been considered. Note that based on the Eyring shear model the increase of the shear stress with shear rate is not limited. To achieve a more realistic representation of the rheological behaviour of the lubricant, a shear stress model based on a limiting shear stress is currently being developed by the authors. However, this is not the focus of this study as it can be easy implemented once it is ready.Table 1
Lubricant and bulk material properties (i=1,2).
Bulk material properties
16MnCr5
E
i in N/mm2
206000
ν
i
0.30
ρ
i in kg/m3
7760
λ
i in W/m/K
44
c
p
,
i in J/kg/K
431
Lubricant properties
MIN100
υ (40°C) in mm2/s
95
υ (100°C) in mm2/s
10.5
ρ (15°C) in kg/m3
885
λ
0 in W/m/K
0.137
c
p (15°C) in J/kg/K
1921
τ
c in N/mm2
6A comparison of measured and simulated coefficients of the friction of lubricant MIN100 with very good correlation is shown in Figure7. The considered sliding velocities vg={0,0.08,0.21,0.42,0.64,0.89,1.41,2.00,2.67}m/s are thereby accompanied with corresponding bulk temperatures of the slower disc of ϑM=42,47,56,69,80,90,109,123,140°C.Figure 7
Simulated and measured coefficients of frictionμ over the sliding velocity vg.


## 5.3. Nonsmooth Surfaces

Surface features have a significant influence on TEHL contact behaviour [30, 31]. As an example, the contact of a deterministically structured surface (upper body) with repetitive protruding bars with a height of 0.2 µm and a length of 50 µm every 30 µm shown in Figure 8 with a smooth surface (lower body) is investigated. The deterministically structured surface is represented by the term R in (7).Figure 8
Schematic representation of the deterministically structured surface.The considered lubricant and material parameters are shown in Table1 and the operation conditions are pH=1200N/mm2, vg=0.89m/s, vΣ=8m/s, and ϑM=89.5°C. The model constants of the considered lubricant MIN100, as well as the macrogeometry of the discs, are the same as in Section 5.2. Figure 9 shows the simulated pressure, film thickness, and temperature distribution, which are significantly influenced by the deterministically structured surface. Very local changes of the pressure, film thickness, and temperature distributions are observed thereby. The maximum temperature rise in the lubricant is simulated to be 4 K higher for the deterministically structured surface than for a TEHL contact with smooth surfaces. The transient calculation time is approximately 5 to 10 hours on a computer with a 3.5 GHz processor. Postprocessing the simulated pressure and shear stress distribution by applying it to the upper and lower solid body with bulk material properties acc. to Table 1 easily delivers the stress distributions in the solid bodies based on linear elastic isotropic material behaviour. Figure 10 shows an example of the von Mises stress distribution in the lower (a) and upper (b) solid body corresponding to the pressure and shear stress distribution in Figure 9. Stresses due to temperature are omitted. Due to the deterministically structured surface, stress maxima close to the surface occur in addition to the well-known stress maxima in the deeper material depth.Figure 9
Temperature (a) and pressure and film thickness distribution (b) of a transient TEHL contact with a deterministically structured surface shown in Figure8.
(a)
(b)Figure 10
von Mises stress in the lower solid body (a) and in the upper solid body (b) correlating to Figure9.
(a)
(b)


## 6. Conclusion

In this paper, guidelines for the implementation of a TEHL model in commercial multiphysics software are presented. The introduced approach can be implemented with moderate effort, and the resulting TEHL model is very easy to extend for various applications. The TEHL model divided in a sequence controller in MathWorks MATLAB and in two self-developed decoupled COMSOL Multiphysics FEM-models is also capable of calculating the transient TEHL contact along the path of contact of gears. Even more challenging tasks, too, such as the extension of the described model to mixed lubrication regimes, are possible. By providing implementation guidelines for a TEHL simulation model in commercial multiphysics software, the authors are confident that the research in computational tribology has been stimulated and accelerated.



---

*Source: 6507203-2016-02-01.xml*