Transient conformal TEHL algorithms for multibody simulation

This article describes aspects of transient thermal elasto-hydrodynamical lubrication (TEHL) contact modelling for conformal contacts. This is to be utilized in a multibody simulation (MBS) framework for engineering purposes. The verification and proof of concept is done by implementation in the tool BEAST (Fritzson et al., 2014) and by comparision to published experiments and simulation results.


Introduction
Modelling and simulation of mechanical systems require often analysis of contacts between surfaces of bodies.This is a demanding task with respect to simulation performance versus accuracy.
Tribological contacts between surfaces of mechanical bodies can be divided into two groups; concentrated contacts and conformal contacts.The concentrated contacts are characterized by: • large space derivatives in gap/intersection, • very small extension compared with surface/body dimension, • very local deformations/stresses, and • high pressures.
while the conformal contacts are characterized by: • small space derivatives in gap/intersection, • extension similar to surface/body dimension, • truncation by the surface edges, • having global deformations/stresses, and • low/medium pressures.
The most common type of contact situation in machine design is a conformal contact.Examples are: press fits between shafts and rings, screw joint contacting surfaces, hydrodynamic bearings (journal bearings, tilting pad bearings, squeeze film dampers, etc.), hydrostatic bearings, bushings, brakes, spline couplings, contacts in cam follower units, see Figure 1, seals, etc. Rolling bearings, see Figure 2, can also have conformal contacts, e.g., cage contacts with inner/outer ring, guide ring contacts with inner/outer ring or cage.In some cases can rolling element cage pocket contacts and rolling element flange contacts be considered conformal.
The most typical examples of concentrated contacts are the rolling element race contacts in rolling bearings.Rolling element flange contacts or cage pocket contacts are mostly concentrated contacts.Other examples are ball/race contacts in ball screws, cam/roller contacts in cam followers, contacts in roller screws, gear mesh contacts, etc.
The scope here is limited to contact models suitable for multi-body system simulation (MBS).Some characteristics of MBS are: • mechanical domain, i.e., motion of system of bodies that interacts, • large motion of bodies with nonlinear interactions, Figure 2: A spherical roller bearing model.It has conformal contacts between cages (blue/yellow) and inner/outer ring, guide ring (cyan) with inner/outer ring or cages.Typical concentrated contacts are rollers with inner/outer race.
• rigid thermal bodies representation, but also simplified flexible and thermal bodies, • use of reduction techniques to reduce number of states for FEM flexible or FEM thermal bodies, • long transient simulations with high accuracy demands, and • use of advanced ODE/DAE system solvers, i.e., variable order, variable step, for efficient solution.
That requires that the contact modelling: • supports all requirements of MBS system solver (ODE/DAE), e.g., high order of continuity, varying time step, Jacobians, iterations, rejected steps, varying order, etc., • does not to reduce the overall efficiency of the system solver, and • does not introduce algebraic constraints to the system solver.
Well defined algebraic equations can be efficiently handled by a DAE solver.However, constraints caused by general contacts are not at all like that and should preferable be avoided.
The main problems to be solved for contact modelling are: • to be able to handle all situations needed for engineering, i.e., with respect to geometry, tribology, boundary conditions, damage, history, etc., • to have sufficient model accuracy for engineering purposes but no more, since more detailed models come at a cost, and • very good computation performance.
If we cannot compute in realistic simulation times, then the contact model will not be used.
The scope of this article is further limited to conformal contacts.The main ideas to be exploited/investigated here for conformal contacts are: • having parametric formulation for "any" geometry, • utilizing the transient simulation situation instead of suffering from it, • continous models over all lubrication domains, i.e., from dry to fully lubricated, and cavitation, • that the deformations are mainly handled by the structural flexible bodies and not the elastic contact model, • reducing numerical stiffness whenever possible, e.g., in friction models, for increased performance, • avoiding algebraic conditions/constraints, use the "soft" model approach instead, • embedded/integrated optimized solvers for the TEHL PDEs in each contact, and do not use system solver for that purpose, and • PDEs in semi-implicit formulation with sub-steps controlled by stability and accuracy criteria.
The verification and proof of concept is done by implementation in tool BEAST (Fritzson et al., 2014) and comparision to published experiments and simulation results, when available.

Basic derivation
The transient Reynolds equation comes from considering the mass balance where the volume flow rate per unit width in x and y direction are The mass flow per unit width in x and y direction are then ṁx = ρ • q x , ṁy = ρ • q y (4) where a indicates the upper surface, and b the lower surface.The variable U is surface velocity in xdirection, and V in y-direction.
It is assumed that the viscosity and density are constant across the film thickness at a specific location (x, y), but they are strongly dependent on pressure and temperature at that location.The fluid properties are evaluated for the local film mean temperature state T m , see Section 6.For fluid pressure dependency, see Section 3.
In our contact model the overall structural/rigid body relative motion is separated from superimposed local elastic deformations that are part of the contact model.Thus we have in the normal direction Applying this to a symmetric contact plane between surface a and surface b, gives Now, as a very simple model for local deformations we use a "brush" model, i.e., where L is a typical deformation length, and E 1 , E 2 are Youngs modulus of the two materials of surface a and b.This adds an elastic transient term to the model.In conformal contacts, the value L is choosen so that the structural contact surface deformations of the bodies are dominating over this local contact deformation term.
A second transient term comes from the time derivative of density in the Reynolds equations, i.e., compression of the lubricant.This term is neglected in quasi-static analysis.The term h ∂ρ ∂t can be written as h ∂ρ ∂p ∂p ∂t .Expanding the mass balance (1) derivatives, applying ( 10) and ( 13), extends Reynolds equation to a transient elastohydrodynamic equation that can be written as Thus we have obtained a transient elastic extension of the Reynolds equations.The smaller the constant C is, the stiffer the equation will be.

