## Abstract

This research focuses on the multiphase oil film tribology between the piston pin and the connecting rod in an internal combustion engine and establishes a new computational approach for thin-film lubrication with unsteady flow channel variation. First, the pin and the connecting rod are considered as rigid bodies, and 3D numerical analysis of the cavitating lubricating oil flow is performed when combustion load is applied to the pin. We find that dynamic pressure does not increase around the connecting rod edge and that pressure is potentially insufficient to support the load. In the second numerical analysis, the pin and the connecting rod are considered to be elastically deformable structures, and coupled 3D multiphase fluid–structure interaction simulation is performed. The boundary lubrication area is detected using a statistical Greenwood–Tripp model as unevenness of the contacted metal surface. The results show that pressure distribution spreads more widely than in the result for rigid bodies and that the film was thicker as well. Also, the pin deformed like a bow, but the deformation of the connecting rod was quite small, suggesting a potential mechanical contact at the edge of the connecting rod with the pin. By comparison with an actual operationally used piston pin, we find that the fluid–structure coupled analysis qualitatively predicted the seizure location.

## 1 Introduction

There is a growing need for further improvements to the efficiency of internal combustion engines due to their impact on environmental and sustainability problems [1]. In internal combustion engines for automobiles, mechanical energy loss accounts for 10–20% of the energy produced, so there are increasing efforts toward its reduction [2]. Furthermore, about 40% of mechanical loss is due to friction in the piston system, so it is important to clarify tribological phenomena there [3]. Reducing the weight and size of parts and reducing the viscosity of lubricating oil effectively realizes low friction, but at the same time, there is an increased risk of internal combustion engine failure due to contact between parts.

In a reciprocating engine, power generated by combustion is extracted as reciprocating motion and converted to rotary motion for driving force [4]. Reciprocating engines have many uses but are primarily used to power automobiles. It is important to efficiently transmit a large amount of power from combustion to the drive system, but high durability is required for parts, which must operate under extremely harsh conditions. The most common cause of failure in a reciprocating engine is seizure of its parts, a phenomenon that occurs when there is a break in the lubricating, thus allowing metal parts to come into contact and causing damage or adhesion (Fig. 1). Seizure generally makes engines inoperative, often causing such extensive damage that the entire engine must be replaced. Lubricating oil plays a major role in preventing seizure, so reciprocating engines must be designed to maintain an oil film between parts.

In a reciprocating engine, reciprocating motion of a piston is converted into rotational motion of a crankshaft through use of a crank mechanism [5]. The rod connecting the piston and the crankshaft is called the connecting rod, and its structure is divided into a small end connected to the piston side and a large end connected to the crankshaft side. Both ends are connected by pins, and relative rotational movement occurs at a sliding part. In particular, the small end is connected by a pin called the piston pin, which is one of the most important lubrication points in a reciprocating engine. The small end supports the explosive pressure generated in the combustion chamber, transmitting that force to the connecting rod and crankshaft. In contrast to the large end, the relative motion of the small end is not rotational motion that keeps sliding in the same direction but oscillating motion that changes the direction of rotation at regular intervals. In other words, at the moment when the direction of rotation switches, the relative velocity becomes zero, and there is a very high possibility that metal parts will come into contact, causing seizure. From the above, lubrication between the piston pin and connecting rod is bearing lubrication under the harshest conditions in an internal combustion engine.

Since the piston system repeats up–down reciprocating motion at high speed, large inertial forces are applied. Weight reductions in the piston system are thus indispensable for improving fuel efficiency and durability. The piston pin is generally made of iron so that it can withstand deformation due to explosive pressure, so its mass is disproportionately large compared with that of the entire piston system. Currently, the piston pins used in reciprocating engines are mostly what is called the fully floating type, in which there is clearance between the piston pin and the connecting rod. Filling this clearance with lubricating oil prevents the seizure of these parts. The piston pin has a hollow cylindrical shape for weight reduction, but its inner diameter is optimized to minimize deformation even under heavy loads. Automakers and others are making various efforts toward reducing the weight of piston pins. Measures to reduce their diameter, for example, include distributing the applied load.

As described above, there is increasing demand for internal combustion engines with higher efficiency, and research on piston system lubrication for reducing mechanical energy loss is actively conducted. It is difficult to visualize and quantitatively analyze lubrication phenomena in experiments, so it is essential to establish numerical analysis methods that can be applied on computers. There have been some interesting studies involving both simulations and experimental methods related to piston pin lubrication problems. Allmaier and Sander [6] found that most piston pins rotate at very slow rotational speeds and even change rotational directions according to operating conditions. They also investigated influencing effects on this dynamic behavior such as clearance and pin surface roughness to see how the lubrication of this crucial part could be improved. Fang et al. [7] developed a numerical model to analyze lubricated full-floating pin bearings in planar multibody systems, which couple a mixed lubrication model with multibody motion equations. To consider the effects of surface roughness and asperity contact, they integrated the average Reynolds equation and the GT asperity contact model to establish the mixed lubrication model. Ferretti et al. [8] performed a tribological study showing that two different asperity contact models have been implemented, the former based on standard Greenwood–Tripp theory and the latter based on a complementarity formulation of the asperity contact problem, and compared the results. They showed that both approaches can provide similar results if the model input data are calibrated.

