Model order reduction by projection applied to the universal Reynolds equation

Abstract

The approach presented in this paper yields a reduced order solution to the universal Reynolds equation for incompressible fluids, which is valid in lubrication as well as in cavitation regions, applied to oil-film lubricated journal bearings in internal combustion engines. The extent of cavitation region poses a free boundary condition to the problem and is determined by an iterative spatial evaluation of a superposed modal solution. Using a Condensed Galerkin and Petrov–Galerkin method, the number of degrees of freedom of the original grid is reduced to obtain a fast but still accurate short-term prediction of the solution. Based on the assumption that a detailed solution of a previous combustion cycle is available, a basis and an optimal test space for the Galerkin method is generated. The resulting reduced order model is efficiently exploited in a time-saving evaluation of the Jacobian matrix describing the elastohydrodynamic coupling in a multi-body dynamics simulation using flexible components. Finally, numerical results are presented for a single crankshaft main bearing of typical dimensions.

Keywords:

• Reynolds equation
• model order reduction
• Galerkin projection
• Jacobian matrix

1. Introduction

Pressure-oiled journal bearings are an essential machine part in various fields of engineering. The mathematical description of hydrodynamic pressure build-up in the oil film dates back to the end of nineteenth century when Osborne Reynolds derived the differential equation named after him [1]. It still remains valid and forms the basis of reams of bearing models proposed during the last century. For a detailed compilation, reference may be made to the relevant literature and only common ones are covered herein.

Basically, when journal bearing lubrication is considered, the oil pressure is determined from clearance gap height, relative sliding speed and oil properties. On the other hand, clearance gap height is a result of the equilibrium condition between journal and shell, which is also influenced by contact forces resulting from the oil pressure. Thus, a highly nonlinear coupling exists and iterative procedures are applied to find an overall equilibrium (e.g. [2,3] and references therein).

Analytical solutions to Reynolds equation can be found using long or short bearing theory. For geometries with very short axial extension, overall oil flow is dominated by the Poiseuille flow component in axial direction, whereas its circumferential component is neglected. The resulting analytical solution is called Ocvirk or short bearing solution [4].

In the Sommerfeld solution for infinitely long bearing geometries [5], the pressure gradient in axial direction disappears and axial flow is neglected.

Another class of solution methods is based on characteristic diagrams. Booker [6] first proposed the mobility method of solution for dynamically loaded bearings. It allows for rapid calculation of the journal orbit, pressure field or required oil flow, but is restricted to rigid shells and simple geometry, e.g. circumferential symmetry.

The most accurate and versatile type of solution for an elastohydrodynamic (EHD) evaluation poses a multi-body dynamic (MBD) simulation using flexible components in conjunction with EHD coupling [2]. A detailed simulation of EHD lubrication leading to fluid–structure interaction can be achieved even for intricate geometries containing various oil grooves or variable boundary conditions, which mainly result from the spatial movement of oil supply features, such as a lubrication bore located on the rotating crank shaft, for example.

In addition to these boundary conditions, there are also free boundary conditions in the description of the occurrence of cavitation in the oil film by cavitation models. As opposed to the analytic approaches, which assume a prescribed cavitation boundary, the extent of cavitation region poses a free boundary condition in a detailed EHD simulation. Thus, regardless of the actual solution approach, an iterative determination of the solution is required.

Subsequently, Galerkin method is used in the model order reduction by projection as the periodic nature of combustion engine operation gives rise to new approaches for the generation of projection space and test space. The solution is composed of a modal superposition resulting from the reduced order model. It is spatially evaluated and cavitation boundaries are determined accordingly. These updated boundary conditions are reapplied to the reduced order approach. In addition, the application of the model to the approximation of the Jacobian matrix describing the EHD coupling in a MBD simulation using flexible components is derived and its accuracy is demonstrated. The geometric set-up, which has been used in this paper, is shown in Figure 1 and will be explained in detail in Section 5.

Figure 1. Representation of the main bearing example under investigation.

As covered in detail in the next section, Elrod–Adams’ universal cavitation algorithm is used in this paper [7,8]. An effective density as a second independent variable is introduced besides the pressure, treating cavitation as free boundary condition to Reynolds equation. It then poses a switching-type partial differential equation.

Various commercial MBD tools, such as, for example, AVL EXCITE [3], implement the approaches just mentioned to perform vibration or stress analysis for internal combustion (IC) engines. Steady state investigations of specific operating points (e.g. 3000 rpm under full load in time domain) are common. Such sophisticated simulations are evaluated at the expense of a very time-consuming calculation.

To overcome this drawback, the approach presented in this paper takes advantage of its application to the main bearings of an IC engine. In the framework of a detailed MBD simulation using EHD coupling, the first of multiple combustion cycles (720° crank angle in case of a four-stroke engine) is simulated with full degrees of freedom (DOF). Since a steady state operating point is investigated, the assumption that patterns in the oil film as well as the shape of cavitation region reappear periodically is made. By using the information gathered during the first cycle, the reduced order model is readily available without much additional effort. Therefore, an anticipation of the solution at a certain crank angle in successive cycles can already be attained.

The procedure described exploits periodicity in form of a detailed result of the reference combustion cycle in two ways: first in Section 3.1.1 for the generation of projection basis used in the Galerkin method and second in the choice of an optimal test space for the Petrov–Galerkin method in Section 4.3.

It is obvious that a reduced order model and the reduction of DOF cannot provide a solution with the same level of details as a simulation on the full grid. It should be mentioned that the approach presented in this paper mainly aims at a short-term prediction of pressure and fill ratio over a few degrees of crank angle only and not at a fully autonomous simulation over a whole combustion cycle.

Nevertheless, a cut down on calculation time can be achieved allowing for efficient application of the approximation in the overall MBD simulation. Its convergence properties can be improved by a more frequent recalculation of the contact Jacobian, which is time consuming in general. To introduce the reduced order model, the Jacobian matrix, which per definition quantifies a contact force change per body force variation, is split into two multiplicative parts. One comprises the compliance of the flexible components and the second specifies oil pressure change, respectively, contact force change per gap height variation. Herein only the approximation of the pressure field is sufficient. Details on this approach are given in Section 4.5 and a first result in Section 5.