Introducing "elastic" pressure
Let us divide the pressure into two superimposed parts; "fluid" presssure p f , and "elastic" pressure p e such as where and gives thus the partial derivative The benefit of this representation in the PDEs, is that we can obtain higher numerical accuracy due to elimination of truncation of floating point numbers in subtractions.The p f can be very small compared to p e which is also computed fully implicit.This representation showed also improved stability properties.
Introducing this into Equation ( 15) gives for ∆ s ≥ 0 where, since p e ≡ 0 For ∆ s < 0 we utilize that into Equation ( 15) to obtain Note, that Equation ( 19) is derived for ∆ s ≥ 0, and Equation ( 24) is derived for ∆ s < 0.
The space derivatives can also be expanded, but the benefits are not obvious there.

Adding material damping
The internal damping of metallic materials is very low but can be important for numerical stability when no damping lubricant film is present.The material damping of steel and similar material has friction-like characteristics, i.e., it is proportional to the sign of the deformation rate and the deformation force.
Using ( 13) and adding a material damping constant C d << 1 we have since the sign function is a constant.The term sign( ∂∆e ∂t ) is a bit cumbersome to compute early in the numerical scheme, so some simplifications are in order.The damping is most important for dry conditions h = 0 which gives ∂∆e ∂t = − ∂∆s ∂t .If assuming constant gap and/or small relative sliding velocity we also have ∂∆s ∂t ≈ ∆W .Thus we have which means that the previous compliance C is now replaced by the expression ∆W ) when including material damping.
In the implementation, the sign change is implemented as "smooth" sign change in order to avoid numerical issues.

Parametric formulation in contact plane
The Equation ( 15) incorporates an elastic contact model, and speeds in a contact plane, i.e., in the middle of the lubricant film.Here the parametric formulation is incorporated into that equation.This so it can be generally applied to parametric surfaces.We now formulate a modified equation (based on Equation ( 15)) for an orthogonal curvilinear grid (u(x, y, z), v(x, y, z), n(x, y, z)) where the direction n defines a contact surface normal.
Assuming that the u and x directions are locally aligned, and directions v and y as well.
We also neglect high order effects in the relation of u and x, and v and y, i.e.

∂ ∂u
Numerical experiments on the parametric representations of interest has shown that this is is a valid assumption.
We then obtain where the velocities are ΣU = U a + U b , ΣV = V a + V b , ∆W = W a −W b , which are expressed in the parametric coordinate directions.

Discretisation
We discretise the (u, v) system with a regular grid with indeces (i, j) and uniform grid sizes δ u and δ v .Metric derivatives are assumed to be known at all grid points and we define the following difference operators As well as the averaging operators We define the upwind difference operators as where s u, frac and s v, frac are smooth switch functions (s frac ∈ [0, 1]) as function of f i,j in a small interval close to zero.This to avoid discontinous behaviour when switching between backward and forward difference.Following the discretisation above for Equation (32) we obtain Equation (33) where K ≡ ρh 3 12η , L U ≡ ρhΣU 2 , and L V ≡ ρhΣV 2 .Dividing by density, and expanding the discretisation operators, we obtain Equation (34).
3 Lubricant pressure dependence 3.1 Density for high/medium pressure Density models that cover high pressure are needed for EHL contact analysis.
There are many ways to define compressibility of the lubricant.It is a rather well investigated area.We selected the Dowson and Higginson expression (DH) (Dowson and Higginson, 1966) where p ≥ 0 is the pressure increase over atmospheric pressure and ρ 0 is the density at local temperature and atmospheric pressure.For p < 0 we set p = 0 in the formula.It is a simple expression that is computational efficient, and widely used.
We can compare these values of bulk modulus to some typical values for other media (at atmospheric pressure when relevant); K diamond = 443 • 10 9 Pa, K steel = 160 • 10 9 Pa, K water = 2.2 • 10 9 Pa, and K air = 1.4 • 10 5 Pa.
The Dowson Higginson expression reaches a constant density for very large pressure.One can question if this behaviour is physical, and if it can cause numeric issues in some situations since the bulk stiffness becomes infinite.