There have also been many papers on piston pin and connecting rod tribology, focusing mainly on mechanisms and multibody dynamics [7], but few studies have focused on multiphase flow behavior inside narrow thin-film channels, which is the subject of this research.

Even in studies of the state between connecting rods and the small end of piston pins, analyses commonly use the Reynolds equation [9,10], which is a simplification of the Navier–Stokes equation, to perform numerically stable analyses. Furthermore, the large loads applied to the piston pin deform the piston pin–connecting rod structure. Conventional lubrication theory underestimates the gap width because it assumes that the structure is a rigid body, so it is necessary to perform analyses based on elastohydrodynamic lubrication (EHL) theory, which treats the structure as an elastic body [11–16]. In general EHL theory, however, only one structure is treated as an elastic body [17,18], so there are limits to considering deformation of both the piston pin and connecting rod, which is the subject of this research.

Ozasa et al. [19] used EHL theory to perform calculations for the large end of the connecting rod, which experiences less harsh lubrication conditions compared with those of small end, confirming the validity of calculation results regarding film thickness and axial locus by comparison with experimental results from an actual machine. Lin et al. [20] performed numerical analyses of lubrication between the small end of the connecting rod and the piston pin. They used the Reynolds equation as the governing equation for lubricating oil and attempted to elucidate seizure phenomena of the piston pin by considering structural deformation. These numerical analyses demonstrated that deformation of the piston pin could cause seizure with the small end of the connecting rod, but satisfactory results were not obtained for lubricating oil flow, so this work serves only for discussion.

Thermal effects in the lubricant are required to formulate the governing equations because temperature distribution due to frictional heat is formed in the thin-film layer where unsteady fluctuating pressure is applied. Also, cavitation occurs with phase change, which influences the behavior and seizure location. In addition, although many previous studies have considered the effects of cavitation in the tribology field [21–27], few studies have evaluated interactions between EHL effects and cavitation.

This research focuses on the multiphase narrow gap flow between the small end of piston pins and connecting rods with the unsteady variation of channel width. Furthermore, it established a numerical method based on the fluid–structure interaction (FSI) model that considers the momentum and energy exchange at the piston pin-lubricant interface and the exchange at the conrod-lubricant interface. To that end, we develop a 3D multiphase fluid–structure coupled analysis that considers elastic deformation of the structure and changes in the flow channel, along with an analysis method that reproduces the complex phenomena of lubricating oil with cavitation under harsh conditions. Our objective is to make proposals for optimizing part shapes and predicting seizure points using the newly developed analysis method. This study also focuses on elastic deformation of piston connecting rods due to flow in a narrow gap with phase change and fluid–structure interactions rather than multibody and kinematic investigations such as the motion of piston pin rotation. We also performed a coupled multiphase fluid–structure analysis with 3D unsteady narrow-channel deformation, which is difficult to analyze with TEHL theory [28], so a simple comparison with TEHL theory is difficult. We perform a comparison of rigid and elastic models to demonstrate the effectiveness of this analysis.

## 2 Rigid-Body Analysis of Lubricating Oil Flow Between the Piston Pin and the Small End of the Connecting Rod

We first perform a coupled analysis using rigid-body approximation, then a comparative study to clarify the effect of elastic deformation on the EHL characteristics. The numerical analyses in this section are used to perform composite fluid analysis of the lubricating oil and rigid bodies, thereby reproducing flow between the piston pin and connecting rod. These numerical analyses consider the connecting rod as a rigid body and apply load to reproduce a flow channel that unsteadily changes. This section describes the fluid model used in the analyses, the mesh morphing method, the analysis conditions, and the analytical results. To consider the numerical analysis results, we also describe an analytical method for pressure distributions in the bearing, considering that lubrication between the piston pin and connecting rod can be simulated by the bearing lubrication model.

### 2.1 Basic Fluid Phase Equations and Numerical Method.

*ρ*is the density, $U$ is the flow velocity vector,

*p*is the pressure, and $\tau $ is the viscous stress tensor represented as

*μ*is the viscosity and

**I**is the unit tensor.

### 2.2 Cavitation Model.

*v*and

*l*subscripts indicate vapor and liquid phases, respectively, and the volume fraction of vapor

*α*is

### 2.3 Lubrication Compressibility.

*ρ*

_{l,0}is the liquid phase density under atmospheric pressure, and coefficients C

_{1}= 5.9 × 10

^{8}(Pa) and $C2=1.34(\u2212)$.

### 2.4 Rheology Model.

Viscosity of the lubricating oil, the target of this numerical analysis, is known to change under conditions where the pressure becomes extremely high. This is called the piezo-viscosity effect, and many studies have aimed to predict the correct viscosity under high shear-rate conditions. As the load increases and pressure in the lubricating oil rises, the density changes due to compressibility, decreasing the kinematic viscosity coefficient. However, as pressures increase further, the viscosity coefficient of the lubricating oil sharply increases due to the piezo-viscosity effect, significantly affecting the change in kinematic viscosity. To predict the viscosity coefficient of lubricating oil under extremely high shear-rate conditions, it is therefore necessary to use an appropriate rheological model [14,15,35,41,42].