In Section 2, the cavitation model as well as the model for rough sliding surfaces are quickly reviewed. Section 3 gives an overview of model order reduction, Galerkin method in general and the generation of ansatz space used in this paper. Three different solution procedures are described in Section 4 and results are given in Section 5.

2. Reynolds equation

The pressure field in the incompressible oil film of a slider bearing is governed by a two-dimensional second-order partial differential equation. This relation is called Reynolds equation. For a Newtonian fluid, it is derived from the Navier–Stokes equations in combination with the mass continuity condition (for details cf. [9]).

In Reynolds equation, the thickness of fluid film is assumed small compared to the radius of curvature of the bodies in contact. This allows for the examination on a two-dimensional cartesian coordinate system and conversion from cylindrical coordinates is straightforward. Time is denoted by t and spatial coordinates by $x=\left(x,y\right)$. The direction of relative surface movement corresponds to x, which is equivalent to the circumferential direction in a bearing and is usually quantified in degree shell angle. The axial coordinate is denoted by y. A typical geometry for main bearings will be covered in Section 5.

For any arbitrary given distribution of the clearance gap height $h\left(x,y,t\right)$, dynamic viscosity of the oil $\eta \left(x,y,t\right)$ as well as its density $\rho \left(x,y,t\right)$ and mean relative sliding speed $u_m\left(t\right)$, the pressure field $p\left(x,y,t\right)$ can be determined by Reynolds equation using Dirichlet boundary conditions:

(1) $\frac{\partial }{\partial x}\left(\rho \frac{h^3}{12\eta }\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial y}\left(\rho \frac{h^3}{12\eta }\frac{\partial p}{\partial y}\right)=\frac{\partial }{\partial t}(\rho h)+u_m\frac{\partial }{\partial x}(\rho h).$(1)

Since Reynolds Equation (1) only holds in fully lubricated regions completely filled with oil, the propagation of cavitation can be prescribed in commonly used boundary conditions such as the Sommerfeld, Gümbel or Reynolds condition.

Elrod and Adams [7,8] proposed the Elrod–Adams cavitation model, which is regarded as a reference model in tribology. A saturation function $\Theta \left(x,y,t\right)$ – subsequently called ‘fill ratio’ – is introduced to describe the local ratio of the liquid phase. For incompressible fluids, $\Theta \left(x,y,t\right)$ equals one in lubrication region, representing a gap completely filled with oil. In cavitation regions, the pressure $p\left(x,y,t\right)$ is assumed at a constant cavitation pressure $p_c$. Only shear and squeeze flow exist to fulfil mass continuity and the solution has to be found with respect to $\Theta \left(x,y,t\right)$, which takes values between zero and one. As the extent of cavitation region is not known a priori and changes with the solution, it poses a free boundary condition. Therefore, an iterative evaluation becomes necessary.

The evaluation of slider bearing used herein considers rough sliding surfaces according to the average flow model proposed by Patir and Cheng [10,11]. Pressure flow factors ${\phi }_x\left(h\right)$ and ${\phi }_y\left(h\right)$ as well as a shear flow factor ${\phi }_s\left(h\right)$ describe the influence of sliding surface orientation and mean roughness $\sigma$ on the pressure and fill ratio field. Found by numerical flow simulation on a mockup bearing model, these flow factors are expressed as empirical relations in terms of a normalized clearance gap height $H\left(x,y,t\right)=h\left(x,y,t\right)/\sigma$ and a factor describing surface characteristics [10,11]. They consider that for small clearance gap heights (i.e. high ratios of clearance gap to the mean roughness), the flow is highly influenced by surface roughness, its orientation or even more a micro-textured surface.

Extending (1) by these models leads to the universal averaged Reynolds equation, which holds in full film lubrication as well as in cavitation regions:

(2a) $\frac{\partial }{\partial x}\left(\Theta {\phi }_x\frac{h^3}{12\eta }\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial y}\left(\Theta {\phi }_y\frac{h^3}{12\eta }\frac{\partial p}{\partial y}\right)=\frac{\partial }{\partial t}(\Theta h)+u_m\frac{\partial }{\partial x}(\Theta (h+\sigma {\phi }_s)),$(2a)

(2b) $\begin{array}{cccc}p>p_c& \wedge & \Theta =1\to & solveforp\\ p=p_c& \wedge & 0<\Theta \le 1\to & solvefor\Theta .\end{array}$(2b)

In lubrication regions, $\Theta$ is equal to one and (2a) solved for p, whereas in cavitation region, p is set to the cavitation pressure $p_c$ and $\Theta$ becomes the independent variable.

3. Model order reduction

A well-known method for approximating the solution of partial differential equations is Galerkin’s method. Consider the general differential equation

$D\left(\gamma \left(x,t\right)\right)=f\left(x,t\right)$

defined on domain $\Omega$, denoting the independent variable by $\gamma$, time by t, spatial coordinates by $x$, a general nonlinear operator by D and a function by f. An approximate solution $\tilde{\gamma }$ can be expressed by linear combination of w linearly independent functions ${\beta }_i\left(x\right)$ forming the basis of a vector space $\{{\beta }_i\left(x\right){\}}_{i=1}^w$ also called ansatz space

$\tilde{\gamma }\left(x,t\right)=\sum _{i=1}^w{c}_i\left(t\right){\beta }_i\left(x\right).$

Galerkin method utilizes a test space $\{{\alpha }_i\left(x\right){\}}_{i=1}^w$ to determine temporal coefficients $c_1\left(t\right),\dots ,c_w\left(t\right)$ by the system of w equations

$⟨D\left(\tilde{\gamma }\left(x,t\right)\right),{\alpha }_k\left(x\right)⟩=⟨f\left(x,t\right),{\alpha }_k\left(x\right)⟩k=1,\dots ,w.$

This is equivalent to the requirement that the residuals, which arise from the projection of the original equation onto the ansatz space, weighted by the test space, become zero. Therefore Galerkin’s method is also called method of weighted residuals. For the ansatz space being equal to the test space – the procedure is sometimes referred to as Bubnov–Galerkin method – this amounts to minimizing the squared sum of residuals. For differing spaces, it is called Petrov–Galerkin method.

3.1. Techniques for ansatz space determination