Density for low pressure -mixing air and oil
In conformal contacts we focus on the low pressure region with special attention on the very low pressures.The Dowson Higginson model does not cover this range.It needs to be addressed.
To have a density pressure model valid for very low pressures is a way to model cavitation.
The low pressure range is of importance in hydrostatic lubrication where also the influence of air en-trainment is an important factor for compressibility (Bassani, 1992).Typical values of dispersed air (air bubbles) in the lubricant are in the range 1.5% to 15% according to (Bassani, 1992).The air bulk modulus is very low so the entrained air affects the resulting bulk modulus.
The density of a lubricant with uniform dispersed air, where V l is the volume of pure lubricant (without air) and V a volume of air, is where the f v = Va V is the volume ratio of air.Derivating Equation ( 38) with respect to pressure and some rearrangements gives finally Thus we have the expression for effective compressibility C which is the inverse of the bulk modulus K.The ∂p is the bulk modulus of pure lubricant (no air bubbles), and K a = 1 ρa ∂ρa ∂p for air.This derivation of Equation ( 39) includes the effect that both volumes and densities are pressure dependent.
Using the mass ratio f m = ma m = f v • ρa ρ we can express Equation (38) as and the Equation (39) as The benefit of the mass ratio f m is that it is pressure independent.

Density of air
We have for an reversibel adiabatic process, the gas (air) law where p a is the absolute pressure.It gives the bulk modulus for air using the notation for specific heat ra- or expressed as compressibility and the relative volume change with pressure and the density For an isothermal process, then we have γ ≡ 1.0.At 40 • C and atmospheric pressure p a,0 = 0.09807 MPa we have γ = 1.40, and it increases slowly with pressure, e.g., γ = 1.51 for 3.923 MPa (Bassani, 1992).
Temperature has also a more direct effect on air density.The ideal gas law gives where p a is the absolute pressure [Pa], T is the absolute temperature [K], and R is the specific gas constant [J/(kg*K)] which for dry air is 287.058[J/(kg*K)].
We neglect the temperature influence for now since the purpose is just to provide a cavitation model where pressure is the important parameter and there are several gas producing/absorbing effects we neglect in this model anyway.For the same reason we simplify by setting γ ≡ 1.0.
Thus we end up with where −p vac is the (tuned) level of "vacum" pressure relative environment pressure we use in our model, p is the difference from environment pressure, and we set ρ a,0 = 1.

Viscosity
The Roelands formula (Roelands, 1966) for isothermal conditions is where p R,0 = 1.96 • 10 8 Pa and p is pressure increase over atmospheric/ambient pressure.This is valid for p ≥ 0 up to medium/high pressure.
We need to develop a relation valid for low pressure p < 0, i.e., in the cavitation zone.One purpose is to obtain correct viscous losses in the cavitation zone.
In the cavitation zone we have the formation of streamers and in between them gas/air.This indicates that we can use a linear combination of viscosity for lubricant and air.The viscosity dependency on relative density for the mixture reflects that a portion of the gap filled by lubricant respective air.Note, if we have a homogenous mix of lubricant and air, i.e., not cavitation, then it is difficult to say what the resulting viscosity will be.
It seems that air does not have much pressure dependency on its viscosity so we treat it as a constant.We select η a = 1.86 • 10 −5 Pa s at 25 • C.
where ρ is the effective density at the pressure, ρ l,0 the density of pure lubricant at zero pressure.The equation will only have effect for very low positive or negative pressures.This since air is very compressible and f density will soon reach unity.
We might limit f l to be within [0, 1], but the difference is small compared to the uncertainty/variation in high pressure viscosity formulas.

Elastic friction
Shear stress on the surface due to the friction will in reality build up some elastic contact shear deformation on the surfaces before gross slip occurs.A model for this "elastic friction" behaviour has been developed for fretting fatigue analysis (Fritzson et al., 2011).According to the model, the rate of surface shear stress at each point on the surface is where v is the relative sliding velocity, k is a shear stiffness, and τ S is the complete "Stribeck" shear stress without any initial elastic behaviour.Note, that τ S switches sign the same time as v.This equation is applied in each of the two parameter directions of the surface.The differential equations are implemented with implicit Euler, which is unconditional stable, and are solved over the whole contact surface in each Right Hand Side (RHS) call.
Compared to a smooth velocity based switching function for shear τ (v), this "elastic friction" has shown very good numerical properties for conformal contacts forces.We believe that the system solver sees a large reduction in numerical stiffness for the system.The "elastic friction" is also more physical correct model than the velocity based switching for friction shear.The downside is more data and more computation.We need to solve the differential equations and store all the states.
Note that a velocity based switching can be made less numerical stiff by reducing the slope ∂τ ∂v .However that makes it less physical accurate and should be tuned in each case study.

Requirements
So far we have assumed that we have lubricant available that completely fills the gap at the inlet of the contact, i.e., a fully lubricated contact.
However, in real situations we often have just a very thin layer (microns) of lubricant on each surface in front of the contact.Thus, for large parts of the contact, the gap is much larger than the film and the inlet miniscus will not fill that gap.Thus we have a "starved" contact situation.
To really investigate what is going on in a starved contact, one needs to do a 3-D (in the plane and through the gap), two phase (air and lubricant) CFD computation of the whole contact including inlet and outlet.There is ongoing research towards this ultimal goal (Bruyere et al., 2012;Tanaka et al., 2017), but the required analysis is not practical today, especially not for large scale investigations.
Here we aim to create an approximate simulation method for the whole lubrication regime, i.e., from fully lubricated gap to dry surfaces, based on the 2-D extended transient Reynolds equations.The requirement is that the method is suitable for engineering purposes.

