Coupled simulation of elastohydrodynamics and multi-flexible body dynamics in piston-lubrication system

Abstract

In this work, we propose a robust modeling and analysis technique of the piston-lubrication system considering fluid–structure interaction. The proposed schemes are based on combining the elastohydrodynamic analysis and multi-flexible body dynamics. In particular, multi-flexible body dynamics analysis can offer highly precise numerical results regarding nonlinear deformation of the piston skirt and cylinder bore, which can lead to more accurate results of film thickness for gaps filled with lubricant and of relative velocity of facing surfaces between the piston skirt and the cylinder block. These dynamic analysis results are also used in the elastohydrodynamic analysis to compute the oil film pressure and asperity contact pressure that are used as external forces to evaluate the dynamic motions of the flexible bodies. A series of processes are repeated to accurately predict the lubrication characteristics such as the clearance and oil film pressure. In addition, the Craig–Bampton modal reduction, which is a standard type of component mode synthesis, is employed to accelerate the computational speed. The performance of the proposed modeling schemes implemented in the RecurDyn™ multi-flexible body dynamics environment is demonstrated using a well-established numerical example, and the proposed simulation methods are also verified with the experimental results in a motor cycle engine (gasoline) which has a four cycle, single cylinder, overhead camshaft (OHC), air cooled.

Introduction

The reduction of automobile fuel consumption and exhaust gas emission is still important in the field of engine technology due to various environmental regulations. Those technologies are directly related to the reduction of power loss in the internal combustion engine.^1,2 The automobile lubrication system is therefore vital because most power loss comes from friction between the piston skirt and the inner surface of the cylinder block, ³ which means that precise predicting of piston movements and appropriate selection of design parameters are prerequisites to design and analysis of the lubrication system.^4,5

The main motion of the piston is the reciprocating translation, but the piston also has tilting motion that is mechanically generated by the rotation of the crank shaft. Although the amount of tilting is relatively smaller than the translation, it can lead to a large moment acting on the cylinder block. Therefore, analysis of the tilting has become a key issue in attempts to improve engine performance and reliability. ⁶

Lubrication analysis therefore allows us to precisely predict the behaviors of the piston skirt and compute the friction loss. In particular, elastohydrodynamic (EHD) lubrication has become popular because it can consider the effect of elastic deformation in computations of the lubrication film thickness. Since the 19th century, many approaches have been taken to tackle EHD lubrication. Zhu and Wang ⁵ summarized the history of EHD lubrication well. In 1886, based on experimental results, Reynolds ⁷ published a milestone for theoretical lubrication analysis in a journal bearing. Early studies have been focused on analytical and experimental methods or have focused on simplified numerical approaches due to the lack of computing power.^8–11 However, practical EHD lubrication has relied heavily on numerical approaches through model-based simulations. Therefore, because of the advancements of the computer technology, numerical EHD lubrication has been a main issue since the 1980s.^12–16 Although EHD lubrication analysis has been dramatically improved over the last 60 years, it has remained an active and attractive research issue. In particular, recent research in EHD lubrication has exhibited a multi-physics point of view, including attempts to solve the fluid–structure interaction (FSI) problem. Since multi-physics based on EHD lubrication directly incurs a heavy computational burden, its numerical efficiency is still a major bottleneck in practical engineering.

In this article, we propose an effective modeling procedure for the piston-lubrication system, a representative FSI problem, coupled with EHD and multi-flexible body dynamics (MFBD) simulation. The piston-lubrication analysis, in which elastic deformations of the piston and the cylinder block as well as overall nonlinear dynamic motion are taken into account, is performed in order to accurately understand the interaction characteristics. Li et al. ¹⁷ showed that lubrication behaviors between the piston rings and cylinder liners can indeed be affected by the dynamic deformations of cylinder structures through their research. Piston-lubrication analysis is based on EHD theory^4,5 which is used to represent the thickness and pressure of the fluid film generated by the relative motions of the objects, which exhibit elastic deformation, and the hydrodynamic pressure, which is computed using the Reynolds equation. ⁷ Greenwood and Tripp’s ¹⁸ asperity contact model is additionally used to cover the direct metal-to-metal contact in dry condition considering the roughness effect. In particular, unlike the typical rigid body approaches, MFBD simulation is used to precisely compute the nonlinear dynamic motions and the local elastic deformations of the piston skirt and the cylinder block. The flexible bodies of the piston and cylinder contain too much degree of freedom. Since the transient coupling analysis of MFBD and EHD requires excessive computation for structural analysis, so the modal reduction technique for the flexible bodies is applied for the appropriate computation and this method is deemed appropriate because the deformation of piston and cylinder due to oil pressure is not large. In a similar study, Kumar et al. ¹⁹ analyzed transient EHD using the modal reduction technique for the engine bearing and obtained stable results according to the appropriate mode selection. JH Schmidt ²⁰ suggested efficient numerical solution procedure for the lubricated impact problem which captures structural dynamic effects of complex geometries using modal reduction technique. Component mode synthesis (CMS), which is a representative modal reduction technique, is used to accelerate the simulation speed in this research.^21–28 We then employ the Craig–Bampton (CB) formulation, which is a standard CMS method. ²¹ The proposed modeling scheme of the piston-lubrication system is fundamentally an FSI problem; as such, it can be effectively tackled by an iterative co-simulation of the EHD and MFBD models, which together make a weakly coupled approach.