In order to reduce the number of DOF, a projection of the infinite dimensional function space of Reynolds equation onto an ansatz space of lower dimension is performed. General spatial coordinates $x$ now become the cartesian coordinates $x=\left(x,y\right)$. Representation of pressure $p\left(x,t\right)$ and fill ratio $\Theta \left(x,t\right)$ as linear combination of m respectively l arbitrary basis-vectors ${\phi }_j^P\left(x\right)$ and ${\phi }_i^F\left(x\right)$ introduces a residual ${\rho }^P\left(x,t\right)$ and ${\rho }^F\left(x,t\right)$. Modal pressure is denoted by ${\pi }_j\left(t\right)$ and modal fill ratio by ${\theta }_i\left(t\right)$:

(3) $\begin{array}{rl}p\left(x,t\right)& =\sum _{j=1}^m{\phi }_j^P\left(x\right){\pi }_j\left(t\right)+{\rho }^P\left(x,t\right),\\ \Theta \left(x,t\right)& =\sum _{i=1}^l{\phi }_i^F\left(x\right){\theta }_i\left(t\right)+{\rho }^F\left(x,t\right).\end{array}$(3)

Prior to performing model order reduction, a spatial discretization using an equidistant rectangular grid with n nodes is applied to (2). A quantity describing field values such as, for example, the pressure $p\left(x,t\right)$ is now represented by a vector $p\left(t\right)\in R^{n\times 1}$ containing one row for each node. The projection (3) can be written as

(4) $\begin{array}{rl}& p(t)=\sum _{j=1}^m{\phi }_j^P{\pi }_j(t)+{\rho }^P(t)={\Phi }_{1,m}^P\pi (t)+{\rho }^P(t)m\ll n,\\ & \Theta (t)=\sum _{i=1}^l{\phi }_i^F{\theta }_i(t)+{\rho }^F(t)={\Phi }_{1,l}^F\theta (t)+{\rho }^F(t)l\ll n.\end{array}$(4)

Indices in the notation of ansatz spaces ${\Phi }^P$ and ${\Phi }^F$ indicate their dimensions, that is, for the pressure ${\Phi }_{1,m}^P=\left\{{\phi }_1^P\cdots {\phi }_m^P\right\}\in R^{n\times m}$.

3.1.1. Choice of ansatz space

Although the selection of ansatz space is arbitrary in general, an appropriate choice should be made to obtain good results with little DOF. Ansatz functions resembling the most common shapes of the solution efficiently lead to a result, which can be represented by a smaller number of them. From a numerical point of view, the best choice is an orthogonal basis. For example, the use of eigenmodes is well established in linear structural dynamics. As an analytic choice for Reynolds equation, Thomas [12] proposes a combination of Fourier series in circumferential direction and Legendre polynomials in axial direction for the pressure field.

As mentioned in Section 1, the assumption that patterns in the oil film as well as the shape of cavitation region reappear periodically is made. In combination with an existing detailed solution for pressure and fill ratio of the $4\pi$ reference combustion cycle, the ansatz space can be generated in an advantageous way by using principal components. This idea is depicted in Figure 2.

Figure 2. Generation of ansatz space by singular value decomposition (SVD) of the reference cycle.

The full pressure solution at N points in time is collected into a matrix $P\in R^{n\times N}$. For example, if the reference cycle is simulated with the time step size being equal to one degree crank angle rotation, $P$ will contain 720 solution snapshots.

Applying singular value decomposition (SVD) to $P$ yields a diagonal matrix $\Sigma \in R^{n\times N}$ containing the singular values of $P$ and a matrix $U\in R^{n\times n}$. Its columns represent the left singular vectors of $P$:

(5) $P={U\Sigma R}^{\ast }.$(5)

The first m singular vectors in $U$ affiliated with the biggest singular values are chosen as orthonormal basis of the ansatz space ${\Phi }_{1,m}^P$. As the equation for the fill ratio only differs in its notation, the same procedure applies to the generation of ansatz space of the fill ratio.

Evaluation of Reynolds equation transforms from a node-by-node examination of the independent variable to the weighted linear combination of principal modes of solution. The basis vectors are equivalent to those eigenvectors, which represent the optimal basis defined in the Karhunen–Loéve transform, also called Proper Orthogonal Decomposition (POD), frequently used in fluid dynamics. In an implementation where $N<n$, the use of POD would lead to exactly the same ansatz functions, but with increased numerical efficiency [13].

3.2. Modal representation and boundary conditions

Substituting the modal representation of pressure and fill ratio (4) in the discretized form of Reynolds Equation (2) yields one of two equations for every node, either for cavitation or lubrication region.

In lubrication region, the substitution of pressure in (2a) evaluated at node i, leads to

(6) $\begin{array}{rl}& \left[\frac{\partial }{\partial x}\left(\frac{h_i^3}{12{\eta }_i}{\phi }_{x,i}\frac{\partial {\Phi }_{1,m}^P\left(i,:\right)}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{h_i^3}{12{\eta }_i}{\phi }_{y,i}\frac{\partial {\Phi }_{1,m}^P\left(i,:\right)}{\partial y}\right)\right]\pi \\ & =\frac{\partial h_i}{\partial t}+u_m\frac{\partial }{\partial x}\left(h_i+\sigma {\phi }_{s,i}\right)+{\rho }_i^P.\end{array}$(6)

The expression $\partial {\Phi }_{1,m}^P(i,:)/\partial x$ denotes the ith line of a matrix containing the spatial derivative of the ansatz functions with respect to x. Depending on their choice, derivatives are found by difference quotient or as analytic derivatives directly.

Since pressure remains constant in cavitation regions, the left hand side of (2a) becomes zero in this case. Applying the product rule and replacing the temporal derivative $\partial {\Theta }_i/\partial t$ with the backward difference quotient ${\Theta }_i(t_k)-{\Theta }_i(t_{k-1})/\Delta t$ [14] denoting discrete points in time $t_k$ and time step size $\Delta t=t_k-t_{k-1}$, the modal form of Reynolds equation can be written as

(7) $\begin{array}{rl}& \left[u_m\left(h_i+\sigma {\phi }_{s,i}\right)\frac{\partial {\Phi }_{1,l}^F\left(i,:\right)}{\partial x}+\left(\frac{h_i}{\Delta t}+\frac{\partial h_i}{\partial t}+u_m\frac{\partial \left(h_i+\sigma {\phi }_{s,i}\right)}{\partial x}\right){\Phi }_{1,l}^F\left(i,;\right)\right]\theta \\ & =\frac{{\Theta }_i\left(t_{k-1}\right)}{\Delta t}{h}_i+{\rho }_i^F.\end{array}$(7)