Varying lubricant mix
Mostly we have not completely dry contacts or not completely filled contact gap, but instead something in between.
So far we have assumed that the lubricant used consists of a constant mixture (mass ratio) of air and pure lubricant, and we draw in lubricant from the contact boundaries with 100% of that mix.
If we instead let the mixture drawn in be determined by the conditions at the contact boundaries, we can model anything from "dry" (100% air) to 100% "wet" (gap completely filled with lubricant).The mixture, or mass fraction air, then becomes a state variable and will vary in the contact and in time.
The mass ratio f m = ma m = f v • ρa ρ , where the volume ratio f v = Va V , is initialized in the contact and outside using the supplied lubricant layers on surface 1 and 2 respectively (h nom,1 , h nom,2 ) and the initial film thickness h.When initializing, the pressure is zero where we have non-zero h.
The initial air content in the lubricant supply, Air Volume Fraction (AVF), is given at p = 0.This needs to be included in the computation.
The resulting volume ratio range is The f m is continuously updated outside the contact according to Equation 58.
How to compute the f m since we have temperature dependency on density?It affects the definition of AVF, i.e., it is given at a known reference temperature, or it is given for a local temperature in the contact.In practice the difference will be small since we treat ρ a,0 as constant, see Section 3.2, and the temperature effects on lubricant density are small.For now we choose the local temperature.
Looking at the mass balance for a discrete element, we arrive to where m is the mass per area unit, and ṁx , ṁy , mass flow per unit with, see Equation (4).
To advance this in time we can use a simple Euler first order method which can be semi-implicit if parts of the expression for ∂m ∂t are computed at the end of the step.
These equations can be applied to air mass m a and total mass m.However, we prefer to work with the state mass fraction f m , and then get Assume that if we at a certain point have a mass ratio, i.e., mix, then the mass flow in the two directions at the same point have the same ratio ṁa,x = f m ṁx , ṁa,y = f m ṁy . (63) Then Equation ( 62) becomes which is used for advance in time An upwind scheme with respect to mass flow is the most suitable since that describes better what is comming into the element.This is similar to the temperature integration.

Inlet conditions for starved contacts
In Section 5.2 we have shown a way to handle varying mix of lubricant and air assuming a resulting one phase flow in the Reynolds equations.
We might utilize this for a starved contact if complemented with an algorithm to select the position of the inlet or where we should start to integrate the contact.This is neccessary since the amount of lubricant and meniscus can be very small compared with the contact gap.There might also be no edges in the contact, e.g., periodic directions of rings, instead the gap just increases.
We add the condition f inside • h lub < ∆, where ∆ is the gap, to determine if the point is outside.This in addition to the edge truncation limits test conditions.To have a continuous border, we need to search the integration limits from outside.
Then these findings are expressed in a formula that fits into the computation framework in a robust way, see Section 5.5.

Verification method
Data is needed for verification.No 3-D two-phase flow CFD simulation results are yet available to us.However, for the concentrated contact analysis we have developed formulas for central film and rolling moment based on a large number of 1-D EHL computations.
The EHL method is originally based on (Venner, 1991) but later extended by us to displacement based input, larger parameter range, and starved calculations.In short the maximum available lubricant layer is given and the integration zone (meniscus position) is varied by a virtual pressure method until that flow condition is fullfilled (i.e., starved condition), or we reached a very large zone (which is in practice a fully flooded condition).

Film factor formula
Evaluating the effect of the parameter f inside , the value that gave thickest film was selected.In fact there is a maximum in film thickness with respect to f inside , first the film increases with increasing f inside , there is a maximum, and then the film starts to decrease.These maximal choice values, see Table 1, were also the values that gave best agreement with the EHL film formula results.
Table 1: Film factor f inside values that gave maximum central film thickness.This considering all three sum velocities [0.4,4.0, 40.0] m/s, and all three line load levels [300, 3000, 30000] N/m.When a range is given for f inside , the simulations in that range give similar results and a factor in the range can be used.

Verification results starvation
There are 216 test cases simulated of which 136 are in various degree of starvation and 80 have reached fully flooded condition.This according to the EHL formula starvation output parameter.
The results for comparision of central film thickness between transient Reynolds and EHL formula are given in Table 2.The following definitions are used in the table: where sub-script R stands for "Reference" which is the EHL formula.
In total for all the 216 cases, the central film relative error mean value is ∆E mean = −0.03431and the relative error standard deviation is ∆E rms = 0.2985 .
When comparing the results one should have the following in mind: • Both the EHL formulas and the transient Reynolds (without side flow) are approximate 1-D theories with different assumptions aimed for different application areas.
• Even if we try to keep the line load low in the studied cases to minimize the elastic deformation, one can discuss if the studied contact is conformal enough.
• When the EHL formulas from the EHL calculations were created, there were curve fitting errors involved.Also, the EHL calculation data used focus on highly loaded concentrated contacts.
• The generated EHL data for the EHL formula was for maximum 5 µm combined lubricant thickness as input.This means that our 10 µm and 100 µm cases are extrapolations.• The thermal film calculations are active for the transient Reynolds.
Considering these factors, the results for central film thickness are satisfactory and can be built upon.
In order to have a really fair verification, we need to have experimental data for starved conformal contacts, or accurate simulation results for starved conformal contacts, e.g., for two phase flow in 3-D EHL.This verification will be revisited when such data is available.
6 Thermal balance