*η*

_{0}is the viscosity coefficient under atmospheric pressure, and

*α*

_{pv}is the pressure–viscosity coefficient, which is a characteristic value specific to the fluid. Lubricating oils are known to show non-Newtonian properties under high shear-rate conditions. In this numerical analysis, we use the Eyring model, usefulness of which has been confirmed in EHL research

*η*

_{Barus}is the fluid viscosity as calculated from Eq. (8), $\gamma \u02d9$ is the shear rate,

*τ*is the shear stress, and

*τ*

_{0}is the Eyring stress. Shear stress

*τ*is calculated as

^{−8}1/s we use Eq. (11) to find the viscosity of the lubricating oil, and Eq. (8) otherwise.

### 2.5 Mesh Morphing Model.

*γ*is the displacement diffusion coefficient and

*d*is the amount of mesh displacement [43]. Displacement diffusivity is a parameter that determines the mesh distribution in the calculation area and is a function of the distance

*r*from the rigid-body interface, given by

### 2.6 Analytical Method for Bearing Lubrication.

In this computation, numerical constraints necessitate adding boundary conditions for oil inflow. An approximate analytical solution can be obtained for such bearing lubrication. Therefore, we obtain an analytical solution by applying the Sommerfeld condition [44] to the Reynolds equation [10,45], which is a simplification of the Navier–Stokes equation.

### 2.7 Half-Sommerfeld Condition.

In the Sommerfeld condition solution shown above, the pressure is negative below atmospheric pressure in the range where the gap width becomes large. In that negative pressure region, the lubricating oil falls below the saturated vapor pressure and evaporates, causing cavitation, so the assumption that the bearing gap is filled with only lubricating oil does not hold under these conditions.

The various available cavitation models are extensively compared in Ref. [46], and both the Reynolds (Swift-Stieber) condition and the Jakobsson–Floberg–Olsson (JFO) cavitation model [47–49] have been used in conventional research to investigate the performance of laser-textured piston rings or sliders [50–52].

However, it still needs to be experimentally validated if multiphase cavity flow with narrow channel deformation can be accurately modeled with a thin-film approximation (e.g., the Reynolds equation) employing a saturation parameter (e.g., the JFO cavitation model) [53–57].

However, mass conservation is not satisfied, so the load capacity is overestimated in cases where mass flux limits the pressure increase. Still, behavior of the asperity parameters was qualitatively identical to that in the JFO model, resulting in the same optimal asperity parameters [58].

The half-Sommerfeld condition proposed by Gümbel [59–61] is a model that considers the occurrence of such cavitation. In this model, the negative pressure part in the pressure distribution obtained under the Sommerfeld condition is set to 0 pressure. For these reasons, we applied the half-Sommerfeld condition in this study.

### 2.8 Analysis Model and Numerical Analysis Conditions.

This section describes the analysis model and numerical analysis conditions used for flow analysis of lubricating oil in the narrow flow channel between the piston pin and connecting rod. We performed this numerical analysis and computing using *cavitatingTribologyFoam*, a customized cavitation solver based on our previously developed solver cavitatingLesInterFoam [33] in the open-source computational fluid dynamics (CFD) code openfoam [62].

Lubricating oil flows in from an oil reservoir at the top of the connecting rod, filling the gap between the piston pin and the connecting rod. In this analysis, therefore, we extract the gap region between the piston pin and the connecting rod and create a fluid mesh. Figure 4 shows the calculation region and the mesh used for analysis, with the *z*-axis showing the axis of rotation, the *x*-axis showing the horizontal radial direction, and the *y*-axis showing the vertical radial direction. Figure 5 shows dimensions in the numerical model.

Both the piston pin and the connecting rod are treated as rigid bodies. The initial clearance width between the two parts is extremely narrow at 11.75 *μ*m, but we created a mesh with seven layers to obtain a sufficiently accurate velocity distribution between the piston pin and connecting rod (Fig. 6). The number of mesh cells is approximately 200,000. We apply a load in the *y*-axis direction and treat the piston pin as a moving rigid body. In an actual reciprocating engine, such as a four-stroke engine, each cycle comprises four steps: intake, compression, combustion, and exhaust. The total inertial force due to vertical reciprocating motion and explosive pressure received from the combustion chamber is applied to the piston pin. Figure 7 shows the load history on a piston pin over the four cycles of a reciprocating engine. Based on actual engine operating conditions, the load history on the piston pin data and the connecting rod velocity data were obtained by flexible multibody dynamic analysis using RecurDyn [63,64]. The highest load of 4672 N occurs during the combustion process. To reduce calculation times in this analysis, instead of unsteady analysis through four cycles, we perform calculations by extraction from the middle of the compression process to the middle of the exhaust process, where the possibility of seizure is the highest. Figure 8 shows the load history applied to the piston pin by numerical analysis.