Again, the expression $\partial {\Phi }_{1,l}^F(i,:)/\partial x$ denotes the ith line of a matrix containing the spatial derivative of the ansatz functions.

Depending on the location of a node, its corresponding Equation (6) or (7) is added as a row to (8a) or (8b):

(8a) $X\pi =y+{\rho }^P,$(8a)

(8b) $V\theta =w+{\rho }^F.$(8b)

Thus, matrices X and V, vectors y and w contain input quantities, their spatial and temporal derivatives as well as that of the ansatz space. Its rows are described by Equations (6) and (7).

Although (8a) and (8b) seem to be independent of each other on the first sight, they are indirectly coupled as they both depend on the extent of the cavitation region as boundary condition. Besides a Dirichlet boundary condition prescribing pressure and fill ratio on the boundary of the domain, an additional condition must be adhered to. For each node lying in the lubrication region, a row is added to (8a) representing the Reynolds equation. At this node, the fill ratio is required to remain one, resulting in a ‘complementary’ condition. A comparable condition applies to the pressure in cavitation region.

All these boundary and complementary conditions are written in matrix form as well. They are contained in the matrices $X_{BC}$ and $V_{BC}$, whereas the fixed pressure values on the right hand side are denoted by $p_{BC}$. Fill ratio values are set to one at the boundary, which is indicated by the vector $1$. The complete set of $2n$ equations under consideration reads

(9a) $\left[\begin{array}{c}X\\ X_{BC}\end{array}\right]\pi =\left[\begin{array}{c}y\\ p_{BC}\end{array}\right]+{\rho }^P,$(9a)

(9b) $\left[\begin{array}{c}V\\ V_{BC}\end{array}\right]\theta =\left[\begin{array}{c}w\\ 1\end{array}\right]+{\rho }^F.$(9b)

4. Solution procedures

To solve the actually difficult model order reduction procedure with free boundary conditions, an iteration is performed using a spatial evaluation of a superposed modal solution. The projection of the original Equation (2) onto the ansatz space using a reduced number of modes leads to the overdetermined system of Equations (9a) and (9b). Spatial node-by-node representation of both independent variables has been transformed into a modal representation without any spatial information left. Independently applying one of the proposed solution procedures to pressure and fill ratio yields their corresponding modal coordinates. By means of the ansatz functions, the spatial form of pressure and fill ratio solution is recovered afterwards. Based on this information, the cavitation boundary is evaluated and the algorithm advances to the next iteration step.

Galerkin method requires the residual to be orthogonal to a chosen test space. This Galerkin orthogonality condition poses the essential property of the method. To find a solution for the modal coordinates, the orthogonality condition is applied to each of the two residuals as they are defined in Equations (9a) and (9b).

For simplicity, the derivation of the following solutions is shown for pressure equations only, as the equivalent procedure with comparable notation applies to the fill ratio.

4.1. Flat Galerkin method

By choosing the test space being equal to the ansatz space and using only the first $\hat{m}$ of m available columns of ${\Phi }_{1,m}^P$ and $X_{1,m}$ such that all $(m-\hat{m})$ higher modes are neglected, the Flat Galerkin method is obtained:

(10) ${\left({\Phi }_{1,\hat{m}}^P\right)}^T{\rho }^P=0.$(10)

Application to the Reynolds Equation (9) yields the solution for the modal pressure

(11) ${\pi }_{FLG}={\left({\left({\Phi }_{1,\hat{m}}^P\right)}^T\left[\begin{array}{c}{X}_{1,\hat{m}}\\ X_{BC;1,\hat{m}}\end{array}\right]\right)}^{-1}{\left({\Phi }_{1,\hat{m}}^P\right)}^T\left[\begin{array}{c}y\\ p_{BC}\end{array}\right].$(11)

In Equation (11), the boundary conditions are equally weighted as nodes lying within the domain. To emphasize or weaken the influence on the solution of some nodes individually, a diagonal weighting matrix $Q$ could be introduced. By doing so, the residua are weighted according to the entries of $Q$ in advance of the projection onto the test space. A special choice of weighting is made in the following approach.

4.2. Condensed Galerkin method

In the theory of dynamical systems, the Nonlinear Galerkin method is strongly connected to centre and inertial manifolds [15]. Basically, the behaviour of a system is reduced using a projection of higher order dynamics onto a sub-manifold of the phase space, which represents the essential dynamics of the system. The existence of an inertial manifold cannot be proved a priori. Only the assumption is made, that it does exist and an approximate inertial manifold (AIM) can be found by various methods (for an example and references see e.g. [16]).

Concerning the independent variables, the Reynolds equation contains a temporal derivative of the fill ratio on the right hand side of (2a). This term is different from zero, only when a change of fill ratio occurs. As this is mainly the case, when cavitation regions develop or lubrication is reestablished, domains, which are influenced by this term, remain very limited. Therefore, in the time discretized version $\partial \Theta /\partial t$ is represented by a backward difference quotient [14]. The temporal derivative of the gap height can be considered as a given input or can alternatively also be found by backward differences. Basically, the concept of an AIM transforms to static condensation. Although it is based on the same basic concept as Nonlinear Galerkin for nonlinear equations, the approach is denoted Condensed Galerkin in this application.

The solution can be found by splitting the ansatz space ${\Phi }_{1,m}^P$ into two subspaces ${\Phi }_{1,\hat{m}}^P$ and ${\Phi }_{\hat{m}+1,m}^P$. Corresponding modal coordinates are denoted by ${\pi }_{CG}$ and ${\tilde{\pi }}_{CG}$ for higher order modes. Combination of Galerkin orthogonality condition (10) and Reynolds Equation (9) leads to

(12) $\left[\begin{array}{c}{\left({\Phi }_{1,\hat{m}}^P\right)}^T\\ {\left({\Phi }_{\hat{m}+1,m}^P\right)}^T\end{array}\right]\left[\left\{\begin{array}{c}{X}_{1,\hat{m}}\\ X_{BC;1,\hat{m}}\end{array}\right\}\left\{\begin{array}{c}{X}_{\hat{m}+1,m}\\ X_{BC;\hat{m}+1,m}\end{array}\right\}\right]\left[\begin{array}{c}{\pi }_{CG}\\ {\tilde{\pi }}_{CG}\end{array}\right]=\left[\begin{array}{c}{\left({\Phi }_{1,\hat{m}}^P\right)}^T\\ {\left({\Phi }_{\hat{m}+1,m}^P\right)}^T\end{array}\right]\left\{\begin{array}{c}y\\ p_{BC}\end{array}\right\}.$(12)