Energy equation
The special form of the energy equation for determing change in temperature is according to (Bird et al., 2007)  This equation needs to be simplified so we can integrate it over the film thickness.The discussion of the viscous dissipation term is in Section 7.
Let us assume that conduction is kept in all directions even if it is largest across the film (z-direction), and convection dominates along the film (x-direction and y-direction).Then Equation (76) becomes Let us assume that the effects of volumetric expansion is small for now, which gives.
The conduction z-term is easy to integrate for any temperature profile.However, the convection terms can be difficult since we have combined velocity profile and temperature profile terms.Let us in the convective terms, and conductive x/y terms, replace the temperature with the mean temperature over the film thickness.Neglect spatial x/y derivatives of k.Thus we obtain which can easily be integrated over the film where q is the volumetric flow rate [m 2 /s] in each direction, i.e., q x = hu m and q y = hv m .The P loss is the integral of dP loss over the film thickness.
Let us assume a parabolic temperature profile, and express it in T m , ΣT = T a + T b , and ∆T = T a − T b , where T a is the upper (z = h/2) surface temperature, and T b is the lower (z = −h/2).It then becomes Using Equation (81) in Equation ( 80) gives where we can see that there is no dependency on the temperature difference ∆T .
The heat flows at the surfaces are which allow us to separate/superimpose the heat flow due to surface temperature difference (∆T ) from the heat flow due to losses in the film.
Rearranging Equation ( 82) we have If the film thickness is zero, we have and all the losses are conducted into the two surfaces.
The effective thermal conductivity, specially in the cavitation zone where we have a mix of gas and pure lubricant streamers, is determined similarly to the effective viscosity The effective specific heat coefficient c p is determined by the mass fractions of air and pure lubricant.However, the air mass fraction is very small ≈ 10 −4 so it can be neglected and we have c p = c p,l .

Numerical scheme
Assuming that the sum of surfaces temperature ΣT change rather slowly, we suggest the following "semiimplicit" integration scheme where ∆t = t 1 − t 0 is the time step, and k h = 12k h .The space derivatives and most of the coefficients are evaluated for T m,0 .The space derivatives of T m must be taken upwind of the mass flows (m x , m y ) to get correct boundary conditions.
Since the thermal equations of the bodies are considered to be slow and only evaluated after the completed timestep in an extra RHS, same scheme is adopted here.
The term including P loss should preferably not be evaluated more than once per step due to computational cost, however it is possible to compute the other parts during sub-steps without much extra costs in order to improve numerical stability.

Thermal contact conductance
We have conduction between the two surfaces if they have different temperatures.
In (ESDU 78029) a formula for thermal contact conductance is given where typical values for the constants are h ac0 = 440 W m 2 K and p ac0 = 0.39 • 10 6 Pa.These data are for steel/steel ring/housing contacts.It is for a dry contact where the asperity contact is where the heat is transfered, the air being a bad conductor.This is the background for many thermal contact conductance formula.
Then we need to add a term for the lubricant thermal contact conductance so it can be applied in HL or EHL contact situations where k is the effective thermal conductivity of the lubricant.The combined surface roughness σ represents the remaining "valleys" when the mean film thickness h goes to zero.Since there is no asperity contacts when we have a thick separating film etc., we utilize the fraction measure θ developed for boundary friction The thermal contact conductance can become quite high.With 3 GPa pressure and h ac0 = 10000 W m 2 K we have h ac = 8.2 • 10 6 W m 2 K .For a typical mineral oil, zero mean film thickness, and minimum roughness 1•10 −8 m we have h fc = 11.6•10 6W m 2 K .It can be that these values are unrealistic high and need to be limited.In reality there are other effects not included in the models, e.g., surface layers like oxides, etc.The highest value we found in literature so far is 200000 W m 2 K .Until better information is available, we use this value 200000 W m 2 K as an upper limit.

Power loss 7.1 Definitions
We have for laminar Newtonian flow in a lubricant film, see Section 2.1: Newton's law of viscosity is where the i and j can take values 1,2,3.The η is the shear viscosity and κ the dilatational viscosity.For an ideal gas we have that κ = 0, and for an incompressible fluid the second term is zero.

Local power loss
We consider only irreversible conversion of kinetic energy to thermal energy, i.e., viscous dissipation, for Newtonian fluids (Bird et al., 2007) Separating the terms with respect to viscosity type gives The dissipation function Φ v for Newtonian fluids with all terms expanded becomes It seems to be common practice (e.g.(Doki-Thonon, 2012;Glavatskikh, 2000)) to neglect all terms but two, i.e., ∂vx ∂z and ∂vy ∂z , where z is in the film thickness direction.Thus we have How large the error is for a somewhat compressible fluid in transient conditions due to the neglected terms, is an open question.
Utilizing the expressions for flow above, and integrating through the film [−h/2, h/2] gives the power loss rate density [W/m 2 ] (103)

Effects due to boundary friction and limiting shear stress
The power loss rate, Equation ( 103), is divided into two type of terms; the first related to shear flow, and the second pressure flow.The pressure flow losses will go to zero when the film thickness goes to zero, i.e., having boundary friction.Nothing much else can be done about this term.
The shear power loss term η h ∆U 2 can be replaced with a more general equivalent expression τ ab,x ∆U which leaves it open how the shear stress is computed.When we have a Newtonian film, it can be the ususal τ ab,x = η h ∆U , but when having boundary friction it is possible to use Coulomb friction τ ab,x = µ•p•sign(∆U ), or some other friction model.It will also be possible to include some non-Newtonian behaviour such as limiting the shear stress.
Thus we have the more general power loss rate density where the only requirement is that the shear stress has same sign as the relative velocity.