Table 1 and Fig. 9 show the boundary conditions used for numerical analysis. Both the upper oil reservoir and the gap between the piston pin and connecting rod are set to free inflow and outflow conditions, and the pressure boundary is set to atmospheric pressure. In other words, lubricating oil flows in or out depending on the pressure difference between the inside and outside of the calculation region. Both the piston pin and the connecting rod wall surface have no-slip boundary conditions. Oscillatory movement of the connecting rod provides velocity according to the swing angle shown in Fig. 10 as the velocity boundary condition on the piston pin wall surface.

Boundary | Pressure | Velocity |
---|---|---|

Inlet | Constant atmospheric (101.325 kPa) | Continuous in–out flow |

Outlet | Constant atmospheric (101.325 kPa) | Continuous in–out flow |

Piston pin | Neumann condition | No slip |

Connecting rod | Neumann condition | No slip |

Boundary | Pressure | Velocity |
---|---|---|

Inlet | Constant atmospheric (101.325 kPa) | Continuous in–out flow |

Outlet | Constant atmospheric (101.325 kPa) | Continuous in–out flow |

Piston pin | Neumann condition | No slip |

Connecting rod | Neumann condition | No slip |

Since the gap between the piston pin and connecting rod changes unsteadily, we apply an adaptive step and perform the unsteady calculation so that the Courant number is one or less. We also assume that the engine is operating at $100\u2218C$, with kinematic viscosity *ν* = 1.0 × 10^{−5} m^{2}/s, and density *ρ* = 823 kg/m^{3} for the lubricating oil. In addition, since the Reynolds number Re is assumed to be about 100 at maximum, the flow is assumed to be laminar. Since this study deals with multiphase flow with cavitation and a free surface based on phase change, a high-precision CFD model is necessary to accurately capture the interfacial flow in a narrow-channel thin-liquid film. Table 2 shows physical characteristics of the lubricating oil used in the analysis. We assume SAE-viscosity grade 10W–30 [65] as physical properties of the oil.

Lubricant | Vapor | Unit | |
---|---|---|---|

Viscosity | 1.0 × 10^{−5} | 8.97 × 10^{−6} | (m^{2}/s) |

Density | 850 | 0.029 | (kg/m^{3}) |

Compressibility | 4.9 × 10^{−7} | 5.7 × 10^{−6} | (1/kPa) |

Saturated vapor pressure | 5000 | – | (Pa) |

Pressure–viscosity index | 0.689 | – | (–) |

Eyring stress | 5.0 × 10^{6} | – | (Pa) |

Lubricant | Vapor | Unit | |
---|---|---|---|

Viscosity | 1.0 × 10^{−5} | 8.97 × 10^{−6} | (m^{2}/s) |

Density | 850 | 0.029 | (kg/m^{3}) |

Compressibility | 4.9 × 10^{−7} | 5.7 × 10^{−6} | (1/kPa) |

Saturated vapor pressure | 5000 | – | (Pa) |

Pressure–viscosity index | 0.689 | – | (–) |

Eyring stress | 5.0 × 10^{6} | – | (Pa) |

### 2.9 Numerical Analysis Calculations and Discussion.

The following describes the results of flow analysis of the lubricating oil between the piston pin and connecting rod. Figures 11–16 respectively show pressure distribution contours (left) and graphs in which the perimeter of the *z* = 0 vertical cross section is plotted on the horizontal axis and the pressure is plotted on the vertical axis (right). Additionally, Figs. 14–16 respectively show the lubricating oil velocity vector distribution (left) and the volume fraction distribution of vapor phase due to cavitation (right) at each time.

In Fig. 11, we can see pressure rising (to a maximum of approximately 650 MPa) in the direction of load application, that is, at the bottom of the analysis model. We adjusted maximum and minimum values for the color legend range to more clearly visualize the pressure distribution in the region (left figure). We can also confirm that because the piston pin moves rapidly, large negative pressure is generated in the upper part of the model, and cavitation occurs. The load applied to the piston pin increases up to *t* = 1.5 × 10^{−3} s, and the wall surface of the piston pin and connecting rod rapidly approach, causing the squeeze effect [66] in the Reynolds equation to appear prominently, increasing pressure in the lubricating oil. The increase in the pressure magnitude is also due to the decreased local film thickness caused by the load increase. We can also see that as the load increases, the pressure distribution increases locally in the central part of the model. This is likely caused by a local increase in the viscosity of the lubricating oil due to the piezo-viscosity effect as modeled by the Barus equation, in which viscosity exponentially increases as the pressure increases. This pressure increase allows the lubricating oil to support the load applied to the piston pin. However, the high-pressure region in the lubricating oil is limited to the center of the lower part of the model, and pressures around the edge are near atmospheric pressure, so no large pressure increases are observed there. As Fig. 14 shows, this occurs because lubricating oil is pushed out by the sudden movement of the piston pin. We can therefore expect the piston pin, which is more easily deformed than the connecting rod, to undergo a bow-like deformation.