In the sections, “Summary of MFBD with modal reduction” and “CB-like formulation,” we first summarize the MFBD formulation with the modal reduction method. The EHD model is introduced in section “EHD lubrication.” The overall process for the numerical method proposed in this article is described in section “Numerical process.” Numerical testing using the proposed coupled modeling process of a simple piston-lubrication system is presented in section “Numerical studies.” In section “Experimental validations,” the proposed simulation scheme is also validated with the experimental results. The conclusions are given in section “Conclusion.”

Summary of MFBD with modal reduction

In this section, we summarize the MFBD formulation using the well-known floating reference frame.^29,30 We also present the coordinate reduction of the flexible body. The equations of motion of the multi-flexible body model on body i are

$${\mathbf{M}}^i{\stackrel{··}{\mathbf{q}}}^i+{\mathbf{K}}^i{\mathbf{q}}^i+{\mathbf{C}}_{{\mathbf{q}}^i}\mathbf{\lambda }={\mathbf{Q}}^i$$

(1)

where M is the mass matrix, K is the stiffness matrix, C is the kinematic constraint equation, q is the generalized coordinate vector, $\mathbf{\lambda }$ is the Lagrange multiplier vector, and Q is the force vector.

The generalized coordinate vector ${\mathbf{q}}^i$ can be defined as follows

$$\begin{array}{c}{\mathbf{q}}^i\approx {\mathbf{q}}_r^i+{\overline{\mathbf{q}}}^i={\mathbf{q}}_r^i+{\overline{\mathbf{T}}}^i{\mathbf{d}}^i,j=1,2,\dots ,N\hfill \\ {\overline{\mathbf{T}}}^i=\left[\begin{array}{cccc}{\mathbf{\phi }}_1\hfill & {\mathbf{\phi }}_2\hfill & \cdots \hfill & {\mathbf{\phi }}_j\hfill \end{array}\right],{\left({\mathbf{d}}^i\right)}^T=\left[\begin{array}{cccc}{a}_1\hfill & a_2\hfill & \cdots \hfill & a_j\hfill \end{array}\right]\hfill \end{array}$$

(2)

where ${\mathbf{q}}_r^i$ represents rigid body movement and ${\overline{\mathbf{q}}}^i$ represents the deformation of flexible body. Flexible body deformation can be approximated using a modal reduction technique. ${\overline{\mathbf{T}}}^i$ is the transformation (or projection) matrix that contains the mathematical basis vector ${\mathbf{\phi }}_j$ , and ${\mathbf{d}}^i$ is the reduced coordinate vector that contains the amplitude $a_j$ of ${\mathbf{\phi }}_j$ . The over bar means the approximated quantity. $N$ is the number of the mathematical basis, which can be used to approximate ${\overline{\mathbf{q}}}^i$ .

Using the transformation matrix ${\overline{\mathbf{T}}}^i$ , the components related to the flexible body in equation (1) can be projected into the reduced coordinates. Then, we can derive the reduced equations of motion from equation (1) as

$$\begin{array}{c}\left[\begin{array}{cc}{\mathbf{M}}_{\mathrm{rr}}^i\hfill & {\overline{\mathbf{M}}}_{\mathrm{rf}}^i\hfill \\ \multicolumn{1}{c}{{\overline{\mathbf{M}}}_{\mathrm{rf}}^{\mathrm{iT}}}& {\overline{\mathbf{M}}}_{\mathrm{ff}}^i\hfill \end{array}\right]\left(\begin{array}{c}{\stackrel{··}{\mathbf{q}}}_r^i\hfill \\ \multicolumn{1}{c}{{\stackrel{··}{\mathbf{d}}}^i}\end{array}\right)+\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\overline{\mathbf{K}}}_{\mathrm{ff}}^i\hfill \end{array}\right]\left(\begin{array}{c}{\mathbf{q}}_r^i\hfill \\ \multicolumn{1}{c}{{\mathbf{d}}^i}\end{array}\right)\hfill \\ +\left[\begin{array}{c}{\mathbf{C}}_{q_r^i}^T\hfill \\ \multicolumn{1}{c}{{\left[{\overline{\mathbf{T}}}^i\right]}^T{\mathbf{C}}_{q_r^i}^T+{\mathbf{C}}_{d^i}^T}\end{array}\right]\mathbf{\lambda }=\left(\begin{array}{c}{({\mathbf{Q}}^i)}_r\hfill \\ \multicolumn{1}{c}{{({\overline{\mathbf{Q}}}^i)}_f}\end{array}\right)\hfill \end{array}$$