Effects due to elastic friction
The work to build up this elastic deformation should be excluded from the losses since it is reversible.
Excluding the "elastic" part of the contact shear/slip speed, the resulting speed is τ τ S v, and we obtain the power loss rate density where τ S is the complete "Stribeck" shear stress without any initial elastic behaviour.

Material hysterisis
Material hysteris is small but can be important for dry contacts which do not have the squeeze damping of the lubricant film.
A main mechanism for material damping in metals is the movement of dislocations in the stressed material volume during elastic deformation.The phenomenon has a very weak dependency on the magnitude of the strain rate, but is dependent on the sign of the strain rate.Thus, in a contact, a normal force opposing the relative movement is created.The magnitude of the damping force is proportional to the stressed volume.As a measure of that volume, the magnitude of the elastic contact force can be used.The most simple model that reflects the behaviour discussed above for the damping force is a dry squeeze damping coefficient, multiplied with the elastic contact force.The sign of the damping force should be chosen so that the damping force opposes the squeeze motion.The switch of the sign has to be done in a smooth way so it does not create numerical issues.The coefficients have been determined for different materials by means of bouncing ball experiments.
In the normal direction, the effect can be implemented by modifying the elastic compressibility of the contact surfaces C e .
where the damping coefficient D hyst includes the correct sign depending on motion.Then we obtain the power loss rate due to material hysteresis.
There is also material hysterisis in the tangential direction due to the elastic deformation in the elastic friction model.However, it is typical a magnitude lower due to the fact that friction coefficients are normally small.It can be interesting to investigate what effect it can have on the efficiency of the solver.It might be beneficial.
8 Numerical methods