Figures 12 and 13 show the analysis results for pressure distributions when load on the piston pin reaches the peak value of 4672 N, and the load gradually decreases. As the load decreases and the movement speed of the piston pin in the *y* direction decreases, the pressure increase effect of the lubricating oil due to the squeeze effect decreases, and the local pressure increase dulls the sharply pointed distribution. Because the eccentricity nearly constantly rotates, the wedge effect becomes stronger, and a pressure distribution similar to that under the half-Sommerfeld condition arises at *t* = 2.5 × 10^{−3} s. Figures 15 and 16 show that the amount of lubricating oil pushed out was small and that the oscillatory motion of the connecting rod formed a velocity field that rotated in the circumferential direction. We can also confirm that the cavitation region occurring in the upper part shrank with the load phenomenon and that cavitation sparsely occurred in the lower part of the model.

Figure 17 shows the initial state (at *t* = 0 s) and the gap width between the piston pin and connecting rod at *t* = 1.5 × 10^{−3} s, where the flow channel was narrowest. The flow channel width was 11.75 *μ*m at the initial state but became as narrow as 1.8 *μ*m due to negative displacement of the piston pin in the *y*-direction under loading. The aspect ratio was approximately 50 at the initial state and increased to approximately 290 at *t* = 1.5 × 10^{−3} s, but we confirmed the possibility of calculations that maintain numerical stability.

## 3 Coupled Fluid–Structure Interaction Analysis of Lubricating Oil Flow Between the Piston and Small End of the Connecting Rod

In the previous section, we performed a rigid-body analysis on lubrication flow between a piston pin and the small end of a connecting rod to predict the location where seizure will occur if the structure does not deform. In Fig. 3, we also showed the validity of the current numerical analysis method through qualitative comparison with analytical solutions by the Reynolds equation. In this section, we extend the solver used in the previous section to perform numerical analysis that considers elastic deformation of the structure. Namely, we conduct coupled FSI analysis [67–69] of lubricating oil flow between the piston and the small end of the connecting rod. By adding an energy equation and implementing a temperature coupling between the structure and the fluid [70], we also consider the thermal effects, which is an essential factor in seizure phenomena [36,44,71–73]. Finally, using the results obtained from this solver, we propose a method for determining piston pin shapes that satisfy the required balance between seizure resistance and weight reduction.

### 3.1 Governing Equation and Numerical Method in the Fluid Model.

*h*is the enthalpy,

*k*is the thermal conductivity,

*T*is the temperature. The second and third terms on the right side of Eq. (16) respectively represent frictional heat (energy transport due to viscous dissipation) and compressive heat. These generated heats are transferred to the solid side via conjugate heat transfer (CHT) [74], resulting in thermal deformation of the solid.

#### 3.1.1 Rheology Model Considering the Thermal Effect.

*η*

_{Barus}is the viscosity according to the Barus equation,

*T*is the temperature,

*T*

_{0}is the basis temperature, and

*β** is defined as

*η*

_{0}is the viscosity under atmospheric pressure,

*Z*is defined as

_{0}is defined as

*α*

_{pv}is the pressure–viscosity coefficient, and

*β*is the temperature–viscosity coefficient, which are characteristic values specific to the fluid.

### 3.2 Governing Equation in the Structural Model.

*s*denotes the structure side.

*ρ*

_{s}is the density of the structure, $D$ is the displacement vector, $\sigma $ is the stress tensor, $fs$ is the body force,

*c*

_{v,s}is the specific heat, and

*k*

_{s}is the thermal conductivity. Using the linear elastic stress $\sigma Le$ and thermal stress $\sigma T$, we can find the stress tensor as

*α*

_{s}is the volume expansion rate,

*T*

_{s}is the solid phase temperature, and

*T*

_{0}is the solid phase reference temperature

*μ*

_{s}and

*λ*

_{s}are Lame’s constants

### 3.3 Fluid–Structure Coupling.

*i*denotes the interfacial value. We also respectively define the exchange of stress and heat flux, $q$ at the interface between the fluid and the structure as

#### 3.3.1 Solution by Strongly Coupled Partitioning.

There are two main methods for finding solutions in fluid–structure coupled analysis: the partitioned method [79] and the monolithic method (Fig. 18). The former method calculates fluids and structures using different solvers, while the latter uses a single matrix to simultaneously solve for fluids and structures. The monolithic method is generally considered to provide higher numeric stability, but there are many issues related to its extension to parallel computing with regional segmentation. In this analysis, the total number of mesh cells for the fluid and the structure is about 330,000, the aspect ratio is large, and the time-step cannot be increased, so computational loads are large. Parallel computing is thus indispensable, so this analysis is performed using the partitioned method.

Partitioned solutions can be further subdivided into weakly and strongly coupled solutions. Weakly coupled analysis is a method that uses values at the boundary surface obtained from fluid analysis as the boundary condition on the structure side. However, this method is disadvantageous in that interactions between the fluid and the structure are considered in only one direction. This is inappropriate for weakly coupled analysis of phenomena such as EHL, where the structure is subjected to large forces. In this analysis, therefore, we attempt numerical analysis using strong coupling.

*Res*

_{fsi}for convergence testing of the fluid–structure interface is

*D*

_{fb}and the amount of displacement at the boundary surface of the solid region $Dfbsolid$ are interpolated on the fluid side. We use the obtained boundary surface residual $Resfsin$ and the immediately preceding residual $Resfsi0$ to find the relaxation factor