(3)

and the component matrices and vector are computed as

$$\begin{array}{c}{\overline{\mathbf{M}}}_{\mathrm{ff}}^i={\left[{\overline{\mathbf{T}}}^i\right]}^T{\mathbf{M}}^i{\overline{\mathbf{T}}}^i=\mathbf{I},{\overline{\mathbf{K}}}_{\mathrm{ff}}^i={\left[{\overline{\mathbf{T}}}^i\right]}^T{\mathbf{K}}^i{\overline{\mathbf{T}}}^i,\hfill \\ {\left({\overline{\mathbf{Q}}}^i\right)}_f={\left[{\overline{\mathbf{T}}}^i\right]}^T{\left({\mathbf{Q}}^i\right)}_r,{\overline{\mathbf{M}}}_{\mathrm{rf}}^i={\left[{\overline{\mathbf{M}}}_{\mathrm{fr}}^i\right]}^T\hfill \end{array}$$

(4)

In addition, ${\mathbf{M}}_{\mathrm{rr}}^i$ is calculated from the total mass of flexible body i. Coupling term ${\overline{\mathbf{M}}}_{\mathrm{rf}}^i$ is a function of terms such as transformation matrix ${\overline{\mathbf{T}}}^i$ and mass matrix $\mathbf{M}$ . Yoo and Haug derived and verified these detailed equations in their research.^31,32

Various projection techniques for coordinate reduction have been proposed since the 1960s; they can be categorized according to which mathematical basis is employed. Among the various techniques, eigen solutions, constraint modes, the Krylov subspace scheme, and singular value decomposition (SVD) are the mainly used basis vectors.^33–41 Since these basis vectors have their own advantages and disadvantages, users have to thoughtfully select an optimal basis vector, or have to combine them to accurately approximate the flexible model.

Finally, the global position vector ${\mathbf{r}}_p^i$ of an arbitrary node p in the i body is computed using the numerical results in equation (3) as shown in Figure 1. It is given as follows

$${\mathbf{r}}_p^i={\mathbf{R}}^i+{\mathbf{A}}^i\left({\mathbf{\Omega }}^i\right)({\left({\mathbf{q}}_0^i\right)}_p^{\prime }+{\left({\overline{\mathbf{T}}}^i{\mathbf{d}}^i\right)}_p)$$

(5)

where ${\mathbf{R}}^i$ is the translation vector from ground inertia to reference frame and ${\mathbf{A}}^i$ is the transformation matrix related to the body reference frame, which determined by the rotation vector ${\mathbf{\Omega }}^i$ that represents rotation of rigid body motion. Both variables determine the coordinate ${\mathbf{q}}_r^i$ associated with the rigid body motion of equation (3). $({\mathbf{q}}_0^i){\prime }_p$ indicates the constant displacement vector from body reference frame to the p node at undeformed state and the deformation of p node is then expressed as the product of the transformation matrix ${\overline{\mathbf{T}}}^i$ and the reduced coordinate vector ${\mathbf{d}}^i$ .

Figure 1.

Nodal position of the flexible body.

The combustion engine consists of four main parts: connecting rod, crank shaft, piston, and cylinder blocks. In particular, the piston and cylinder block are the main components in the typical piston-lubrication system; they should be modeled as flexible bodies for EHD lubrication: Figure 2 shows a schematic diagram of the piston and the cylinder block in the piston-lubrication system. The each reference frame of the piston and the cylinder block are denoted using the subscripts p and c at a certain time t.

Figure 2.

Schematic of the piston and cylinder block structure with reference frame.

The arbitrary position vector of flexible piston and cylinder $\mathbf{r}$ and its velocity vector $\stackrel{·}{\mathbf{r}}$ are used in the EHD analysis. However, the external force $({\mathbf{Q}}^i{)}_r,({\overline{\mathbf{Q}}}^i{)}_f$ in equation (3) is determined by the results for the hydrodynamic pressure in the EHD analysis. This issue is presented in section “EHD lubrication.”

CB-like formulation

The CMS method, combining the eigen solutions and constraint modes, is popular and widely used for coordinate reduction of a linear elastic body in the multibody dynamics (MBD) field. In particular, the CB-like formulation, ²¹ which is a standard CMS method, is employed in this work. In the CB method, each body is connected by the joints and there are some parts which force is applied to. The degrees-of-freedom (DOFs) associated with these parts is set as the interface and are called the interface DOF. The other DOFs in the configuration of the flexible body are defined as the internal DOFs. With an assumption of free vibration, we then rewrite the component matrices and vectors in equation (1) as

