An extended Elrod-Adams model to account for backpressure and blow-by inception

TL;DR

This paper extends the Elrod-Adams model for piston-ring-liner systems to include backpressure effects and blow-by prediction, improving the accuracy of hydrodynamic and friction force assessments in engines.

Contribution

It introduces a novel model that incorporates combustion chamber backpressure and a criterion for blow-by inception into the Elrod-Adams framework.

Findings

Backpressure significantly affects predicted friction and MFT.

The extended model accurately predicts blow-by inception.

Inclusion of pressure difference improves hydrodynamic assessments.

Abstract

The Piston-Ring-Liner system is the main tribological component of internal combustion engines. The Elrod-Adams model is customarily used to numerically assess the hydrodynamics of different ring designs and liner surface treatments. However, that model does not incorporate the backpressure boundary condition, which in this case corresponds to the combustion chamber pressure and is quite significant. In this contribution a model is proposed that imposes the combustion-chamber pressure in a mass-conserving way, together with an effective algorithm for its numerical approximation. The new model incorporates the pressure difference across the ring, which is shown to have a substantial effect on the predicted friction force and MFT. The model is further elaborated so as to provide a criterion for predicting blow-by inception.

Figures

Tables

Table 1. Table 1: Fundamental scales.

Symbol	Value	Units	Description
$H$	$1\times {10}^{-6}$	m	Gap thickness
$L$	$1\times {10}^{-3}$	m	Scale on the $x_1$-$x_2$ plane
$\mu$	$4\times {10}^{-3}$	Pa$\cdot$s	Lubricant viscosity
$U$	$10$	m/s	Reference speed along the $x_1$ axis

Table 2. Table 2: Nomenclature.

Quantity	Description	Scale
$x_1$	Longitudinal coordinate	$L$
$x_2$	Transverse coordinate	$L$
$t$	Time	$L/U$
$h,h_{\text{ring}},h_{\text{wear}}^{\delta }$	Functions defining the surfaces gap	$H$
$h_{\text{feed}}$	Film thickness at the boundaries	$H$
$u$	Liner’s speed along $x_1$	$U$
$Z$	Minimum gap in time	$H$
$L_{\text{ring}}$	Ring’s length along $x_1$	$L$
$B_{\text{r}}$	Bore radius	$L$
$p$	Hydrodynamic pressure	$6\mu UL/H^2$
$W^{\text{h}},W^{\text{a}},W^{\text{c}}$	Radial forces per unit width	$6\mu U{L}^2/H^2$
$m$	Ring’s mass per unit width	$6\mu L^4/\left(H^{3}U\right)$
${\mu }_{\text{c}}$	Boundary friction coefficient ($=0.11$)	-

Table 3. Table 3: Minimum value of δ 𝛿 \delta for which the simulations fail for every δ ≥ δ max 𝛿 subscript 𝛿 max \delta\geq\delta_{\textnormal{\scriptsize{max}}} .

	$\Delta h$($R$)
	0.5(250)	1.0(125)	1.5(83.3)	2.0(62.5)	2.5(50)	3.0(41.7)	3.5(35.7)	3.91(32.0)	4.0-8.0(31.25-15.6)
${\delta }^{\text{max}}$	0.025	0.035	0.040	0.045	0.050	0.055	0.025	0.005	0

Equations

$$h(x_1,x_2,t)=Z(t)+h_{\text{ring}}^R(x_1)+h_{\text{wear}}^{\delta }(x_2),$$

$$h_{\text{ring}}^R(x_1,x_2)=\frac{L}{H}\frac{{\left(x-0.5\right)}^2}{2R},$$

$$h_{\text{wear}}^{\delta }(x_2)=\delta \exp \left(-(x_2-\pi B_{\text{r}}{)}^2/c^2\right)x_2\in [0,2\pi B_{\text{r}}].$$

$$m\frac{d^{2}Z}{d{t}^2}=W^{\text{h}}(t)+W^{\text{ps}}+W^{\text{con}}(h)+W^{\text{cc}}(t).$$

$$W^{\text{h}}(t)=\frac{1}{B_r}{\int }_0^{2\pi B_{\text{r}}}{\int }_a^{b}p(x_1,x_2,t)d{x}_{1}d{x}_2,$$

$$W^{\text{con}}(h)=\frac{1}{B_r}{\int }_0^{2\pi B_{\text{r}}}{\int }_a^b{p}^{\text{con}}(h(x_1,x_2,t))d{x}_{1}d{x}_2,$$

$$p^{\text{con}}(h)=\pi \frac{16\sqrt{2}}{15}{\left(\eta \beta \sigma \right)}^2\sqrt{\sigma /\beta }{E}^{\prime }{F}_{5/2}\left(h/\sigma \right),$$

$$W^{\text{cc}}(t)=-\gamma p_{\text{cc}}(t)L_{\text{ring}}$$