Contact model with embedded specialized numerical methods
The PDEs of the contact analysis are very different from the ODE (or of the MBS system.Therefore we choose to employ tailored numerical methods for each special task.In this way we have optimal numerical methods for each type of equations.The system solver will not see these embedded states, the contact model is a black box to the system solver.
In order to make this strategy to succeed, one has to take care that the embedded numerical methods have high precision, have smooth solutions, can handle all the situations the system solver requires, and ensure that the combined system is numerical stable.
The main embedded numerical method is for the PDE for pressure, discussed below, but there are also embedded numerical methods for elastic friction, material removal, and thermal.

Backward Differential Formula (BDF) methods
To advance solution in time for the PDEs we employ BDF methods.
In the examples below we solve for the pressure p as function of time t.The sub-script 1 indicates current time, 0 previous time step, and −1 the one before that.The BDF1 or implicit Euler formula is The formula for BDF2 or implicit two-step BDF method has been derived several times independently.For c = 1 we have An alternative expression of (109) directly in time increments has less issues for δt 1 = 0.

Practical implementationssemi-implicit schemes
A main reason why we look at implicit schemes is that we have the relative motion (gap, velocities, etc.) at time t 1 which are the main input for the computation in the contact.
Introducing some more specific contact variables, the gap ∆, the film thickness h we have in principle the following scheme (example BDF1) for every RHS and i, j (omitted below, i.e. scalar expressions) If the film thickness prediction is negative h 1 < 0 then we compute a dry contribution only otherwise we compute the predicted pressure as As seen we have p 1 on the right hand side (RHS) which calls for iterations.If Newton iterations, then we need a Jacobian, etc.
An alternative way is to use a semi-implicit method (use p for previous sub-time step on right hand side) with sub-steps controlled by a stability criterion, see Section 8.4.
where subscript 1 refers to current sub-time step and subscript 0 to previous sub-time step.The ∆ s and C are at the current system time step and kept constant during the sub-stepping.
If the film thickness prediction is negative h 1 < 0 then we compute a dry contribution only otherwise we compute the predicted pressure as ∂p ∂t 1 = RHS(p 0 , h 1 ), ( 118) We have tried both first order (BDF1) and second order (BDF2) variant extensively, with some preference for the second order variant.
To conclude, the numerical method for pressure currently used is the second order BDF semi-implicit scheme using sub-stepping with a stability criterion and additional semi-empirical accuracy criteria.
At the end of a successful system time step, we need to update and save state p, which we have computed.

Stability and sub-step in time
Assuming that it is the Poiseuille term in Equation (32) that determines the stability limit, our first estimate is that the maximum allowed time step is (or is proportional) to δt stb which is defined as This is valid for each point, so the minimum value needs to be used.
When h goes to zero, the expression goes to infinity, so implementation has to take care of this case.
The stability criterion is evaluated at each sub-step since there can be significant changes during a system step.

Infinite journal bearingverification
The test case is a journal bearing where the shaft is displaced in the vertical direction.The radius of the shaft is 0.1 m.The radial clearance is 10 • 10 −6 m, and the radial displacement 5 • 10 −6 m.The rotation speed of the shaft is 200 rpm, or 20.944 rad/s.The lubricant temperature is 40 • C, at which the kinematic viscosity is 6.8 • 10 −5 m 2 /s.The density is 876 kg/m 3 .
The flow is allowed in v-direction only, i.e., we have the solution of an infinite journal bearing, even if the bearing shown is finite.
The grid is 17x670 where the larger number is in the circumferential direction.
The simulation is run for 10 revolutions where the solution has stabilized.
The full Sommerfeld stationary solution for incompressible lubricant gives a maximal pressure of 4.65 • 10 8 Pa, at location 2.30 rad.The pressure distribution is anti-symmetric.
It is not fully possible to represent this case in the transient formulation since we need a pressure time derivative, however we can get close by deactivating the elastic deformations and density change.The obtained solution is shown in Figure 6.As seen, we obtain the Sommerfeld solution with high degree of accuracy.The small differences seen can be accounted to the grid resolution, since data is recorded at discrete grid points only.
10 Finite immersed 360 • journal bearing -verification The test case is a 360 • journal bearing totally immersed in lubricant.This according to experiment 2 on page 66, with results shown in Figure 70.1 on page 70, in (Jakobsson and Floberg, 1957).The radius of the shaft is 0.05 m, the width is 0.1 m, the radial clearance is 1.82 • 10 −4 m.The load is F = 620 N. The aim was to test the bearing for nondimensional eccentricity = 0.6.The angular speed of the shaft is 48.6 rad/s.The viscosity is 0.0153 Pa/s.The density is assumed to be 876 kg/m 3 .
These values are compared in (Jakobsson and Floberg, 1957) with computed values for nondimensional eccentricity = 0.6 which gave a load of F = 634 N with direction 40 • .The power loss was not measured, but computed to 17.2 W.
The grid is 63x168 where the larger number is in the circumferential direction.The elastic contact deformation is activated, but has no significant effect at these low pressures.The simulation is run for 3 revolutions where the solution has stabilized.
Here, we simulated a power loss of 18.0 W, and a resulting load is F = 681 N, with direction 37.5 • .The pressure distribution including experimental data points are shown in Figure 9.
The agreement of our simulation model and the data given in (Jakobsson and Floberg, 1957) is fully satisfactory.11 Finite 360 • journal bearing with oil groove at 90 • before loadverification In (Floberg, 1959) a comprehensive study on journal bearing power loss has been conducted.Here we utilize those experiments to verify our simulation model.
The test bearing has a journal diameter of 99.968 mm, a non-dimensional clearance of 0.00462, and a bearing width of 100 mm.The lubricant density was assumed to be 900 kg/m 3 .
In the experiments, the load, angular speed, viscosity (through temperature), were varied, and the power loss was measured.See Table 12.1 on page 12 in (Floberg, 1959).
These experiments were simulated using our model, and the results presented in Figure 10.The transient simulation was for 10 revolutions, for dispersed air cavitation model, with Dowson-Higginson (DH) pure oil compressibility combined with gas compressibility (10% at p = 0), and elastic surfaces.
The agreement of our simulation model and the measured data is good.Probably, most of the differences are due to the fact that thermal effects are not yet included in the simulations.
The cases have also been simulated with the steady state model by J. Ståhl (Ståhl, 2002).The results agree well with the results from our model, see Figure 10.
Note, that the calculations by (Floberg, 1959) shown in Figure 10 are for a 180 • journal bearing which means that the power loss in 180 • (cavitation zone) is not included.This is why those values are on the low side.(Floberg, 1959).This for 360 • journal bearing with oil groove at 90 • before the load.
There is good agreement between the simulations and the result by Floberg (Floberg, 1969).
13 Tilting-pad thrust bearingthermal verification Verification of the thermal hydrodynamic model against experiments means that the test apparatus must be modelled with regard to thermal effects as well, and not just thermal effects in the film build up.This increases the complexity significantly, and thus the sources of errors.
Here we choose to utilize and follow the extensive experimental research done at Luleå University of Technology, Sweden, regarding investigation of a spherical pivoted hydrodynamic thrust bearing subjected to various conditions and lubricants (Glavatskikh, 2000;Glavatskikh and Larsson, 2000).

Experimental setup
The mechanical test rig, see Figure 14 and Figure 15, includes a housing, containing two identical thrust bearings to be tested and a main shaft, supported by two journal bearings.The principal problem in designing a thrust bearing test rig is the necessity of measuring a small value of friction while accommodating a substantial axial load on the bearing.This is achieved by means of a balanced pair of test bearings, and installing the whole bearing assembly on a low friction bearing.
The test bearing is an industrial "catalogue type" thrust bearing with steel-backed and babbit-faced tilting pads.The bearing characteristics are given in Table 3.The spherical pivot gives the ability to tilt not only in the direction of rotation but also in the radial direction.
Table 3: Test bearing characteristics (Glavatskikh, 2000;Glavatskikh and Larsson, 2000).however most of the analysis in (Glavatskikh, 2000) has been done for a synthetic oil ISO VG 46 of poly-αolefin type at 50 • C inlet temperature, therefore those experiments are choosen for the verification.The data for the lubricant is given in Table 4.

Simulation model
The BEAST model is shown in Figure 18.It has three bodies; a pad, a collar plus shaft, and a pivot support.The temperature of the lubricant flowing through the bearing is not known other at the inlet and outlet.We choose to set the environment (lubricant) temperature to the mean value of the temperatures at the inlet and outlet.
The heat transfer coefficients of the surfaces to the environment are not known either.We did a parameter Figure 17: Sensor location and their mounting in the collar and shaft (Glavatskikh, 2000;Glavatskikh and Larsson, 2000).
study changing those coefficients while mainly looking at the errors in power loss and temperatures.With one set of heat transfer coefficients for all test cases, we ended up choosing nearly the same set of coefficients as choosen in (Glavatskikh, 2000;Glavatskikh and Larsson, 2000), i.e., 1000 W m 2 K for the shaft and collar surfaces, and 500 W m 2 K for pad surfaces.Due to symmetry (there are two bearings) the left end surface of the shaft has zero heat transfer coefficient.
In the contact the heat transfer is between the contacting surfaces and the lubricant in the contact, thus the heat transfer to the environment is zero for those surfaces.
Since only one pad is modelled, i.e., we have 1/6 of the power loss, the heat transfer coefficients to the environment for the shaft/collar are reduced with the same factor.
The model is fully thermo-mechanical and transient.The thermo-mechanical effect is important since the film thickness is affected by the thermal deformations, especially for that of the pad.With accelerated thermal equations, 60 revolutions are needed to reach steady state conditions.

Results
The results are summarized in Table 5 and Table 6.
The following definitions are used where sub-script m stands for "measured" experimental data.The power losses are underestimated, but in reality we have viscous losses outside the contacts as well which are not included in the model.Using different heat transfer coefficients for the different speeds might

Figure 1 :
Figure 1: A cam follower unit model.It has nine conformal contacts and one concentrated contact.
. It is used to compute both with the standard contact model (using EHL formulas) and with the transient Reynolds equations.The 1-D flow behaviour is enabled by disabling the axial flow (u-direction).The discs are made of steel whereby assuming a Youngs modulus of 2.03 • 10 11 Pa, and Poisson ratio of 0.29.

Figure 3 :
Figure 3: The two disk model.Both discs rotate with same speed in different directions, i.e., perfect rolling conditions.One disc is loaded in the contact normal direction, the other degrees of freedom are prescribed or locked.The axial width is 10 mm.The simulations are run until steady state conditions occur.

Figure 4 :
Figure 4: The film factors which gave maximum film for D = 0.12 m.The developed function is also shown.

Figure 5 :
Figure 5: The film factors which gave maximum film for D = 0.24 m.The developed function is also shown.
lubricant c p is the specific heat coefficient [J/(kg K)], k is the thermal conductivity [W/(m K)], dP loss the irreversible viscous dissipation [W/m 3 ] (negative value, see Section 7), and β T = − 1 ρ ∂ρ ∂T p is the volumetric expansion coefficient [1/K].The equation above assumes Newtonian fluid and Fourier's law for conductivity.The Df Dt = ∂f ∂t + u ∂f ∂x + v ∂f ∂y + w ∂f ∂z is the substantial (or material) time derivative.

Figure 6 :
Figure 6: Pressure distribution after transient simulation of 10 revolutions, for incompressible lubricant and rigid surfaces.

Figure 7 :
Figure 7: Pressure distribution after transient simulation of 3 revolutions, for dispersed air cavitation model, with Dowson-Higginson (DH) pure oil compressibility combined with gas compressibility (10% at p = 0), and elastic surfaces.Markers outline the location z = 10 mm.

Figure 8 :
Figure 8: Pressure distribution after transient simulation of 3 revolutions, for dispersed air cavitation model, with Dowson-Higginson (DH) pure oil compressibility combined with gas compressibility (10% at p = 0), and elastic surfaces.Markers outline the location z = 30 mm.

Figure 9 :
Figure 9: Comparsion between experiments and simulation for z = 10 mm and z = 30 mm.
4 .The data in the tables are given as non-dimensional values.For the simulations below, the viscosity and density are taken from calculation examples in(Floberg, 1969), The speed and minimum separation are selected to give realistic pressure values for the application.Thus we have: η = 0.049 Pas, ρ = 900 kg m 3 , h min = 20 • 10 −6 m, and n = 100 rpm.The variables are defined as; n p number of pads/sectors, b 0 = b R non-dimensional pad width, k nondimensional height difference at mean radius, P 0 = F h min 2 ηU R non-dimensional pad load per unit width, E 0 = Eh min ηU 2 R non-dimensional power loss per unit width for one pad, x 0P = x P R non-dimensional load x-coordinate, y 0P = y P R non-dimensional load y-coordinate, where R is the mean radius, the speed U = R • ω, and F is the pad load per unit width.
sensors for temperature, film thickness, and pressure, see Figures16-17.The lubricant flow rate is 15 liters per minute in all tests, and the lubricant inlet and outlet temperature are measured as well.

Figure 18 :
Figure 18: The BEAST model of the tilting pad bearing.It consist of three parts; the shaft with collar, the pad, and the pivot support.The surface temperatures are shown.This for a case with 2 MPa load and 3000 rpm.

Table 2 :
Absolute and relative differences in central film thickness between transient Reynolds and EHL formula.This in each case for all three sum velocities [0.4,4.0,40.0]m/s, and all three line load levels[300, 3000, 30000] N/m.In total there are 216 test cases.The Shell type lubricant is at 40 • C.
, test case and results according to table 32.3 and table 34.