$$\begin{array}{c}{\mathbf{M}}^i=\left[\begin{array}{cc}{\mathbf{M}}_{11}^i\hfill & {\mathbf{M}}_{12}^i\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}_{21}^i}& {\mathbf{M}}_{22}^i\hfill \end{array}\right],{\mathbf{K}}^i=\left[\begin{array}{cc}{\mathbf{K}}_{11}^i\hfill & {\mathbf{K}}_{12}^i\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{21}^i}& {\mathbf{K}}_{22}^i\hfill \end{array}\right]\hfill \\ {\mathbf{M}}_{12}^i={\left[{\mathbf{M}}_{21}^i\right]}^T,{\mathbf{K}}_{12}^i={\left[{\mathbf{K}}_{21}^i\right]}^T,{\mathbf{q}}^i=\left(\begin{array}{c}{\mathbf{q}}_1^i\hfill \\ \multicolumn{1}{c}{{\mathbf{q}}_2^i}\end{array}\right)\hfill \end{array}$$

(6)

Transformation matrix ${\overline{\mathbf{T}}}^i$ in the CB method is defined by combining the fixed interface normal modes ${\mathbf{\zeta }}^i$ and the constraint modes ${\mathbf{ѱ}}^i$ , as follows

$$\begin{array}{c}{\overline{\mathbf{T}}}^i=\left[\begin{array}{cc}{\mathbf{\Psi }}^i\hfill & {\mathrm{Z}}^i\hfill \end{array}\right]\hfill \\ {\mathbf{\Psi }}^i=\left[\begin{array}{cccc}{\mathbf{ѱ}}_1\hfill & {\mathbf{ѱ}}_2\hfill & \cdots \hfill & {\mathbf{ѱ}}_j\hfill \end{array}\right],j=1,2,\dots ,N_1\hfill \\ {\mathbf{Z}}^i=\left[\begin{array}{cccc}{\mathbf{\zeta }}_1\hfill & {\mathbf{\zeta }}_2\hfill & \cdots \hfill & {\mathbf{\zeta }}_k\hfill \end{array}\right],j=1,2,\dots ,N_d\hfill \end{array}$$

(7)

where $N_1$ and $N_d$ are the numbers of the constraint modes and the fixed interface normal modes, respectively. ${\mathbf{\Psi }}^i$ and ${\mathbf{Z}}^i$ can be computed as

$${\mathbf{\Psi }}^i=\left[\begin{array}{c}\mathbf{I}\hfill \\ \multicolumn{1}{c}{-{[{\mathbf{K}}_{22}^i]}^{-1}{\mathbf{K}}_{21}^i}\end{array}\right],{\mathbf{Z}}^i=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{{\mathbf{X}}_d^i}\end{array}\right]$$

(8)

in which $\mathbf{I}$ is an identity matrix and ${\mathbf{X}}_d^i$ is the dominant eigenvector matrix, calculated by the following eigenvalue problem

$$[{\mathbf{K}}_{22}^i-{\left({\omega }_k^i\right)}^2{\mathbf{M}}_{22}^i]{\mathbf{\chi }}_k^i=0,k=1,2,\dots ,N_2$$

(9)

where ${\omega }_k^i$ and ${\mathbf{\chi }}_k^i$ are the kth natural frequency and the mode shape of the internal DOFs, respectively. $N_2$ is the number of the internal DOFs. ${\mathbf{\chi }}_k^i$ then satisfies the mass orthonormal condition. The constraint mode matrix ${\mathbf{\Psi }}^i$ in equation (8) is derived in the static condition only, which is known as the static correction. Substituting equation (7) into equation (4), we can obtain the details of the reduced equations of motion in equation (3) using the CB approach. The frequency cut-off rule, employed by ranking the natural frequency ${\omega }_k^i$ , is a standard way to select the dominant $N_d$ modes among the entire $N_2$ modes that can be computed in equation (9) $(N_d=N_2)$ . However, the frequency cut-off rule should be carefully used in the MFBD simulation because ${\omega }_k^i$ is just a factor that can show the contribution of the kth mode in the flexible body i itself, which means that it has no physical meaning in the entire multibody system and the system natural frequency. ²⁸ Consequently, the mode shapes of ${\mathbf{\Psi }}^i$ and ${\mathbf{Z}}^i$ play important roles in the accuracy of the reduced flexible model. This problem may be used to tackle the high-fidelity mode selection criteria based on the modal contribution. This issue is well summarized in references.^42,43 In particular, for high precision problems, robust dynamic correction with residual flexibility may be required for better approximation of the reduced system.^26–28 In this research, the piston-lubrication system focuses on the EHD analysis of the piston skirt and cylinder inner surface. Therefore, the piston skirt and the cylinder inner surface are set as interface nodes, and the mode shapes in the radial direction perpendicular to the axial direction of the piston and cylinder are selected for more accurate analysis.