#### 3.3.2 Boundary Lubrication Model.

When no load is applied to the piston pin, the wall surfaces of the piston pin and connecting rod are sufficiently separated, with the two parts separated by only lubricating oil. As the load increases, however, the gap between the piston pin and the connecting rod narrows. Then, microscale asperities of the metal surface begin to cause mechanical contact; this lubrication state is called boundary lubrication (Fig. 20). Mixed lubrication may also occur in the system, wherein hydrodynamic and asperity contact coexist. Predicting boundary lubrication is vital for predicting seizure in lubrication where large loads are generated, as in this analysis.

#### 3.3.3 The Greenwood–Tripp Model.

*p*

_{c}generated by mechanical contact at the rough surface of the piston pin and connecting rod is calculated as

*K*is given by

*E*

_{G}is given by

_{asp}is the number of asperities per unit area,

*β*

_{asp}is the radius of asperity curvature, $\sigma *$ is calculated as [80,84]

_{1}and σ

_{2}are standard deviations for asperity heights on the opposing metal surfaces. Assuming the height distribution of the asperities follows a Gaussian distribution,

*F*(

*λ*) in Eq. ( 37) can be calculated as [28,85]

*h*

_{mt}is the distance between the two metal surfaces [85]. In other words, when the gap width becomes smaller and

*λ*becomes 4 or less, we can judge that the asperities are in contact and the boundary lubrication state is reached. The current model does not allow analysis of full asperity contact and requires a thin lubricant film area. Numerical methods for conditions of complete asperity contact are currently under investigation.

### 3.4 Analysis Model and Numerical Analysis Conditions.

This section describes the analysis model and analysis conditions used for numerical analyses of EHL for piston pins and connecting rods. To compute the present system, we newly developed a custom solver named *cavitatingFsiEHLFoam* using the open source finite volume CFD code openfoam [86], which is in turn based on an extention of *solids4Foam*, an under-development openfoam solver [87,88]. This analysis requires solving for structural deformation, increasing computational loads over those for the analysis performed in the previous section, so we use a half model to reduce calculation times.

### 3.5 Fluid Analysis Model.

The fluid model used in this numerical analysis is very similar to the one described in the previous section, but as mentioned above, we use the half model (about 100,000 mesh cells) shown in Fig. 21 to reduce calculation times. Pressure and velocity boundary conditions are the same as in the analysis in Sec. 3. Also, this analysis considers temperatures, so we set the initial temperature condition as 350 K and apply the Neumann condition to the entire boundary (Table 3).

Boundary | Pressure | Velocity | Temperature |
---|---|---|---|

Inlet | Constant atmospheric | Continuous in–out flow | Neumann condition |

Outlet | Constant atmospheric | Continuous in–out flow | Neumann condition |

Piston pin | Neumann condition | No slip | Neumann condition |

Connecting rod | Neumann condition | No slip | Neumann condition |

Symmetry plane | Symmetry condition | Symmetry condition | Neumann condition |

Boundary | Pressure | Velocity | Temperature |
---|---|---|---|

Inlet | Constant atmospheric | Continuous in–out flow | Neumann condition |

Outlet | Constant atmospheric | Continuous in–out flow | Neumann condition |

Piston pin | Neumann condition | No slip | Neumann condition |

Connecting rod | Neumann condition | No slip | Neumann condition |

Symmetry plane | Symmetry condition | Symmetry condition | Neumann condition |

Table 4 shows physical property values for the lubricating oil used in the rheology model by Eyring–Houpert. SAE-Viscosity grade 10W–30 [65] was assumed for the physical properties of the oil.

Lubricant | Vapor | Unit | |
---|---|---|---|

Viscosity | 1.0 × 10^{−5} | 8.97 × 10^{−6} | (m^{2}/s) |

Density | 850 | 0.029 | (kg/m^{3}) |

Compressibility | 4.9 × 10^{−7} | 5.7 × 10^{−6} | (1/kPa) |

Saturated vapor pressure | 5000 | – | (Pa) |

Pressure–viscosity index | 0.689 | – | (–) |

Thermo-viscosity index | 0.0476 | – | (–) |

Eyring stress | 5.0 × 10^{6} | – | (Pa) |

Thermal conductivity | 0.15 | 0.025 | [W/(m · K)] |