$$F(t)=\frac{1}{B_r}{\int }_0^{2\pi B_{\text{r}}}{\int }_a^b\left(\frac{\mu u}{h}g(\theta )-\frac{h}{2}\frac{\partial p}{\partial x_1}-p\frac{\partial h_{\text{ring}}}{\partial x_1}+{\mu }_{\text{c}}{p}^{\text{con}}\right)d{x}_{1}d{x}_2,$$

$$g(\theta )=\{\begin{array}{ll}\theta & \text{mbox}if\theta >{\theta }_{\text{s}}\\ 0& \text{mbox}otherwise\end{array},$$

$$\nabla \cdot \left(\frac{h^3}{2}\nabla p\right)=\frac{u(t)}{2}\frac{\partial h\theta }{\partial x_1}+\frac{\partial h\theta }{\partial t},$$

$$p\ge T(\theta )\text{mbox}and0\le \theta \le 1$$

$$\{\begin{array}{ll}p(0,x_2)=0,& p(1,x_2)=p_{\text{cc}}(t)\\ \theta (0,x_2)=\frac{h_{\text{feed}}}{h(0,x_2)},& \theta (1,x_2)=\frac{h_{\text{feed}}}{h(1,x_2)}\end{array}\text{mbox}for{x}_2\in [0,2\pi B_{\text{r}}]$$

$$\left(p-T(\theta )\right)\left(1-\theta \right)$$

$$T(\theta )(\mathbf{x})=\{\begin{array}{ll}{p}_{\text{cc}}(t)& \text{mbox}for\mathbf{x}\in {\Omega }_0^{\text{r}}\\ 0& \text{mbox}otherwise\end{array}.$$

$$\frac{\partial p}{\partial \mathbf{n}}=2\frac{(1-{\theta }_{-})}{h^2}\left(\frac{U}{2}{\mathbf{e}}_{\text{c}}-{\mathbf{V}}^{\text{i}}\right)\cdot \mathbf{n},$$

$$\{\begin{array}{ll}\frac{\partial p}{\partial \mathbf{n}}\ge 0& \text{mbox}atrupturepoints\\ \frac{\partial p}{\partial \mathbf{n}}\le 0& \text{mbox}atreformationpoints\end{array},$$

$$\frac{\partial p}{\partial \mathbf{n}}\le 0\text{mbox}in\partial {\Omega }_{+}.$$

$$-\frac{h^3}{2}\frac{\partial p}{\partial x_1}+\frac{u}{2}h\theta \simeq -\frac{1}{2}\frac{(h_{i-1,j}^n{)}^3+(h_{i,j}^n{)}^3}{2}\frac{p_{i,j}^n-p_{i-1,j}^n}{\Delta x_1}+\frac{u}{2}{h}_{i-1,j}^n{\theta }_{i-1,j}^n,$$

$$a_{i,j}^{00}{p}_{i,j}^n+e_{i,j}^{00}{\theta }_{i,j}^n$$

$$p_{i,j}^n$$

$$0\le {\theta }_{i,j}^n$$

$$\left(p_{i,j}^n-T{\left({\mathbf{\theta }}^n\right)}_{i,j}\right)\left(1-{\theta }_{i,j}^n\right)$$

$$C_{i,j}({\mathbf{p}}^n,{\mathbf{\theta }}^n)=$$

$$\begin{array}[]{ll}a_{i,j}^{00}=s_{i+1,j}^{n}+s_{i-1,j}^{n}+(\Delta x_{1}/\Delta x_{2})^{2}\left(s_{i,j+1}^{n}+s_{i,j-1}^{n}\right),&e_{i,j}^{00}=\left(u\,\Delta x_{1}+2\Delta x_{1}^{2}/\Delta t\right)h^{n}_{i,j},\\ a_{i,j}^{+0}=-s^{n}_{i+1,j},&a_{i,j}^{-0}=-s^{n}_{i-1,j},\\ a_{i,j}^{0-}=-(\Delta x_{1}/\Delta x_{2})^{2}\,s^{n}_{i,j-1},&a_{i,j}^{0+}=-(\Delta x_{1}/\Delta x_{2})^{2}\,s^{n}_{i,j+1},\\ e_{i,j}^{-0}=-u\,\Delta x_{1}\,h^{n}_{i-1,j},&s^{n}_{i\pm 1,j\pm 1}=\frac{1}{2}\left((h_{i,j}^{n})^{3}+(h_{i\pm 1,j\pm 1}^{n})^{3}\right),\\ f_{i,j}^{n}=2\Delta x_{1}^{2}/\Delta t\,h_{i,j}^{n-1}\theta_{i,j}^{n-1}.&\end{array}$$

$$B_T(\mathbf{p},\mathbf{\theta }{)}_{i,j}=\{\begin{array}{ll}\left(\frac{C_{i,j}(\mathbf{p},\mathbf{\theta })-e_{i,j}^{00}}{a_{i,j}^{00}(T\mathbf{\theta }{)}_{i,j}},1\right)& \text{if}\frac{C_{i,j}(\mathbf{p},\mathbf{\theta })-e_{i,j}^{00}}{a_{i,j}^{00}(T\mathbf{\theta }{)}_{i,j}}\ge 0\\ \left(T(\mathbf{\theta }{)}_{i,j},\frac{C_{i,j}(\mathbf{p},\mathbf{\theta })-a_{i,j}^{00}T(\mathbf{\theta }{)}_{i,j}}{e_{i,j}^{00}}\right)& \text{if}\frac{C_{i,j}(\mathbf{p},\mathbf{\theta })-e_{i,j}^{00}}{a_{i,j}^{00}(T\mathbf{\theta }{)}_{i,j}}<0\end{array},$$

$$B_T(\mathbf{p},\mathbf{\theta })=(\mathbf{p},\mathbf{\theta }).$$

$$C_{i,j}(\mathbf{p},\mathbf{\theta })-a_{i,j}^{00}T(\mathbf{\theta }{)}_{i,j}=f_{i,j}+S\Delta x_1{h}_{i-1,j}{\theta }_{i-1,j}+{\Delta }_{i,j}^{\prime \prime },$$

$${\Delta }_{i,j}^{\prime \prime }$$

$$+q^2\left\{s_{i,j-1}\left(p_{i,j-1}-T(\mathbf{\theta }{)}_{i,j}\right)+s_{i,j+1}\left(p_{i,j+1}-T(\mathbf{\theta }{)}_{i,j}\right)\right\},$$

$$T_{ϵ}(\mathbf{\theta }{)}_{i,j}=\{\begin{array}{ll}{p}_{\text{cc}}(t)& \text{ if }(i,j)\in I_r\cup V^{ϵ}\\ 0& \text{ if }(i,j)\notin I_r\cup I^{ϵ}\end{array},$$

$$\frac{{\partial }^{2}p}{\partial x_1^2}\le 0\text{mbox}andp>0$$

Full text

Keywords: Piston-Ring-Liner, Combustion Chamber Pressure, Elrod-Adams model, Blow-by.

1 Introduction

Current hydrodynamical models allow to estimate load carrying capacity, friction losses and minimum film thickness (MFT) of tribological devices that involve challenging physical phenomena such as cavitation, which is a non-linear rupture of the fluid continuity [1]. Regarding the Piston-Ring-Liner (PRL) mechanism, many studies have addressed its tribological performance setting the Combustion Chamber Pressure (CCP) as a boundary condition along the Reynolds cavitation model for pressure [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. For instance, Kligerman et al [10] compared the predicted friction force and MFT between two models: The “Real engine simulation” where the CCP was considered as a boundary condition, and the “Approximate solution”, based on a quasi-static simulation where a boundary condition equal to the ambient pressure was imposed. The authors found that the MFT was strongly affected by the CCP time variation, and that the “Approximate solution” underestimated up to 30% the time-averaged friction force.

However, mass-conservation is a mandatory property of numerical schemes when considering, as is the current trend, the microscopic features of liner surfaces [14, 15]. Recently, many studies on the PRL performance have adopted the mass-conserving Elrod-Adams model [16, 17, 18, 19, 20, 21, 22, 23]. Unfortunately, that model assumes the boundary condition for pressure as being equal to the cavitation pressure, which is invalid for the PRL mechanism since the CCP can assume values in the order of 100 atm.

In this work, an extension of the Elrod-Adams model for non-homogeneous boundary conditions is proposed based on the numerical model presented in [24, 25]. Friction and MFT for a textured ring are computed with the extended model and a numerical comparison with both the Elrod-Adams and the Reynolds models is performed.

A particular characteristic of the PRL operating conditions is that gas leakage can take place through the gaps between piston and ring and/or between ring and liner. This undesirable effect is known as blow-by and is associated to power-loss and air pollution [26]. Performing two-dimensional simulations it is shown here that the extended model can be used to predict blow-by inception, which manifests itself as the development of channels of gas/vapor through the gap between the ring and the liner.

This paper’s structure is as follows: In Section 2 the extended Elrod-Adams model is described in terms of a restriction on the pressure gradient at the rupture boundary, and typical models for the PRL geometry and the ring dynamics are introduced. The numerical method to solve the extended model coupled to the ring dynamics is presented in Section 3. A comparison of the resulting tribological performance between the extended model and both the Elrod-Adams and the Reynolds model is performed in Section 4. In that section we also expose predictions of an encounter between the reformation and the rupture boundaries for rings of different degree of conformity, which is put forward as a criterion for blow-by. Concluding remarks are drawn in the final section.

2 Mathematical models

The main components of the PRL system are depicted in Fig. 1, here a simplified modeling of that system is assumed as described in the next section. The fundamental scales are given in Table 1 and the nomenclature for the involved quantities is given in Table 2.

2.1 Geometrical model and ring dynamics

Both surfaces (ring and liner) are considered to be smooth (no surface roughness) for simplicity. The gap between the surfaces is modeled by the function

[TABLE]

where $Z$ is the dynamic variable accounting for the expansion/contraction of the ring, $h_{\textnormal{\scriptsize{ring}}}$ is the function describing the ring’s shape, defined as satisfying $\min_{x\in[0,1]}h^{R}_{\textnormal{\scriptsize{ring}}}(x_{1})=0$, depicted in Fig. 2 and written as

[TABLE]

and $h^{\delta}_{\textnormal{\scriptsize{wear}}}$ is a function modeling a non-homogeneous wear of the surfaces along the circumferential direction, also satisfying $\min_{x\in[0,2\pi B_{\textnormal{\scriptsize{r}}}]}h^{\delta}_{\textnormal{\scriptsize{wear}}}(x_{2})=0$ and assumed for simplicity to have a Gaussian profile centered at the circumferential angle equal to $\pi$:

[TABLE]

The ring follows a reciprocal motion along the $x_{1}$-axis and its speed is denoted by $u(t)$. The dynamics of $Z(t)$ is determined by Newton’s equation:

[TABLE]

The first term on the left-hand side of this equation corresponds to the hydrodynamical force per unit width (in the circumferential direction), computed as

[TABLE]

the second and third terms correspond to the applied load (datum) and the contact load, respectively, the latter being computed by means of the Greenwood-Tripp model [27] as

[TABLE]

with the contact pressure $p^{\textnormal{\scriptsize{con}}}$ given by

[TABLE]

where $E^{\prime}$ is the composite modulus of elasticity, and $\eta$, $\beta$ and $\sigma$ are the asperity density, radius of curvature and composite standard deviation of the asperities height. The value of $F_{5/2}\left(h/\sigma\right)$ corresponds to an integral that is approximated as done by Panayi-Schock in [28]. The last term in Eq. (4) corresponds to the pressure exerted on the back of the ring, pushing it towards the liner and modeled by

[TABLE]

where $\gamma$ is the Back-pressure factor and $L_{\textnormal{\scriptsize{ring}}}$ is the ring’s length along the $x_{1}$ direction.

The friction force exerted by the fluid on the ring at time $t$ along the $x_{1}$-axis is given by [29]

[TABLE]

where

[TABLE]

and $\mu_{\textnormal{\scriptsize{c}}}$ is the boundary friction coefficient. The threshold $\theta_{\textnormal{\scriptsize{s}}}$ is interpreted as the minimum oil fraction for shear forces to be transmitted between the surfaces and it was set to 0.95 (e.g., [18]).

2.2 Hydrodynamics and cavitation modeling

The hydrodynamical pressure $p$ and the saturation field $\theta$ are modeled by an extension of the Elrod-Adams cavitation model that reads

[TABLE]

along the restrictions for the unknown fields:

[TABLE]

the boundary conditions:

[TABLE]

and the complementarity conditions:

[TABLE]

The novelty of this model lies in the operator $T$. When $T(\theta)\equiv 0$ one recovers the classical Elrod-Adams model, which essentially assumes a zero combustion-chamber pressure. Let us denote $\Omega_{0}=\{\mathbf{x}\in\Omega:\theta(\mathbf{x})<1\}$, the region of the domain where the film is not complete, which in general has many connected components. The basic idea of the model proposed here is that the connected component of $\Omega_{0}$ that touches the right boundary must have $T$ equal to the combustion-chamber pressure, while in the others $T\equiv 0$. More specifically, denoting by $\Omega_{0}^{\textnormal{\scriptsize{r}}}$ the connected component of $\Omega_{0}$ that intersects the set $\{(x_{1},x_{2})\in\Omega:x_{1}=1\}$, we define

[TABLE]

An example of this is shown in Fig. 3. Notice that if the right and left sides of the system are connected by a “channel” of imcomplete film, i.e., if the intersection $\Omega_{0}\cap\Omega_{0}^{\textnormal{\scriptsize{r}}}$ is non-empty and $p_{\textnormal{\scriptsize{cc}}}(t)>0$, the model becomes undefined. This condition intuitively determines the blow-by of combustion chamber gases, and will be so considered hereafter.

As in the case of the standard Elrod-Adams model, the extended model consists of a transport equation on $\Omega_{0}$ and an elliptic equation on the full-film region. Though the standard model (interpreted in the distribution sense) enforces mass conservation throughout the domain, when $p_{\textnormal{\scriptsize{cc}}}>0$ a supplementary boundary condition must be supplied at rupture points. In fact, on the interface $\partial\Omega_{+}$ the mass conservation condition can be written, assuming $h$ to be continuous, as

[TABLE]

where $\theta_{-}$ is the value of $\theta$ at the cavitated/separated side, $\mathbf{e}_{\textnormal{\scriptsize{c}}}$ is the unit vector along the Couette flux, $\mathbf{V}^{\textnormal{\scriptsize{i}}}$ is the interface velocity and $\mathbf{n}$ is the outward normal to $\partial\Omega_{+}$. This last condition implies

[TABLE]

which when $T\equiv 0$ corresponds to the null-gradient condition at the rupture points but, when $p_{\textnormal{\scriptsize{cc}}}>0$, is not enough to assure the well-posedness of the extended model [30, 31, 32]. To tackle this issue, in this work the condition $\frac{\partial p}{\partial\mathbf{n}}=0$ is imposed at the rupture points. This can be summarized as

[TABLE]

3 Numerical Method

The domain $\Omega$ is divided in $N_{x_{1}}\times N_{x_{2}}$ cells and a finite volume scheme for Eq. (11) is used. For this, the $x_{1}$ flux component going from node $(i-1,j)$ to node $(i,j)$ is discretized by

[TABLE]

where an upwind approximation is used for the Couette term. Balancing the fluxes on each cell and discretizing time along an implicit scheme for the temporal term $\frac{\partial h\theta}{\partial t}$ one gets the system of equations:

[TABLE]

where

[TABLE]

with

[TABLE]

Notice that $a^{00}_{i,j},e^{00}_{i,j},f_{i,j}^{n}$ are non-negative, and each term $a^{+0}_{i,j}$, $a^{-0}_{i,j}$, $a^{0+}_{i,j}$, $a^{0-}_{i,j}$, $e_{i,j}^{-0}$ is non-positive, thus $C_{i,j}(\bm{p}^{n},\bm{\theta}^{n})\geq 0$. This type of discrete system was already studied for the case $T\equiv 0$ (standard Elrod-Adams algorithm) in [33] by means of a fixed-point method and in [34] where the convergence of some Gauss-Seidel-like algorithms was proved.

Dropping the superscript $n$ and denoting by $N$ the number of unknowns, the system of equations (19)-(22) can be written as a fixed-point problem by defining the operator $B_{T}:\mathbb{R}^{N}\times\mathbb{R}^{N}\mapsto\mathbb{R}^{N}\times\mathbb{R}^{N}$ as

[TABLE]

This way, Equation (19) can be rewritten as the fixed point problem

[TABLE]

Let us denote by $I=I_{0}\overset{\circ}{\cup}I_{\textnormal{r}}$ the set of cells indexes for which $\bm{p}$ and $\bm{\theta}$ are unknowns, where $I_{\textnormal{r}}$ is the subset of cells in the region $\Omega_{0}^{\textnormal{r}}$ and $I_{0}=I\setminus I_{\textnormal{r}}$. A solution to Eq. (25) accomplishes automatically the conditions (20) and (22) for any $(i,j)\in I$. However, the only immediate condition accomplished by $\theta_{i,j}$ is that $\theta_{i,j}\leq 1$. In fact, we have

[TABLE]

where $f_{i,j}+S\Delta x_{1}h_{i-1,j}\theta_{i-1,j}\geq 0$ and

[TABLE]

with $q$ a constant not depending on $\bm{p}$ nor $\bm{\theta}$. Thus, $\theta_{i,j}$ may be negative in the cells of $I_{\textnormal{r}}$ for which some of its neighbors belongs to the set $I_{0}$. This issue can be overcome by defining $T$ in such a way that $\Delta^{\prime\prime}_{i,j}\geq 0$.

The discrete version of the operator $T=T_{\epsilon}$ is thus defined as

[TABLE]

Where $I^{\epsilon}$, $\epsilon\in\{0,1\}$, corresponds to the indexes $\left(i,j\right)$ of the cells contained in $\Omega_{+}$ such that some of its neighbor cells $\{(i\pm 1,j),(i,j\pm 1)\})$ belong to $\Omega^{\textnormal{r}}_{0}$. By extending the value of $T$ to the neighboring cells the zero normal derivative at rupture points is automatically imposed and, furthermore, the discrete values of $\theta$ become automatically non-negative. Figure 4 illustrates the action of the discrete operator $T_{\epsilon}$. In [33] an algorithm to solve Eq. (25) for the case $T_{\epsilon}\equiv 0$ is proposed. That numerical procedure is presented in Algorithm 1, which is also based on an adaptation of the numerical strategy proposed by Ausas et al [35] to solve the coupling of Eq. (25) with a Newmark scheme for the Newton’s equation (4). The computation of $T(\bm{\theta})$ required at each iteration is performed by a flooding algorithm that starts from the right boundary (the combustion chamber side).

4 Numerical examples

The nonlinearity of the proposed model is quite strong. At each time, the rightmost connected component of the cavitated/separated region must be identified to compute the operator $T$. Further, as a way to impose the zero normal pressure gradient at rupture points, $T$ is extended by one cell outside the numerical boundary of this connected component. Though the model is consistent with the physics, it may well be ill-posed, which would result in a chaotic behavior of the numerical solution as the mesh and the time step are refined. Convergence analyses in non-trivial cases are thus crucial to assess the model and its numerical implementation. Such analyses, with increasing complexity, are reported in what follows.

For the following simulations the bore radius is set to $B_{\textnormal{\scriptsize{r}}}=41$ (4.2 cm), the ring’s length along $x_{1}$ to $L_{\textnormal{\scriptsize{ring}}}=1$ (1 mm), the ring mass to $m=1.25\times 10^{-5}$ per unit width ($0.03$ kg/m), the applied load to $W^{\textnormal{\scriptsize{a}}}=-1.666\times 10^{-4}$ per unit width ($-40$ N/m) and the Back-pressure factor to $\gamma=0.9$.

4.1 One-dimensional stationary solutions

We begin by illustrating the model’s behavior through a series of stationary solutions assuming different values for $p_{\textnormal{\scriptsize{cc}}}$ (constant in time). The ring-to-liner distance $Z$ is also fixed to unity. The solutions are obtained by means of transient simulations for a time length $T$ big enough such that no variation on $\bm{p}$ nor $\bm{\theta}$ is observed for $t>T$. Based in a convergence analysis exposed in the next section, the cell size and time step length are chosen as $\Delta x_{1}=0.005$ and $\Delta t=0.01$.

The pressure and saturation fields obtained for different values of $p_{\textnormal{\scriptsize{cc}}}$ are shown in Fig. 5 for the lower surface moving to the right ($u=1$), thus the point $\beta\in\partial\Omega_{0}^{r}$ corresponds to a rupture point. As indicated in that figure, the higher the value of $p_{\textnormal{\scriptsize{cc}}}$ the lower the fluid film thickness that exits the domain. Moreover, there exists an upper limit on $p_{\textnormal{\scriptsize{cc}}}$ above which the algorithm does not converge. This coincides with the value of $p_{\textnormal{\scriptsize{cc}}}$ for which the rupture point reaches $x_{1}=0.5$, which is around $120$ atm. For higher values of $p_{\textnormal{\scriptsize{cc}}}$ the simulation fails in the sense that a maximum number of iterations (=100,000) is reached without the numerical variable “change” becoming lower than “tol” in Algorithm 1. Since stationary solutions with a rupture point placed at the left side of $x_{1}=0.5$ have no physical meaning, that value can be interpreted as the maximum CCP such that the ring is able to seal the combustion chamber and prevent blow-by. Based on these observations, the condition

[TABLE]

is identified as a blow-by prediction condition at rupture points.

The pressure and saturation fields obtained for different values of $p_{\textnormal{\scriptsize{cc}}}$ are shown in Fig. 6 for the lower surface moving to the left ($u=-1$) so that the cavitation boundary adjacent to the combustion chamber is a reformation point. Observe that for $p_{\textnormal{\scriptsize{cc}}}>60$ atm the pressure profiles coincide to the left of $\beta$ (that depends on $p_{\textnormal{\scriptsize{cc}}}$), only changing to the right of that point (where $p=p_{\textnormal{\scriptsize{cc}}}$). For these cases the pressure gradient is not necessarily null (recall that $\theta_{-}$ can be lower than one in (16)), thus, from the observed behavior the condition

[TABLE]

is identified as a blow-by prediction condition at reformation points.

4.2 Rings with textures

Consider now simulations with a texture consisting of elliptically shaped dimples on the ring. The dimples are arranged at a distance equal to 0.1 along the $x_{2}$ axis, its depth is set to 1 micron and their size along the $x_{1}$ and $x_{2}$ axis is fixed to 80 and 60 microns, respectively. The resulting ring’s profile ($h_{\textnormal{\scriptsize{ring}}}(x_{1},x_{2})$) is detailed in Fig. 7. This is a sufficiently complex case to numerically assess the well-posedness of the proposed model.

4.2.1 Convergence in space

The radial dynamics of the ring is solved as exposed in Algorithm 1, with an initial value $Z(t=0)=0.8$. The cell sizes are set by taking $\Delta x_{1}=\Delta x_{2}$ and by comparing stationary solutions obtained when fixing $p_{\textnormal{\scriptsize{cc}}}$ to 50 atm. These stationary solutions consist of time-converged transient simulations for a period of 0.001 s, long enough for the time derivatives to be negligible. The time step is set such that the Courant-Friedrichs-Levy (CFL) number is equal to the unit.

Figure 8 shows pressure profiles obtained along the centerline of the computational domain (i.e., fixing $x_{2}=0.05$). Spatial convergence of the proposed algorithm is observed as the mesh is refined. Next, time convergence is addressed keeping the number of cells along the $x_{1}$-axis fixed to 200.

4.2.2 Convergence in time and comparison with classical cavitation models

In this transient simulations the CCP is modeled as a Gaussian pulse in time of duration $\simeq 0.02$ s with its amplitude denoted by $A_{\textnormal{\scriptsize{cc}}}$. In Fig. 9 it is shown that the Gaussian pulse indeed approximates a measured CCP (taken from [36]) for $A_{\textnormal{\scriptsize{cc}}}=50$ atm. This example is quite challenging, not just because of the texture, but also because $p_{\textnormal{\scriptsize{cc}}}(t)$ is non-monotonous, making $\Omega_{0}^{r}$ to grow and shrink dynamically.

A convergence analysis was performed by computing the relative difference on friction in time according to

[TABLE]

where $F^{\Delta}$ is a discretization of (9) and $\|\cdot\|_{1}$ is a discretization of $\int_{0}^{600}|\cdot|\,dt$. Figure 10 shows $\Delta\bar{F}^{\Delta t}$ vs. $\Delta t$, considering both $u=-1$ (reformation boundary on the combustion-chamber side of the ring) and $u=1$ (rupture boundary). A convergence rate $\approx\Delta t^{0.7}$ is observed. An analogous analysis on the MFT dependence in $\Delta t$ showed a convergence rate $\approx\Delta t^{0.6}$. These results suggest that the model is indeed well-posed and stably approximated, and it is thus interesting to compare it with the models previously available in the literature, that either consider $p_{\textnormal{\scriptsize{cc}}}=0$ for the hydrodynamic problem (standard Elrod-Adams model) or are not mass-conserving (Reynolds cavitation model). For this comparison, the number of cells along $x_{1}$ is set to 200 and the time step is set to $2\Delta x$.

The comparative results for friction and MFT are shown in Fig. 11. Notice that the simulated time ($0.06$ s) corresponds to a whole engine cycle running at 2000 rpm. To analyze these results the extended model is set as a reference. When comparing with the Elrod-Adams model, a difference of $\simeq 0.5$ microns is found for the MFT while the explosion is taking place, which translates into a 17% of difference in the time-averaged friction. On the other hand, the Reynolds model predicts a much higher value of the MFT ($Z\simeq 6.8$) when $p_{\textnormal{\scriptsize{cc}}}$ is negligible, while for $280<t<330$ the MFT is quite similar to the one predicted by the extended model. Overall, this implies a difference of 25% in the time-averaged friction. The pressure profiles along $x_{2}=\pi B_{\textnormal{\scriptsize{r}}}$ are shown in Fig. 12 for different time instants. While the Elrod-Adams model is not able to impose $p=p_{\textnormal{\scriptsize{cc}}}$ at $x_{1}=1$, the Reynolds model fails to predict the cavitation boundary correctly. Both these issues are resolved by the extended model, which provides a physically sound solution.

Remark: The friction formula used for the Reynolds model is similar to (9) but replacing the Couette term $\mu u/h\,g(\theta)$ by

[TABLE]

4.3 Simulation of the full engine cycle considering wear

In this section both surfaces are assumed untextured, but a non-uniform wear is incorporated on the ring by taking $\delta\geq 0$ in Eq. (1). Notice that the wear can also be interpreted as taking place on the cylinder bore. Based on a convergence analysis analogous to the one presented in the previous section the cell size is set to $\Delta x_{1}=1/200$ and $\Delta x_{2}=2\pi B_{\textnormal{\scriptsize{r}}}/40$, while the Courant number is set to $1$ during most of the simulation. The film thickness at the fluid entrance is set to $h_{\textnormal{\scriptsize{feed}}}=1.5$. The liner’s speed $u$ (relative to the ring) and the normalized variation of $p_{\textnormal{\scriptsize{cc}}}(t)$ are shown in Fig. 13. Some convergence difficulties that arise when $p_{\textnormal{\scriptsize{cc}}}(t)$ is near its maximum and the speed of the ring is small (at the half of the cycle) were overcome by setting the Courant number to $0.05$ for $290<t<314$.

4.3.1 Blow-by inception for different ring’s curvatures

Here the height of the ring’s profile at the boundaries is denoted by $\Delta h\overset{\scriptsize\mbox{def}}{=}h^{R}_{\textnormal{\scriptsize{ring}}}(x_{1}=1)$. A series of simulations are performed by varying $\Delta h$ as $0.5,1.0,1.5,2,\ldots,8.0$ and $\delta$ as $0.005,0.010,0.015,\ldots,0.060$. From the results it is observed that for each $\Delta h$ there exists a value of $\delta$, denoted by $\delta^{\textnormal{\scriptsize{max}}}$, for which the simulation fails at a certain time step before the whole cycle is completed, and the same is observed for every $\delta\geq\delta^{\textnormal{\scriptsize{max}}}$, while for $\delta<\delta^{\textnormal{\scriptsize{max}}}$ the simulation is completed. These results on $\delta_{\textnormal{\scriptsize{max}}}$ are summarized in Table 3. It is worth to notice that $\delta_{\textnormal{\scriptsize{max}}}$ is maximum (as a function of $\Delta h$) for $\Delta h$ in the interval $2.5-3.0$, or $R$ in $41.6-50$, which corresponds to typical values for the compression ring’s curvature [5].

The time step at which the simulations fail varies with $\Delta h$ and $\delta$. In general, for $\Delta h\geq 3.91$ all the failures take place at $t\approx 300$, where $p_{\textnormal{\scriptsize{cc}}}$ is maximum, while for $\Delta h<3.91$ the time at which the simulations fail diminishes and the value of $p_{\textnormal{\scriptsize{cc}}}$ at that instant is much lower. Next, we give more details on this, showing the pressure profiles around the failure time for different regimes.

Setting first $\Delta h=2$ ($R=62.5$) and the wear amplitude to $\delta=0.040$ one observes that the simulation is completed. Five time instants ($t_{\textnormal{\scriptsize{A}}}$, $t_{\textnormal{\scriptsize{B}}}$, …) are chosen as indicated in Fig. 13 to show the pressure and saturation fields in Fig. 14. Notice that the full-film region (in black in the upper part of the figures for each time) separates two incomplete-film regions, one on the top that is connected to the combustion chamber, and one on the bottom where the pressure value is equal to 1 atm. The effect of the wear appears as a thinning of the full-film region around the circumferential position $x_{2}=\pi B_{\textnormal{r}}$, which corresponds to the sector of maximum wear.

Increasing $\delta$ to $0.045$ for the same ring curvature ($R=62.5$), one observes that as $t$ tends to $t=203$ the full-film region under the ring tends to break up, as shown in Fig. 15 . After $t=203.18$ the Algorithm 1 fails to converge. At that simulation time, the CCP is equal to $0.035$ atm and our interpretation of the algorithm’s failure is that the proposed model looses validity because blow-by takes place. Details on the pressure profiles and film-thickness along $x_{2}=0.5$ are given in Fig. 16. Notice the drop of the pressure build-up in the full-film region $\{\theta=1\}$ as the value of $p_{\textnormal{\scriptsize{cc}}}$ augments in time.

When $\Delta h$ is augmented to $3.91$ ($R=32$) the same phenomenon takes place even for $\delta$ as small as 0.005, while for $\Delta h=4$ ($R=31.25$) the simulations fails for every $\delta\geq 0$ tested. For this last case some of the pressure profiles and film-thickness along $x_{2}=0.5$ in time are detailed in Fig. 17. As for the previous case, the pressure build-up in the full-film region decreases and the rupture points (for $t>300$) advance to the left.

4.3.2 A blow-by inception criterion

In this section we focus on simulations that ended because of lack of convergence of Algorithm 1 (defining a maximum number of iterations equal to $100,000$). We show that this lack of convergence can be related to blow-by inception as identified in the previous section. Further, and encompassing condition (27) (for rupture boundary) and condition (28) (for reformation boundary) which are hard to test numerically, we propose the condition

[TABLE]

where $\Omega_{<}$ is the region of the domain where $p<p_{\textnormal{\scriptsize{cc}}}$, while $d\left(\Omega_{0}^{\textnormal{\scriptsize{r}}},\Omega_{<}\right)$ denotes the distance between the set $\Omega_{0}^{\textnormal{\scriptsize{r}}}$ (where $p=p_{\textnormal{\scriptsize{cc}}}$) and $\Omega_{<}$. Essentially, this condition expresses that if the full-film region separating the gas region in which $p=p_{\textnormal{\scriptsize{cc}}}$ from the region where $p<p_{\textnormal{\scriptsize{cc}}}$ is narrow enough, blow-by will happen. In this work we take $\epsilon_{\textnormal{\scriptsize{b}}}=1/50=0.02$ and $N_{\textnormal{\scriptsize{b}}}=4$, but the results are not much sensitive to this specific choice. In fact, a convergence analysis showed that $d\left(\Omega_{0}^{\textnormal{\scriptsize{r}}},\Omega_{<}\right)$, as a function of time, becomes mesh-independent once the mesh is fine enough. Figure 18 shows $d\left(\Omega_{0}^{\textnormal{\scriptsize{r}}},\Omega_{<}\right)$ for different configurations. For the ones where the simulation failed, it is observed that the condition (29) is reached some time steps before the end of the simulation. While for the case where the whole cycle is completed the condition (29) is never reached. This qualifies (29) as a criterion to numerically predict blow-by.

5 Conclusions

An extension of the Elrod-Adams model that accommodates non-homogeneous boundary conditions for pressure has been presented and it has been applied to the PRL system. A mass-conserving algorithm that imposes a non-negative pressure gradient condition at the rupture boundary is used, which exhibits convergence in both time and space. To our knowledge, these are the first mass-conserving transient simulations for the PRL mechanism that predict a rupture boundary with a non-zero value of $p_{\textnormal{\scriptsize{cc}}}$ enforced at it.

Full-engine cycle simulations were also presented, which showed significant differences in the friction computed by means of the extended model when compared to the one obtained with either the standard Elrod-Adams or the Reynolds models, none of which is able to correctly predict the position of the pressurized regions.

A perturbation on the ring’s surface along the circumferential direction was assumed to model a non-uniform ring wear. For a big-enough perturbation the separated region where the pressure value is equal to the CCP (or $p_{\textnormal{\scriptsize{cc}}}(t)$) reaches the separated region where the pressure is zero, which is interpreted as leading to blow-by. Based on observations of stationary pressure profiles, two analytic conditions on the derivatives of $p$ at the rupture/reformation points were proposed to predict blow-by inception. Furthermore, an equivalent and easily implementable numerical criterion (observed to be mesh-independent) was proposed, and it was shown that users can rely on it to predict blow-by from numerical simulations of PRL systems.

Acknowledgments

The authors thank the financial support provided by CAPES (grant PROEX-8434433/D), FAPESP (grants 2013/07375-0 and 2018/08752-5) and CNPq (grant 305599/2017-8).