EHD lubrication

EHD lubrication is a mode of fluid-film lubrication in which hydrodynamic behavior is significantly increased by surface elastic deformation; the viscosity is also increased due to high pressure. The dynamic behavior of a lubricant film between a piston and a cylinder block is described by a Reynolds equation.^4,5 In addition, Greenwood and Tripp’s asperity model, which is a well-known contact theory that considers the roughness effect, is used to compute the contact pressure and force in dry condition. ¹⁸

Reynolds equation for flexible piston-lubrication system

In order to perform EHD lubrication, the lubrication area and the coordinate system for lubricant behavior should be defined first because the governing equation for the lubricant is then formulated based on this coordinate system. The coordinate system is located at the center position of the cylinder in its initial state and its y-axis, which is cross product of z- and x-axis in Figure 3, is set parallel to the axial direction of the cylinder. It is a fixed coordinate system for defining variables for the lubrication analysis. Figure 3 shows an arbitrary cross section of the piston-lubrication model. h is the lubrication film thickness which represents clearance between cylinder and piston in radial direction on the $x-z$ plane, $r_p$ is the piston radius, $v_{\theta }$ is the angular velocity of the piston, and $v_z$ is the translational velocity of the piston. Using a Newtonian fluid and neglecting the gravity force by fluid mass, the general form of the Reynolds equation is then written in the cylindrical coordinate system as book by Martin ¹⁰

$$\begin{array}{c}\frac{1}{r_p^2}\frac{\partial }{\partial \theta }\left(\frac{{\rho }_L{h}^3}{\mu }\frac{\partial p_H}{\partial \theta }\right)+\frac{\partial }{\partial z}\left(\frac{{\rho }_L{h}^3}{\mu }\frac{\partial p_H}{\partial z}\right)\hfill \\ =6v_{\theta }\frac{\partial \left({\rho }_{L}h\right)}{\partial \theta }+6v_z\frac{\partial \left({\rho }_{L}h\right)}{\partial z}+12\frac{\partial \left({\rho }_{L}h\right)}{\partial t}\hfill \end{array}$$

(10)

where $p_H$ is the hydrodynamic pressure, $\mu$ is the dynamic lubricant viscosity, ${\rho }_L$ is the density of the lubricant, and $t$ is the time variable. In the Barus ⁴⁴ pressure–viscosity relationship, $\mu$ is computed as

$$\mu ={\mu }_0{e}^{\alpha p-\gamma \mathrm{\Delta }T}$$

(11)

in which ${\mu }_0$ is the dynamic lubricant viscosity measure at atmospheric pressure, $\alpha$ is the viscosity–pressure coefficient, $\gamma$ is the thermal coefficient, and T denotes the temperature. The density is commonly defined as the Dowson–Higginson pressure–density function, ⁴⁵ as follows

$${\rho }_L={\rho }_{L_0}+{\rho }_{L_0}\frac{0.6\times {10}^{-9}{p}_H}{1+1.7\times {10}^{-9}{p}_H}+0.0007{\rho }_{L_0}\mathrm{\Delta }T$$

(12)

where ${\rho }_{L_0}$ is the fluid density measured at the normal temperature. Regarding an arbitrary cross section along the y-axis, equation (10) can be simplified in the Cartesian coordinate system as

$$\begin{array}{c}\frac{\partial }{\partial x}\left(\frac{{\rho }_L{h}^3}{\mu }\frac{\partial p_H}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{{\rho }_L{h}^3}{\mu }\frac{\partial p_H}{\partial y}\right)\hfill \\ =6v_x\frac{\partial \left({\rho }_{L}h\right)}{\partial x}+6v_y\frac{\partial \left({\rho }_{L}h\right)}{\partial y}+12\frac{\partial \left({\rho }_{L}h\right)}{\partial t}\hfill \end{array}$$

(13)

where x denotes the circumferential direction coordinate and y denotes the axial direction coordinate. Considering the finite volume discretization in equation (13), the film thickness and the velocity vectors can be defined by the nodal position vector $\mathbf{r}$ and its velocity $\stackrel{·}{\mathbf{r}}$ , which are the results of the MFBD simulation in equation (3), as follows

$$h=h\left(\mathbf{r}\right),v_x=v_x\left(\stackrel{·}{\mathbf{r}}\right),v_y=v_y\left(\stackrel{·}{\mathbf{r}}\right)$$

(14)

Figure 3.

Cross-section of the piston-lubrication model.

Then, the following numerical system of equations for the pressure $p_H$ can be evaluated using the film thickness and the velocity components.

Asperity contact model