The relation representing condensed modes can be found by collecting ${\tilde{\pi }}_{CG}$ in the second row of (12)

(13) ${\tilde{\pi }}_{CG}={\left({\left({\Phi }_{\hat{m}+1,m}^P\right)}^T\left\{\begin{array}{c}{X}_{\hat{m}+1,m}\\ X_{BC;\hat{m}+1,m}\end{array}\right\}\right)}^{-1}{\left({\Phi }_{\hat{m}+1,m}^P\right)}^T\left(\left\{\begin{array}{c}y\\ p_{BC}\end{array}\right\}-\left\{\begin{array}{c}{X}_{1,\hat{m}}\\ X_{BC;1,\hat{m}}\end{array}\right\}{\pi }_{CG}\right).$(13)

Substituting higher order modes (13) in the first row of (12) yields

(14) $\begin{array}{rl}Z& =I-\left\{\begin{array}{c}{X}_{\hat{m}+1,m}\\ X_{BC;\hat{m}+1,m}\end{array}\right\}{\left({\left({\Phi }_{\hat{m}+1,m}^P\right)}^T\left\{\begin{array}{c}{X}_{\hat{m}+1,m}\\ X_{BC;\hat{m}+1,m}\end{array}\right\}\right)}^{-1}{\left({\Phi }_{\hat{m}+1,m}^P\right)}^T,\\ {\pi }_{CG}& ={\left({\left({\Phi }_{1,\hat{m}}^P\right)}^{T}Z\left\{\begin{array}{c}{X}_{1,\hat{m}}\\ X_{BC;1,\hat{m}}\end{array}\right\}\right)}^{-1}{\left({\Phi }_{1,\hat{m}}^P\right)}^{T}Z\left\{\begin{array}{c}y\\ p_{BC}\end{array}\right\}.\end{array}$(14)

Compared to the flat procedure result (11), Condensed Galerkin introduces a variable weighting matrix $Z$, which depends on the neglected modes. Since it also contains inputs of the current point in time, it has to be updated online.

4.3. Petrov–Galerkin method

In contrast to the Galerkin procedures described above, Petrov–Galerkin chooses the test space ${\Gamma }^P$ being different from the ansatz space ${\Phi }^P$. Based on the flat approach (11), the solution now is of the general form

(15) ${\pi }_{FLPG}={\left(({\Gamma }_{1,\hat{m}}^P{)}^T\left[\begin{array}{c}{X}_{1,\hat{m}}\\ X_{BC;1,\hat{m}}\end{array}\right]\right)}^{-1}({\Gamma }_{1,\hat{m}}^P{)}^T\left[\begin{array}{c}y\\ p_{BC}\end{array}\right].$(15)

The main focus in this approach lies in the optimal choice of test space. In application to an IC engine, the available information, gathered in a previous combustion cycle by a detailed MBD simulation with EHD coupling, can be exploited by utilizing periodicity. That pressure solution, which arises in an equilibrium condition of the solids involved, is taken as reference. Optimally, the Petrov–Galerkin solution ${\pi }_{PG}$ is equal to that found by projection of such a reference solution $p_{ref}$ and the ansatz functions ${\Phi }_{1,\hat{m}}^P$:

(16) ${\pi }_{PG}={\left(({\Phi }_{1,\hat{m}}^P{)}^T{\Phi }_{1,\hat{m}}^P\right)}^{-1}({\Phi }_{1,\hat{m}}^P{)}^T{p}_{ref}.$(16)

By equating (15) and (16), the definition of an optimal test space is obtained

(17) ${\left({\Gamma }^P\right)}^T\left[y-\underset{z}{\underbrace{X_{1,\hat{m}}{\left(({\Phi }_{1,\hat{m}}^P{)}^T{\Phi }_{1,\hat{m}}^P\right)}^{-1}({\Phi }_{1,\hat{m}}^P{)}^T{p}_{ref}}}\right]=0,$(17)

which can also be expressed by using the substitution $z$ and null space

(18) ${\Gamma }^P\in null(y^T-z^T).$(18)

The dimension of the null space corresponds to the number of nodes $n_R$ at which the Reynolds equation has been evaluated. According to the rank-nullity theorem, null $(y^T-z^T)\in R^{n_R\times (n_R-rank(y-z))}$ holds.

From (15), the additional condition, that the expression $X_{1,\hat{m}}^T{\Gamma }^P$ has to be a $\hat{m}\times \hat{m}$ matrix and invertible, is derived. Using an arbitrary real matrix $\Psi$, the test space ${\Gamma }^P$ can be obtained by

(19) ${\Gamma }^P=null(y^T-z^T)\Psi .$(19)

Two intuitionally chosen possibilities for the determination of $\Psi$ denoting unity matrix $I$ could be

(20) $X_{1,\hat{m}}^{T}null(y^T-z^T)\Psi =I,$(20)

or

(21) $X_{1,\hat{m}}^{T}null(y^T-z^T)\Psi =X_{1,\hat{m}}^T{X}_{1,\hat{m}}.$(21)

The Petrov–Galerkin procedure described above did not consider boundary conditions yet. All matrices under consideration only include nodes in the lubrication region where Reynolds equation applies. The boundary conditions describing pressure in cavitation region and at nodes held at supply or ambient pressure are included as constraints into the optimization.

Solutions showed that best results could be achieved by optimization of the boundary conditions using Reynolds equation as constraint. Since this poses an optimization problem under equality constraints, Reynolds equation is implemented by Lagrange multipliers $\Lambda$ and the objective function is the Lagrange expression

(22) $L=\parallel X_{BC}{\pi }_{PG}-p_{BC}{\parallel }_2^2+2{\Lambda }^T\left({\Gamma }^P{)}^T\right(X_{1,\hat{m}}{\pi }_{PG}-y).$(22)

A possible choice of reference solutions in the Petrov–Galerkin method is shown in Figure 3. A combination of solutions of the previous time steps and the previous combustion cycle can be used.