Heat capacity | 2300 | 1800 | [(J/(kg · K)] |

Lubricant | Vapor | Unit | |
---|---|---|---|

Viscosity | 1.0 × 10^{−5} | 8.97 × 10^{−6} | (m^{2}/s) |

Density | 850 | 0.029 | (kg/m^{3}) |

Compressibility | 4.9 × 10^{−7} | 5.7 × 10^{−6} | (1/kPa) |

Saturated vapor pressure | 5000 | – | (Pa) |

Pressure–viscosity index | 0.689 | – | (–) |

Thermo-viscosity index | 0.0476 | – | (–) |

Eyring stress | 5.0 × 10^{6} | – | (Pa) |

Thermal conductivity | 0.15 | 0.025 | [W/(m · K)] |

Heat capacity | 2300 | 1800 | [(J/(kg · K)] |

### 3.6 Structural Analysis Model.

Figure 22 shows a numerical model of the piston pin–connecting rod structure, and Fig. 5 shows dimensions in the analysis model. There is a total of 224,000 mesh cells. Figure 23 shows the boundary names for both parts. The load applied to the piston pin in the combustion process described in Sec. 3 is applied to the wallTraction boundary, which connects the piston head and the piston pin. Table 5 shows structural displacement and the temperature boundary conditions. Both the piston pin and connecting rod are made of iron, and Table 6 shows material property values for the structure.

Boundary | Displacement | Temperature |
---|---|---|

Connrod | Zero traction | Neumann condition |

Pin | Zero traction | Neumann condition |

ConnrodWalls | Zero traction | Neumann condition |

ConnrodWallBottom | Zero displacement | Neumann condition |

PinWalls | Zero traction | Neumann condition |

ConnrodSym | Symmetry condition | Symmetry condition |

PinSym | Symmetry condition | Symmetry condition |

WallTraction | Dynamic traction | Neumann condition |

Boundary | Displacement | Temperature |
---|---|---|

Connrod | Zero traction | Neumann condition |

Pin | Zero traction | Neumann condition |

ConnrodWalls | Zero traction | Neumann condition |

ConnrodWallBottom | Zero displacement | Neumann condition |

PinWalls | Zero traction | Neumann condition |

ConnrodSym | Symmetry condition | Symmetry condition |

PinSym | Symmetry condition | Symmetry condition |

WallTraction | Dynamic traction | Neumann condition |

Property | Value | Unit |
---|---|---|

Density | 7850 | (kg/m^{3}) |

Young’s modulus | 209 × 10^{9} | (Pa) |

Poisson’s ratio | 0.3 | (–) |

Thermal conductivity | 50 | (W/(m · K)) |

Thermal elasticity | 1.0 × 10^{−6} | (1/K) |

Heat capacity | 500 | [(J/(kg · K)] |

Property | Value | Unit |
---|---|---|

Density | 7850 | (kg/m^{3}) |

Young’s modulus | 209 × 10^{9} | (Pa) |

Poisson’s ratio | 0.3 | (–) |

Thermal conductivity | 50 | (W/(m · K)) |

Thermal elasticity | 1.0 × 10^{−6} | (1/K) |

Heat capacity | 500 | [(J/(kg · K)] |

Table 7 shows information related to surface roughness in the Greenwood–Tripp model. Values in that table are based on measured data for the surface roughness of piston pin and conrod contact surfaces.

### 3.7 Fluid-Side Analysis Results.

Figures 24–26 show analysis results for the volume fraction of the vapor phase due to cavitation at each of the respective times above. Large negative pressure is generated in the upper part of the analysis region due to displacement of the piston pin in the negative direction of the *y* axis. In some cases, we can see that pressure in the lubricating oil was lower than the saturated vapor pressure, so cavitation occurred. Cavitation occurred only at the boundary on the piston pin side; no vapor phase was observed at the boundary on the connecting rod side. This vapor phase distribution was likely caused by sudden displacement of the piston pin due to a large load. In addition, cavitation can be seen at the connecting end of the oil supply hole at the upper part of the analysis region. Because cavitation occurs at the connecting end on the side of rotational direction of the oscillating motion, it is likely because pressure sharply drops when lubricating oil is drawn into the gap from the oil supply hole. As the cavitation region increases, viscosity of the lubricating oil decreases, and the compression rate increases, so there is a possibility that sufficient film thickness cannot be realized when the combustion process completes and an upward force is applied to the piston pin (Figs. 27–29).

Figure 30 compares gap widths at the initial state and at *t* = 1.5 × 10^{−3} s, where load is maximized. In the rigid-body analysis, the minimum gap width narrowed to 1.8 *μ*m and 9.95 *μ*m as shown in Fig. 17, but in this analysis, the reduction was only 0.5 *μ*m, for a minimum gap width of 11.25 *μ*m. This is likely because elastic deformation of the structure suppresses the pressure increase, confirming that the gap width is overestimated by the rigid-body analysis alone, without considering elastic deformation.

### 3.8 Structure-Side Analysis Results.

Figures 31–33 show piston pin deformations and stress distributions. For ease of viewing, deformation amounts are magnified 100 times. The figures confirm that the piston pin exhibits a bow-like deformation as the load increases. The highest stress, approximately 500 MPa, is at the lower center of the piston pin. We can see that pressure in the lubricating oil was high at the lower center of the piston due to the squeeze effect, so the film thickness did not significantly change even when a large load was applied. Around the connecting rod edge, however, pushing the lubricating oil out suppresses pressure increase in the lubricating oil, and elastic deformation of the piston pin increases as the load increases. This is similar to the behavior predicted from the rigid-body analysis in Sec. 2.

### 3.9 Discussion of the Piston Pin.

Figure 34 shows deformation of the piston pin and the connecting rod at *t* = 1.5 × 10^{−3} s. Deformation in both parts is magnified 100 times for ease of viewing.

The film thickness is smallest at the lower end of the connecting rod when the piston pin exhibits a bow-like deformation. This part has the highest probability of seizure, and we can see that this is the same kind of seizure as the actual one shown in Fig. 1. We can thus say that this analysis qualitatively predicted the seizure location by calculation.

The amount of elastic deformation is much smaller in the connecting rod than in the piston pin, so deformation of the piston pin has the most significant effect on reducing the film thickness. In other words, minimizing bow-like deformation of the piston pin will reduce mechanical contact at the connecting rod end. An effective way of achieving this is to reduce the inner diameter of the piston pin. Piston pins generally have a hollow, cylindrical shape to reduce their weight, but while a larger inner radius makes the piston pin lighter, its rigidity against load is lessened. The analysis method developed in this research is expected to be very useful for numerical predictions of the minimal piston pin inner diameter that prevents seizure, allowing numerical analysis design using supercomputing without requiring high costs and long times for durability tests.

## 4 Conclusions

We conducted multiphase numerical analyses and investigations of thin-film fluid lubrication flow between a piston pin and connecting rod in an internal combustion engine, where there are particularly severe cavitating lubrication conditions at sliding parts with unsteady flow channel variation. We performed two analyses: a composite analysis with a rigid-body piston pin and a connecting rod, and a fluid–structure coupled analysis with a structure comprising both parts that can be elastically deformed. From these analyses, we obtained the following findings:

The analysis results for pressure distribution obtained from the composite rigid-body analysis were similar to those under the Sommerfeld condition, confirming the validity of this numerical analysis method and the analysis results. In the first half of the combustion process, when the load applied to the piston pin increases, pressure in the lubricating oil locally increases due to a squeeze effect resulting from rapid movement of the piston pin. At this time, movement of the piston pin caused large negative pressures on the side opposite of the high pressure, and cavitation was observed in some areas. Pressure increase in the lubricating oil was limited to the model’s central region, with no pressure increase observed around the edge.

We obtained the following findings from FSI analysis of elastic deformation of the piston pin and connecting rod. As in the rigid-body analysis, pressure in the lubricating oil increased in the loading direction due to the load applied to the piston pin, but the distribution was smoother than that in the rigid-body analysis. This is likely due to the effect of elastic deformation suppressing the pressure increase.

Fluid–structure coupled analysis showed bow-like deformation of the piston pin when a load is applied. We found that deformation of the connecting rod was smaller than that of the piston pin, and that the amount of piston pin deformation had a significant effect on the gap width between the two parts. We also showed that there is a high possibility of mechanical contact occurring at the connecting rod edge due to the bow-like deformation of the piston pin.

By comparison with a piston pin after actual use, we found that the fluid–structure coupled analysis qualitatively predicted the seizure location. One way to prevent piston pin seizure is to reduce its inner diameter to increase its rigidity. However, reducing the inner diameter leads to increased piston pin weight, so it is essential to use this numerical analysis method to select the optimal inner diameter. We showed that this numerical analysis method could predict seizure phenomena, which have been considered difficult to predict. This will allow some costly and time-consuming durability tests to be replaced with numerical analyses.

## Acknowledgment

The authors would like to thank Dr. Zheng Zhang and Mr. Kiyoshi Yoshinaga (Softflow Co., Ltd., Japan) for their helpful discussions. This work was partly supported by the GCORE Program and the Collaborative Research Project of the IFS, Tohoku University. Part of the numerical results in this research was obtained using supercomputing resources at the Cyberscience Center, Tohoku University. The numerical simulations were partially performed using FUJITSU AFI-NITY scalar parallel computers at the Institute of Fluid Science, Tohoku University.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*h*=specific enthalpy, J/kg

*k*=thermal conductivity, W/(m · K)

- $n$ =
unit vector normal to the surface

*p*=pressure, Pa

- $q$ =
heat flux vector, W/m

^{2}*r*=distance from the rigid-body, m

*D*=displacement, m

- $D$ =
displacement vector, m

- E =
Young’s modulus, Pa

**I**=unit tensor

*T*=temperature, K

- $U$ =
flow velocity vector, m/s

- $fs$ =
body force, N/m

^{3}- N
_{asp}= number of asperities per unit area

*D*/*Dt*=substantial derivative, 1/s

## Greek Symbols

*α*_{pv}=pressure–viscosity coefficient

*α*_{s}=volume expansion rate 1/K

*β*=temperature–viscosity coefficient

*β*_{asp}=radius of asperity curvature

- $\u03f5$ =
strain tensor

- $\gamma \u02d9$ =
shear rate, Pa

*η*=viscosity coefficient, Pa · s

*γ*=displacement diffusion coefficient

*μ*_{s},*λ*_{s}=Lame’s constants

*μ*=viscosity, Pa · s

- $\nu (poi)$ =
Poisson’s ratio

*ρ*=density, kg/m

^{3}- $\sigma Le$ =
linear elastic stress, Pa

- $\sigma $ =
stress tensor, Pa

- $\sigma T$ =
thermal stress, Pa

*σ*=normal stress, Pa

- $\tau $ =
viscous stress tensor, Pa

*τ*=shear stress, Pa

*τ*_{0}=Eyring stress, Pa