When the film thickness is very thin, direct metal-to-metal contact can occur, which means that massive pressure might be added to the hydrodynamic pressure computed by the Reynolds equation. In particular, we here use the well-known Greenwood and Tripp ¹⁸ asperity contact model to consider the surface roughness in the piston skirt and inner surface of the cylinder block. The asperity contact pressure can be calculated as follows

$$\begin{array}{c}{p}_A\left(h_r\right)=k_A{E}_A{q}_A(h_r),h_r=\frac{h}{{\sigma }_s}\hfill \\ q_A\left(h_r\right)=\{\begin{array}{c}4.4086\times {10}^{-5}{(4-h_r)}^{6.804},h_r<4.0\hfill \\ \multicolumn{1}{c}{0,h_r\ge 4.0}\end{array}\hfill \end{array}$$

(15)

where $k_A$ is the generalized stiffness parameter, $E_A$ is the elastic modulus, and ${\sigma }_s$ is the surface roughness. In the discretized mesh grid, the asperity contact pressure $p_A$ can be defined, and then the total pressure p, which is a result of the EHD model, is obtained as

$$p=p_H+p_A$$

(16)

The calculated total pressure p is calculated as the nodal force by multiplying the characteristic area for each node of the flexible body. This nodal force acts in the centripetal direction of the coordinate system for lubrication analysis and applies to $({\mathbf{Q}}^i{)}_r$ in equation (3). In addition, this nodal force is converted to a modal force by multiplying the transformation matrix ${\overline{\mathbf{T}}}^i$ and then applied to $({\overline{\mathbf{Q}}}^i{)}_f$ .

Modal force by gas combustion pressure

In the piston-lubrication system, pressure resulting from gas combustion should be considered to check the characteristic of piston-lubrication system because it is a main factor to drive system. In order to apply gas force to modal reduced flexible body, first we have to compute nodal force resulting from the gas pressure. The nodal force ${\mathbf{f}}^i$ for arbitrary node i can be calculated as follows

$${\mathbf{f}}^i=\frac{\left(\sum _{j=1}^{n^i}{\left(p_{\mathrm{gas}}\right)}_j{s}_j{\hat{\mathbf{k}}}_j\right)}{n^i}$$

(17)

where $(p_{\mathrm{gas}}{)}_j$ is the pressure by gas combustion on element patch j to which node i belongs and $s_j$ is the area of it. ${\hat{\mathbf{k}}}_j$ is unit normal vector to the patch j and $n^i$ is the number of patch that shares node i. Finally, this nodal force ${\mathbf{f}}^i$ is also employed to define $({\mathbf{Q}}^i{)}_r,({\overline{\mathbf{Q}}}^i{)}_f$ for the MFBD simulation (see equation (3)).

Numerical process

In this work, we focus on developing a co-simulation procedure of EHD and MFBD solvers, which is highly effective from a computational point of view, and which is a kind of weakly coupled simulation. The numerical process is specifically described in Figure 4.

Figure 4.

Numerical process of piston-lubrication model.

First, transformation matrix ${\overline{\mathbf{T}}}^i$ of piston-lubrication system is computed using mass and stiffness matrix of them through CB method. And then, reduced MFBD solver calculates nodal position $\mathbf{r}$ and velocity $\stackrel{·}{\mathbf{r}}$ . From this information, EHD solver calculates the boundary conditions for such as $h,v_x,v_y$ and calculates the Reynolds equation for hydrodynamic pressure $p_H$ and asperity contact pressure $p_A$ . External force vector $({\mathbf{Q}}^i{)}_r,({\overline{\mathbf{Q}}}^i{)}_f$ can be evaluated from this pressure information and it is applied to reduced MFBD solver again. Above steps are repeated until residual equation of implicit integrator converges. And if it converges, time step is increased and this procedure is repeated.

Numerical studies

The proposed coupled modeling process of the piston-lubrication system is implemented in the RecurDyn™ MFBD environment. ⁴⁶ We here consider a piston-lubrication example to evaluate the performance of the proposed numerical process (see Figure 5). In the example model, the rod and crank shaft are connected by revolute joint and they are modeled as rigid bodies. Driving force of piston-lubrication system is derived from rotation of crank shaft. And, the rod and piston are connected by revolute joint in point A at which master node are connected to the slave nodes on flexible piston. The piston and cylinder block, which are related to the lubrication region, are modeled as flexible bodies. Parabolic tetrahedral solid elements are used for the piston and the cylinder block. The axial length of inner space of cylinder block is 95 mm and 40 by 40 fluid mesh grids are distributed in the lubrication area between the piston skirt and the inner surface of the cylinder block. In order to calculate the Reynolds equation for the lubrication area, first define the coordinate system for the Reynolds equation as shown in Figure 5. Here, the circumferential direction is defined as the direction of rotation from the z-axis to the x-axis, the axial direction is the y-axis, and the radial direction for defining the film thickness is defined on the $x-z$ plane. If we know the position and velocity of the nodes on the piston skirt and the cylinder inner surface, we can interpolate to find the relative velocity in the circumferential direction, the relative velocity in the axial direction, and the film thickness. Based on this, lubrication analysis is performed to find the pressure distribution and then convert it into a force vector, and apply the force vector to the nodes on piston skirt and cylinder inner surface (Figure 6). The numerical parameters for EHD of the example model are presented in Table 1. In order to consider the effect of deformation on the lubrication film, the mode shapes are selected so that deformation in the direction that affects the film thickness is mainly set.

Figure 5.

Numerical example of piston-lubrication system.

Figure 6.

Examples of selected mode shapes.

Table 1.

Asperity contact parameters.

Parameter	Notation	Unit	Value
Roughness	${\sigma }_s$	mm	6.3617e–04
Composite elastic modulus	$E_A$	MPa	94,747.6
Generalized stiffness parameter	$k_A$	–	6.714e–02
Dynamic viscosity	$\mu$	Pa s	2.28e–03

For the numerical study, we consider a crank rotation as $V_{\mathrm{rpm}}=7000$ , which generates the driving force. For numerical case 1, 15 piston modes and 15-cylinder block modes are used to coordinate the reduction for the two flexible bodies; 60 piston modes and 60-cylinder block modes are used in numerical case 2. We then compare the numerical results of the MBD and MFBD approaches. In the MBD simulation, every component in piston-lubrication system is regarded as a rigid body.

For the two numerical cases, we plot the lateral motion and the rocking angle of point A (see Figure 5), which is the center point of the piston connected with the rigid body elements (RBEs). Here, lateral motion indicates movement of point A in x direction and rocking angle which is the piston’s secondary motion indicates the rotation of point A in z direction which is cross product of x- and y-axis. They can be seen in Figures 7 and 8. Numerical results show that the lateral motions and rocking angle of the piston, considering the MFBD simulation, might be different from the results obtained using MBD, although all results show the same tendency. This implies that the elastic behaviors of the piston and the cylinder block can lead to significant effects in the EHD lubrication model. It also implies that a great enough number of modes that can represent the elastic behaviors of the flexible body well might be useful in the modal reduction model.

Figure 7.

Lateral motion of point A.

Figure 8.

Rocking angle of point A.

We also investigate the pressure and thickness distribution according to the circumferential angle. Those results clearly show different phenomenon in the MBD and MFBD approaches. We found that high pressure occurs when the simulation time is 0.018 s, as can be seen in Figure 9. It is also seen that the contours of the MBD and MFBD approaches can be clearly distinguished from each other. As shown in Figure 10, in the MBD case, except for the elastic modal effect, the maximum pressure occurs at the center of the piston skirt; however, in the MFBD case, the maximum pressure occurs at the circumferential end of the piston skirt. This phenomenon can be understood through the thickness curve plotted in Figure 11. It can be seen that in the case of the MFBD results, the piston and the cylinder block can represent elastic deformation in which the piston skirt becomes bent inward. Thus, the oil film thickness increases around the center of the piston skirt. This causes the hydrodynamic pressure to decrease around it. Meanwhile, there is no increase of film thickness in the case in which only the rigid body motion is considered; then, the distribution of the hydrodynamic pressure is different from that in the MFBD results.

Figure 9.

Pressure contours of the MBD and MFBD approaches in case 1: MBD (left) and MFBD (right).

Figure 10.

Pressures according to angle on the left-hand side of the piston skirt.

Figure 11.

Film thicknesses according to angle on the left-hand side of the piston skirt.

Figures 12 and 13 show minimum thickness and friction force in axial direction for the rigid case and a case with 60 modes during the operation time. In the case of friction force, the shear stress in the Newtonian fluid due to the viscosity and the friction force due to the asperity contact are combined. In the upstroke situation, the friction force is smaller when the flexible body (a case with 60 modes) is shown, which may be related to the larger minimum thickness of the flexible body under similar circumstances.

Figure 12.

Minimum thickness.

Figure 13.

Friction force.

We here examine the computational costs of increasing the number of modes in the MFBD analysis. The results are summarized in Table 2. To theoretically compute the operation counts, we here consider two major computational processes: number of Jacobian evaluations and number of residuals. It can be seen that the computational cost dramatically increases as the number of modes increases. The reason for this is that the operation of Jacobian evaluations, which are derivatives of the pressure force by the system coordinates to meet the force equilibrium in an implicit integrator, increases dramatically. It should also be noted that the present piston-lubrication model is quite simple; it is modeled for the feasibility study only. This means that the practical model might be much more complicated and have around a million DOFs; this allows us to conclude that if the interpretation is made through a full nodal analysis rather than through the modal model reduction approach, its computational burden will not be acceptable in engineering practice.

Table 2.

Comparison of computational costs: number of Jacobian evaluations (NJAC) and number of residual evaluations (NRES).

Number of modes		Operation counts		Time (s)
Piston	Cylinder block	NJAC	NRES
15	15	108,565	322,108	2229
30	30	279,576	819,395	13,529
60	60	349,796	1,022,366	26,991

Experimental validations

We here consider another piston-lubrication example (see Figure 14) to check how valid the numerical results are. This example model is compared with the experimental results by Honda Motors for motor cycle engine (gasoline) which has four cycles, single cylinder, overhead camshaft (OHC), air cooled to verify the numerical results. The measurement results such as lateral motion and rocking angle are compared to the numerical analysis results. In the example model, the rod and crank shaft are connected by revolute joint and they are modeled as rigid bodies. Driving force of piston-lubrication system is derived from rotation of crank shaft. The rod and piston are then connected by revolute joint in point A at which master nodes (see Figure 14) are connected to the slave nodes on flexible piston. The piston and cylinder block, which are related to the lubrication region which is formed between piston skirt and inner surface of cylinder block, are modeled as flexible bodies. Parabolic tetrahedral solid elements are used for the piston and the cylinder block. The axial length of inner space of cylinder block is 103.55 mm and 40 by 40 fluid mesh grids are distributed in the lubrication area between the piston skirt and the inner surface of the cylinder block. As in the previous test case, we define a coordinate system for calculating the Reynolds equation and obtained the boundary conditions accordingly. In order to obtain this, the position and velocity information of the nodes of the piston skirt and the cylinder inner surface nodes were interpolated. The numerical parameters for EHD of the example model are presented in Table 3. As in the test case of section “Numerical studies,” in order to consider the effect of deformation on the lubrication film, the mode shapes are selected so that deformation in the direction that affects the film thickness is mainly set (Figure 15).

Figure 14.

Validation model of piston-lubrication system.

Table 3.

Asperity contact parameters.

Parameter	Notation	Unit	Value
Roughness	${\sigma }_s$	mm	6.3617e–04
Composite elastic modulus	$E_A$	MPa	94,747.6
Generalized stiffness parameter	$k_A$	–	6.714e–02
Dynamic viscosity	$\mu$	Pa s	2.28e-03

Figure 15.

Selected mode shape for validation model.

For the numerical study, we consider a crank rotation as $V_{\mathrm{rpm}}=7000$ , which generates the driving force. In total, 383 piston modes and 1022-cylinder block modes are used. The reason for using many modes is to realize accurate deformation.

In order to reflect the effect of combustion pressure on the piston top, measurement combustion pressure according to crank rotation is applied to the piston top, see Figure 16. Combustion pressure is maximized when crank angle reached to 10° and this type of combustion pressure is repeatedly applied to the piston during the entire simulation. This combustion pressure is employed to define nodal force for reduced flexible body (see equation (17)), and the resulting nodal force is also employed to define $({\mathbf{Q}}^i{)}_r,({\overline{\mathbf{Q}}}^i{)}_f$ for the MFBD simulation (see equation (3)) with the oil film and the asperity contact pressure.

Figure 16.

Application of combustion pressure to piston top.

For the validation test with experimental result, we plot the lateral motion and the rocking angle of point A (see Figure 14), which is the center point of the piston connected to the RBEs. Here, the definition of lateral motion and rocking angle is the same as in the test case in section “Numerical studies.” They can be seen in Figures 17 and 18. Numerical results expressed according to crank rotation during two cycles show similar tendency to the experimental results. In particular, it shows similarities in overall amplitude variation.

Figure 17.

Lateral motion at point A according to rotation of crank shaft for the validation test model (third and fourth stroke cycle).

Figure 18.

Rocking angle at point A according to rotation of crank shaft for the validation test model (third and fourth stroke cycle).

Conclusion

In this study, we proposed an efficient coupling method for the EHD and MFBD analyses for a piston-lubrication system. In order to convince ourselves of the computational efficiency, the modal reduction model with standard CB formulation is employed. In order to analyze the lubrication characteristics of the piston skirt and the cylinder inner surface, the lubrication area is set to the interface node, and the mode shape that influences the lubrication analysis is selected. In the lubrication analysis, to cover metal-to-metal contact in dry condition, we also consider the asperity contact model. In the proposed modeling procedures, EHD and MFBD simulations are repeated by updating the elastic deformation and the hydrodynamic pressures. By developing the efficient coupling method between MFBD and EHD, the users can analyze the system which includes piston-lubrication problem in multibody system. The numerical results show that the elastic behaviors of the flexible bodies can play significant role in the computation of the film thickness, which implies that EHD lubrication with MFBD simulation might be considered to precisely predict piston-lubrication behaviors. The numerical results also demonstrate that coordinate reduction in the MFBD simulation is necessary to effectively handle this problem in engineering practice, and the proposed numerical method is verified to some extent through comparing to experimental results.