Figure 3. One possible choice of reference solutions for the generation of Petrov–Galerkin test space ${\Gamma }^P$.

4.4. Cavitation detection

As the extent of cavitation region poses a free boundary condition, it has to be evaluated each time after a new solution has been found. In a full scale evaluation, the algorithm for detecting cavitation and solving Reynolds equation node by node is well established. A description of this algorithm with an application example for micro-textured bearings can be found in [17]. Due to its highly nonlinear nature, the convergence of such a solution approach is analytically hard to prove, if feasible at all. In [18], a detailed study is presented that discusses the convergence behaviour of the iterative procedure. Though no formal proof is given for the reasons stated above, the paper gives a clear indication on the robustness of the algorithm.

Therefore, the cavitation detection procedure for the reduced order model is based on the same principal ideas outlined in [17,18]. To do so, the modal coordinates have to be transformed into a spatial representation of pressure and fill ratio by means of the ansatz functions first. Based on these elements, the occurrence of cavitation has to be reassessed node by node before the next iteration step can be performed with now possibly different boundary conditions.

After a cavitation detection algorithm has decided which nodes belong to the cavitation region, the appropriate Equation (9a) or (9b) is applied in the next iteration step at the corresponding node. At the start of the iteration, this is done based on initial values, which can be taken from the final pressure distribution obtained in the preceding time step, for example.

Due to the approximation character of the solution as well as the optimization procedure, evidence for cavitation is not always distinct. Especially in domains where the transition to cavitation is not sharp, i.e. gradients are small and the pressure only approaches cavitation pressure slowly, detection can become problematic. Here, the pressure and fill ratio results can even lead to contradicting assertions.

To overcome such ambiguities, the fact is exploited, that opposed to the full scale evaluation in the modal evaluation a calculated pressure as well as a calculated fill ratio value exist at the same time for each node. As a minor extension to the known cavitation detection algorithm, basically two fill ratio thresholds are defined. For fill ratio values below a lower bound slightly below one, cavitation is assumed; above an upper bound, lubrication is assumed. If a fill ratio result lies in between these thresholds, additionally the ratio of pressure solution and cavitation pressure will be taken into account to decide on the onset of cavitation. Additional measures could include neighbouring nodes or time history to improve detection rate.

The application of the presented approach in various crank angle positions gained the experience that after typically three iterations the number of nodes in the cavitation region remains unchanged. Usually a convergent behaviour with respect to the number of modes used can be observed. Figure 4 shows a quantitative comparison of the fill ratio result at crank angle ${1240}^{\circ }$ found by Condensed Galerkin with 50 and 75 modes, respectively, compared to the full model to illustrate the usually convergent behaviour with an increasing number of modes.

Figure 4. Cavitation regions obtained by the reduced order model using 50 and 75 modes respectively.

4.5. Approximation of Jacobian matrix using the reduced order model

In the MBD simulation, flexible solids are represented by discretized models. Multiple bodies such as crankshaft and bearing shell are coupled via joint forces, which are determined by the oil pressure in the clearance gap. To achieve an overall equilibrium, external nodal forces calculated for the solids and joint forces have to be equal. Body forces $f$ are considered as predicted and joint forces $r$ as calculated. The latter indirectly also depend on body forces, therefore joint forces are a function of body forces and $r\left(f\right)$ holds [19].

Equilibrium is established when the residual force $ϵ\left(f\right)$ becomes zero:

(23) $ϵ\left(f\right)=r\left(f\right)-f=0.$(23)

In the modified Newton–Raphson approach [19], which is applied in the MBD simulation procedure, the Jacobian matrix $J\left(f\right)$, describing the coupling, arises. Denoting the identity matrix by $I$, it is defined as

(24) $J\left(f\right)=\frac{\partial ϵ\left(f\right)}{\partial f}=\frac{\partial r\left(f\right)}{\partial f}-I.$(24)

Therefore, the calculation of the Jacobian matrix can be reduced to the evaluation of joint force changes per body force variation. In other words, a variation $\Delta f_j$ of the body force at node j (i.e. one entry in $f$) is applied at a time and after solving the joint equation using this varied input, one column of the Jacobian matrix results. Thus, for a complete evaluation, Reynolds equation has to be solved as many times as there exist DOFs of the flexible bodies. As the finite volume grid usually consists of substantially more DOFs than the flexible bodies, a costly solution of Reynolds equation has to be found for each variation anew. Such a repeated evaluation of Reynolds equation to find the Jacobian matrix optimally in every time step poses the bottleneck in the MBD simulation.

The reduced order procedure is intended to decrease calculation time by using an approximation of the pressure solution in the calculation of the Jacobian matrix. The principles of the solution approach are based on analytic differentiation of the reduced order pressure result with respect to the clearance gap height and shall be described in the sequel.

Expanding the first addend on the right hand side of (24) into two parts, leads to

(25) $\frac{\partial r\left(f\right)}{\partial f}=\frac{\partial r\left(f\right)}{\partial h}\frac{\partial h}{\partial f}.$(25)

Herein, the first multiplier describes the joint forces change per gap height variation, whereas the second term can be seen as compliance matrix of the flexible bodies describing the gap height change per body force variation. The latter can be taken from the MBD simulation: however, the former can readily be found using the procedure described in this paper, thus exploiting the reduction of DOFs multiple times per evaluation of the Jacobian matrix.

The result of Reynolds equation is a pressure distribution, which is mapped to the joint forces by a constant transformation matrix $T$ according to $r\left(f\right)=Tp\left(h(f)\right)$. Substituting, for example, the Flat Galerkin reduced order model of the pressure (11) for the first multiplier in (25) yields, with boundary condition entries omitted:

(26) $\frac{\partial r\left(f\right)}{\partial h}=T\frac{\partial p\left(h\left(f\right)\right)}{\partial h}=T\frac{\partial }{\partial h}\left[{\Phi }_{1,\hat{m}}^P{\left(\underset{A}{\underbrace{{\left({\Phi }_{1,\hat{m}}^P\right)}^T{X}_{1,\hat{m}}}}\right)}^{-1}{\left({\Phi }_{1,\hat{m}}^P\right)}^{T}y\right].$(26)

Using abbreviation $A$ in the well-known identity

(27) $\frac{\partial A^{-1}}{\partial h_i}=-A^{-1}\frac{\partial A}{\partial h_i}{A}^{-1}$(27)

one row in (26) can be rewritten as

(28) $\frac{\partial r\left(f\right)}{\partial h_i}=T{\Phi }_{1,\hat{m}}^P{A}^{-1}{\left({\Phi }_{1,\hat{m}}^P\right)}^T\left[\frac{\partial y}{\partial h_i}-\frac{\partial X_{1,\hat{m}}}{\partial h_i}{A}^{-1}{\left({\Phi }_{1,\hat{m}}^P\right)}^{T}y\right].$(28)

The same procedure can be applied to Condensed Galerkin or Petrov–Galerkin solutions and leads to more complicated terms, which are not covered in this paper, though.

Using the scheme (28) in conjunction with the Galerkin reduced order model and the compliance matrix, the Jacobian matrix can be determined at a drastically reduced computational cost. The evaluation of a first academic example in Section 5.1 yields promising results, which encourage further research into this application.

5. Numerical results

The results presented in this section have been found for a single crankshaft main bearing for an inline IC engine of typical dimensions as it is depicted in Figure 6. Its diameter is $60mm$, its width $19mm$. The bearing shell is divided into two segments. On the upper half of the shell – in relation to the orientation of the cylinder – an oil groove is located at the middle plain extending over ${180}^{\circ }$ in circumferential direction. The split line between upper and lower shell also serves as oil supply over the whole bearing width. Both are indicated by dashed lines in the following figures. The oil (SAE 5-W30) is assumed to remain at a constant temperature of ${105}^{\circ }C$, which corresponds to a dynamic viscosity of $8.0625\times {10}^{-3}Pas$. Surface roughness is assumed isotropic (i.e. ${\phi }_x={\phi }_y$) with a composite roughness of 0.67 μm.

Figure 5. Decay of normalized singular values of pressure and fill ratio ansatz functions.

Figure 6. Bearing geometry with pressure contour in lower shell area and circumferential cut position.

As starting point for reduced order solutions and the generation of ansatz space, a detailed MBD simulation with EHD coupling at $5000rpm$ using the method, which is proposed in [2], has been performed over ${1440}^{\circ }$ crank angle. A 3D volume model of the bearing shell with 336 nodes is coupled with a structured beam mass model of the journal with 7 nodes (Figure 1). In particular, 336 radial flexible shell DOFs are coupled with 2 DOFs per connected journal node, which results in 14 flexible journal DOFs in total. The EHD calculation is done using a finite volume grid with 1885 cells. The coupling of the MBD and the EHD domain is done using a modified Newton–Raphson approach [19].

Of these two combustion cycles, the first one is used to overcome transient effects due to missing initial conditions before returning valid results in the full procedure. The solutions in this initialization phase, containing 720 different crank angle positions of the first combustion cycle is used as collection of snapshots for the generation of ansatz functions for pressure and fill ratio. Here, the periodicity of the overall simulation is exploited as after the initialization cycle a reduced order model approximation is already available for subsequent cycles.

Figure 5 shows the decay of the normalized singular values of pressure and fill ratio ansatz functions. Subsequent results have been found by using bases containing 150 modes. In this case, this equals a decay of the normalized singular values below $5\times {10}^{-4}$.

The reduced order model utilizes the detailed result at crank angle ${970}^{\circ }$ as initial value and results are presented at ${971}^{\circ }$. To consider a change in the extent of cavitation region, two iterations are performed. When more iterations are applied, the number of cavitation nodes remains constant, thus indicating that the algorithm converged. The generation of Petrov–Galerkin test space is based on the full solution at crank angles $\left[{250}^{\circ },{251}^{\circ },{252}^{\circ }\right]$.

In Figure 6, the bearing geometry with oil supply groove and split line is shown. The pressure solution at crank angle ${971}^{\circ }$ is indicated by contour lines, cut position of Figures 9 and 10 as crosses.

Figure 7. Comparison of pressure solutions of the full model (b) and Condensed Galerkin (a).

Figure 8. Comparison of fill ratio solutions of the full model (b) and Condensed Galerkin (a).

Figure 9. Pressure solutions found by Flat Galerkin, Condensed Galerkin and Petrov–Galerkin compared to the full solution shown as circumferential cut.

Figure 10. Fill ratio solutions found by Flat Galerkin, Condensed Galerkin and Petrov–Galerkin compared to the full solution shown as circumferential cut.

Bearing geometry is unrolled into a two dimensional representation in Figure 7, comparing the pressure solution of the full model to the reduced order model solution, which has been found by Condensed Galerkin. Peak pressure amounts to 9.2 MPa. Both solutions are of very similar shape and differ only slightly.

The same type of comparison for the fill ratio is shown in Figure 8. The intensity of grey corresponds to the strength of cavitation. Minimum fill ratio turns out to be 0.0156. Although the results are very similar in general, larger deviations can be observed than in the pressure solution.

A more detailed graphical evaluation of results offers a circumferential cut, whose position relative to the bearing can be seen in Figure 6. Each marker in Figure 9 represents the pressure solution at one of the nodes with same axial coordinate.

In relation to the full solution, Flat Galerkin with 150 modes in the ansatz space leads to a unusable result. The two other procedures reproduce the full solution well, thereby Condensed Galerkin slightly outperforms Petrov–Galerkin. A general trend in the review of various solutions is the fact, that an influence of the reference solution used in the generation of test space for Petrov–Galerkin on its reduced order solution is identifiable.

The same conclusion can be made for the fill ratio solution depicted in Figure 10. In contrast to the pressure solution, Flat Galerkin leads to a good result in this evaluation although this is not the case in general. The tendency of Petrov–Galerkin reproducing reference solutions is again evident in the evaluation of fill ratio.

For a quantitative comparison of solutions, the sum of absolute errors ${\sum }_i|e_i|={\sum }_i|p_i-p_{full,i}|$; the average absolute error per node ${\sum }_i|e_i|/n$ as well as the overall sum of pressure ${\sum }_i{p}_i$ and fill ratio values ${\sum }_i{\Theta }_i$ are listed in Table 1. It can be observed that Condensed Galerkin leads to the lowest absolute and relative errors.

Table 1. Quantitative comparison of solution methods for pressure and fill ratio results.

	Pressure (MPa)			Fill ratio (–)
	${\sum }_i|e_i|$	${\sum }_i|e_i|/n$	${\sum }_i{p}_i$	${\sum }_i|e_i|$	${\sum }_i|e_i|/n$	${\sum }_i{\Theta }_i$
Full	−	−	2165.9	−	−	1685.8
Condensed	95.141	0.051	2115	10.646	5.648e-3	1681.2
Petrov–Galerkin	194.76	0.103	2270	13.439	7.129e-3	1690.1
Flat	1053.1	0.559	1459.5	17.138	9.092e-3	1679.2

5.1. Approximation of Jacobian matrix

While in the preceding section the approximation capabilities of the Galerkin reduced order model have been discussed and illustrated, the next step will now be to demonstrate its applicability in the fast determination of the Jacobian matrix. Note that the result, which is shown below, still has to be multiplied with the compliance matrix to yield the final Jacobian matrix.

In this context, lines of action are defined for nine contact forces $F_1$ to $F_9$, which are unevenly distributed around the circumference. They are depicted in Figure 11 in addition to a simplified representation of geometry, which is considered rigid in this case. Moreover, supplemental contact forces $F_{1\perp }$ to $F_{9\perp }$ are defined, whose lines of action are passing through the origin of $\left({\zeta }_1,{\eta }_1\right)$ and which are at the same time perpendicular to those corresponding to $F_1$ to $F_9$. That is, $F_{1\perp }$ is pointing into the ${\eta }_1$ direction, for example.

Figure 11. Axial view of the bearing with lines of action of the contact forces shown as dotted lines.

The intensity of each of these 18 forces is determined by integrating the pressure, which is projected onto the corresponding direction, over the circumference. As such a definition of forces is usually not used in a MBD simulation, the conducted comparison has to be considered as academic result to basically show the ability of the method.

Further, nine variations of gap height are applied. For that, the journal is translated as a whole by $\Delta h$ into directions ${\alpha }_1$ to ${\alpha }_9$, wherein $\alpha$ describes the nine lines of action of $F_1$ to $F_9$. Corresponding variations are denoted as $h_{\alpha 1}$ to $h_{\alpha 9}$. In this example, $\Delta h$ is chosen as $5\%$ of the smallest gap height, which amounts to $5\times {10}^{-8}$ m. The first part in the Jacobian matrix describing the contact force change per gap height variation is determined in two ways.

First, a pressure solution using the full basis is found for each gap height variation. The expression $\partial r\left(f\right)/\partial h\in R^{18\times 9}$ is calculated by numerical differentiation, i.e. subtracting the original reference solution without variation from the currently found solution with variation and then dividing by the gap height variation node by node. Thus, a numerically derived matrix $J_{\Delta }$ is obtained based on full solutions.

Second, the same expression is determined by using the reduced order model according to (28). This solution, denoted by $J_R$, therefore can be found with reduced effort and is based on $\hat{m}$ modes.

To compare $J_{\Delta }$ and $J_R$, an SVD is conducted. Both matrices should optimally describe the change of contact forces in their intensity and principal direction equally well. Singular values ${\sigma }_{\Delta ,j}$ should match with ${\sigma }_{R,j}$ and right-singular vectors should point into the same direction, that is their inner product $⟨{\text{silon}}_{\Delta ,j},{\text{silon}}_{R,j}⟩$ should be one or minus one.

Table 2 shows such a comparison for the configuration described. The first two rows show the normalized singular values using ${\sigma }_{\Delta ,1}$. Therein, the first two values nearly match exactly. The others come out to be several magnitudes smaller, although they all still lie closely together. In the last row, the absolute values of the inner product of right-singular vectors $⟨{\text{silon}}_{\Delta ,j},{\text{silon}}_{R,j}⟩$ are shown. For the dominant singular values, their direction is again basically equal (i.e. their inner product is very close to one).

Table 2. Comparison of singular values and inner products of right-singular vectors of the two matrices $J_{\Delta }$ and $J_R$ describing contact force changes per gap height variation.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
${\sigma }_{\Delta ,j}/{\sigma }_{\Delta ,1}$	1.000	0.448	8.314e-17	2.844e-17	1.574e-17
${\sigma }_{R,j}/{\sigma }_{\Delta ,1}$	0.964	0.436	8.06e-17	3.706e-17	1.782e-17
$|⟨{\text{silon}}_{\Delta },{\text{silon}}_R⟩|$	0.9985	0.9857	0.9362	0.9594	0.9349

Considering these first results, the future application of a reduced order evaluation in the Jacobian matrix approximation seems promising.

6. Outlook

The methods proposed in this paper can reduce the number of DOF in the solution of Reynolds equation by approximately up to $90\%$ depending on the operating conditions and point in time under consideration. At the same time, they preserve good approximation quality characteristics in short-term calculations. Basically, the algorithm is applicable at every load level without adaption or non-physical input, except the dimension of ansatz space to be used.

In the multi-body simulation of IC engines, different contact models can be selected depending on the requirements. Compared to a powerful solver on today’s hardware, replacement of the full EHD evaluation of radial slider bearings under consideration of local effects and cavitation is not feasible with a reduced order model based on projection. However, switching from full to reduced order procedure and backwards can readily be implemented and a fast approximation of solution is available, which can be used as intermediate step in a Runge–Kutta method, for example. Application of the proposed method is directly possible to axial slider bearings or lubricated piston liner contacts as well.

As it has been shown, the approximation can beneficially be exploited in the calculation of the Jacobian matrix for coupling of MBD and EHD lubrication, which poses a key factor in the performance of the overall multi-body simulation procedure. To avoid time-consuming recalculations in the current algorithm, the Jacobian is kept constant over several time steps until a criterion like an alteration of damping coefficients or the exceeding of a maximum number of iterations is violated. Acceleration of the recalculation procedure of the Jacobian, which describes the coupling, using a reduced order model solution of the Reynolds equation could allow for more frequent updates and subsequently for better overall convergence properties. Since the pressure field is readily available, force derivatives can directly be applied to the pressure, leading to a reduced effort in their evaluation.

Another possible field of application is noise, vibration and harshness simulation, where a detailed resolution of cavitation region or fill ratio is not necessary and could be neglected therefore. In addition, these results could be used in an identification of coefficients needed in a spring-damper representation of the bearing.

Although further investigation especially of the properties of Petrov–Galerkin and a possibility for error diagnosis has to be performed, the methods proposed in this paper proved well in their designated implementation.