© R. Hamel et al., Published by EDP Sciences 2023

1 Introduction

Automotive industry requirements are currently of increasing importance. Each new generation of thermal engines requires a formulation of more efficient lubricants. The current trend in the design of modern engines is to optimize their operating characteristics in terms of power, fuel consumption, and compliance with increasingly stringent environmental standards. A consequence of these requirements is the severity of the stresses to which the various engine components are subjected, which results in an increase of undesirable friction. It is well known that friction is one of the main causes of loss of efficiency in an internal combustion engine. Indeed, from 10 to 20% of the power contained in the engine fuel is absorbed by the mechanical power which is dissipated by friction [1]. In this regard, petroleum refiners and engine lubricant manufacturers must be able to formulate lubricating oils that meet the increasingly stringent specifications such as environmental standards, but also efficient in terms of reducing friction and wear of moving parts.

To reduce friction in an engine, two technological solutions can be considered. The first solution consists of treating the active contact surfaces by depositing coatings such as thin layers of amorphous carbon (DLC: Diamond-Like Carbon), lead or tin-based white metals (Babbitts), and polyetheretherketone (PEEK composites). The properties of these materials to reduce friction and wear have turned out to be very interesting [2]. The major drawback of such a solution lies in the difficulty, or even the impossibility, of depositing such coatings in certain geometrically complex areas of the engine.

The second solution concerns the choice of a more efficient lubricant with physicochemical properties other than its viscosity. It is then a question of formulating lubricants by adding additives to the base lubricant which can act when the lubricating film is no longer thick enough to avoid the destructive risks of direct metal-to-metal contact between the rough surfaces. To this end, the inclusion of anti-wear additives to the base lubricant constitutes a technological solution to enable a reduction of friction and wear rates, especially in the mixed and boundary lubrication regimes.

It should be noted that the additives with tribological action such as anti-wear additives currently used in the formulation of traditional motor oils are organic compounds based on sulfur and phosphorus. These two compounds are the cause of the deterioration of vehicle gas treatment systems and are called into question by anti-pollution standards. The substitution of these molecules by lubricating nanoparticles will allow the formulation of lubricating nanofluids (i.e., nanolubricants) that are more efficient and more respectful of the increasingly strict environmental standards.

The harmfulness of sulfur and phosphorus emissions of which the two molecules MoDTC and ZnDTP are at the origin, as well as the ineffectiveness of these two molecules at cold temperatures, have recently mobilized significant research on the development of new alternatives to these molecules. Such new materials should respect the environmental standards and have tribological properties at least equivalent to those of the additives being used currently. This research has led to a completely revolutionary new nanotechnology approach in the field of lubricating additives through the development of a new family of nanoparticle (NP) additives. It has been shown that some of these nanomaterials, when used in dispersion as lubricating additives, have, sometimes, clearly superior lubrication performance compared to those of traditional additives [3–5]. In addition, these nanoparticles are found to be more chemically inert than traditional molecular additives and therefore potentially less harmful to the environment. Accordingly, the use of nanotechnologies in the context of tribological applications (lubrication) is currently seeing increasing interest. This interest has resulted in a substantial development of nanomaterial synthesis techniques allowing the development of nano-objects of different structures (fullerenes, nanotubes, etc.) and of different chemical nature (metallic, oxides, carbonaceous, etc.). Numerous studies on the tribological properties (friction and wear) of these nano-objects dispersed in a lubricating base have shown substantial reductions of wear as nanoparticles were added to oil bases. According to Chen et al. [6], the multi-layered carbon nanotubes dispersed in oil reduce friction by 10% and wear by 30% compared to the lubricating base.

On the other hand, the inorganic fullerenes IF-MoS₂ and IF-WS₂, synthetized for the first time in 1992 by Tenne et al. [7], can reduce friction and wear by up to 70% compared to the lubricating base [8,9].

From the tribological point of view, fullerenes are then the most interesting systems and they can be used as potential substitutes to currently used molecular lubricating additives since they are environmental pollutants and less efficient in service under certain specific conditions. Indeed, the MoS₂ and the WS₂ act like a third body in lubricated rough contacts and essentially behave like tiny microscopic ball bearings thus reducing the undesirable friction between the two metal surfaces in contact as illustrated in Figure 1.

Many key questions were addressed such as the involved lubrication mechanisms as well as the influence of the concentration, size, structure, morphology, and aggregation of nanoparticles on the nanolubricant rheological behavior [10]. Experimental rheological studies demonstrated that nanolubricants (base oil + nano-additives) have a non-Newtonian rheological behavior and their flow cannot be described by the classical continuum mechanics theory that neglects the size of particles [11]. This requires the use of alternative continuum theories that accurately describe the hydrodynamic flow of nanolubricants in the bearing clearance space, considering the effects of both concentration and size of NPs. Among them, the Vijay Kumar Stokes micro-continuum theory [12–14] was successfully applied to this kind of fluids and the theoretical analyses revealed an increase in load carrying capacity and an enhancement of dynamic stability of journal bearings [15–20].

As far as the authors know, very few research investigations have been carried out to study the journal bearing dynamic behavior operating with NP based lubricants. Nair et al. [21] showed, for a given eccentricity ratio, that the addition of NPs on commercially available lubricant may enhance the viscosity of lubricant and, in turn, change the steady-state and dynamic performance characteristics such as the load capacity, attitude angle, end leakage, the stiffness and damping coefficients, and the whirl frequency.

Nair et al. [22] studied the static performance characteristics in terms of load capacity, friction force, end leakage, and attitude angle using iso-viscous and thermo-viscous approaches for different values of eccentricity ratios and lubricants containing different nanoparticles (CuO, CeO₂, and Al₂O₃). They found that for the thermoviscous case at high values of eccentricity, addition of nanoparticles produces a significant effect on the performance characteristics. The load capacity and friction force of the bearing increase with the increase of concentration of nanoparticles for both isoviscous and thermoviscous cases. They also found that the end leakage and attitude angle decrease with increase of concentration of nanoparticles in both isoviscous lubricants and lubricants that contain nanoparticles.

Shenoy et al. [23] analyzed the effects of CuO, TiO₂ and Nano-diamond nanoparticles additives in SAE 30 engine oil on static characteristics of an externally adjustable fluid-film bearing. The obtained results showed, in the case of a bearing having negative radial and negative tilt adjustments, and operating with SAE engine oil added with TiO₂ nanoparticles showed that the carrying load capacity was higher by about 23% and 35% than that obtained for SAE engine oil without nano-additives and base oil, respectively. Furthermore, the friction force of an externally adjustable fluid-film bearing increased by about 25% for the SAE engine oil without nanoparticles compared to the same oil blended with TiO₂ nanoparticle additives.

Kalakada et al. [24] investigated the static performance characteristics of thermohydrodynamic journal bearings operating under nanolubricants (lubricants containing per cent weight concentration of nanoparticles). The lubricant and the nanoparticles used for the present work are SAE 15W40, copper oxide (CuO), cerium oxide (CeO₂), and aluminum oxide (Al₂O₃). Their study was focused on the distributions of pressure and temperature throughout the lubricant film as well as the static performance characteristics. They found that addition of nanoparticles increases the load capacity of journal bearings at any eccentricity ratio, especially in the thermo-viscous case and at high values of the eccentricity ratio. They also found that the friction force increases with an increase in concentration of nanoparticles for both isothermal and thermo-viscous cases. However, at any eccentricity ratio, both end leakage flow and attitude angle decrease with the increase in concentration of nanoparticles in both cases especially at higher eccentricity ratios.

Binu et al. [25] presented a theoretical study on static characteristics of a journal bearing operating with TiO₂ based nanolubricant considering a variable viscosity model. The predicted shear viscosities of TiO₂ based nanolubricant at different volume fractions and aggregate particle sizes were obtained using modified Krieger–Dougherty model and were found to be in good agreement with experimental shear viscosities. The modified Reynolds equation considered the Krieger–Dougherty viscosities and couple-stress effects of TiO₂ nanoparticle additives at different volume fractions and particle sizes. The obtained results revealed that an increase in load carrying capacity of 45%, in comparison with plain engine oils, for journal bearings operating with nanolubricants containing TiO₂ nanoparticle additives at a concentration of 0.01 volume fraction and couple stress parameter, ${\displaystyle \tilde{l}}=\frac{l}{C}$, of 0.03108, which corresponds to an average TiO₂ nanoparticle aggregate size of 777 nm. The analysis also revealed a 40% increase in friction force for journal bearings using TiO₂ based nanolubricants as compared to plain engine oil.

Sadabadi et al. [26] employed the CFD simulation approach to model the effect of nanofluid lubricant on load carrying capacity of highly loaded industrial journal bearings. Inorganic fullerene-like WS₂ nanoparticles were used as an additive to the base synthetic industrial oil with different weight fraction percentages. The obtained results showed that the load carrying capacity of the journal bearing can be improved up to about 20% using a practical 5.0 wt.% IF-WS₂ NP additive. In addition, the SEM images showed that the high aggregation tendency of WS₂ NPs limits the application of nanolubricant at weight fractions above 6.0 wt.%.

More recently, Kadhim et al. [27] have analyzed the TEHD performance of journal bearings with pure and Al₂O₃ nanolubricant with different volume fractions by using fluid-structure interaction with CFD technique. The relevant results obtained show that the use of Al₂O₃ nanolubricant with 5.0% of NPs increases the oil film thickness by about 18%, and the presence of the NPs in the base oil lightly decreases the maximum oil film temperature.

Dang et al [28] studied the static thermal performances of a pure elliptical bearing lubricated with NPs based mineral oils at different eccentricity ratios and journal speeds. Two types of NPs CuO and TiO₂ with 0.5, 1.0, and 2.0 wt.% concentrations have been separately added in three different viscosity grades of oils. They found that the peak pressure and load capacity increase with addition of NPs and this increase was more pronounced at higher concentrations of NPs and at higher viscosity grade oils.

Gundarneeya and Vakharia [29] conducted an excellent experimental work by using a journal bearing test rig to investigate the performance analysis of a journal bearing operating on nanolubricants formulated by adding TiO₂, CuO, and Al₂O₃ NPs as lubricant additives in ISO VG 46 oil. The experimental results revealed an increase in maximum pressure and load carrying capacity of the journal bearing by adding NPs to the base oil especially for the higher volume fraction of NPs.

Bangotra and Sharma [30] investigated the impact of surface waviness on the static performance parameters of hydrodynamic journal bearings operating with lubricants containing copper oxide (Cu0) and cerium oxide (CeO₂) NPs. The static performance parameters were calculated for different waviness numbers and various wave amplitudes with variable viscosities of nanolubricants using a relationship between the relative viscosity, temperature, and weight fraction of NPs developed from experimental results. They found that the addition of such NPs to the lubricant enhanced its viscosity which further improved the steady-state parameters of the wave bearing, namely the carrying- load capacity and the friction coefficient.

Byotra and Sharma [31] developed a finite element model to analyze the performance improvement of a journal bearing by applying the arc-shaped textures on various regions of bearing operating with and without nanolubricants. They observed that the addition of Al₂O₃ and CuO NPs increased the carrying-load capacity.

In the existing technical literature, the studies of the dynamic behavior of the poroelastic internal combustion engine crankshaft journal bearings subjected to an arbitrary force torsor and operating with non-Newtonian nanolubricants are scarce. Therefore, it is felt that there is a need to carry out the dynamic performance characteristics of such journal bearings using the micro-continuum theory of Vijay Kumar Stokes [12,13] which considers the characteristic size of nanoparticles or agglomerates of nanoparticles to the extent that the latter is neglected when using classical or traditional continuum theory. The use of this theory, also known in literature as the theory of non-Newtonian polar fluids, involves stress couples and volume couples in addition to surface and volume forces in the transport equations. The non-Newtonian character of such fluids is particularly due to the skew-symmetry of the stress tensor.

The assumed incompressible and non-Newtonian nanolubricant is then characterized using three physical constants, namely: the density ρ_nf, the dynamic or absolute viscosity μ_nf, and a constant responsible for the presence of stress couples in the lubricant denoted η whose the dimensional equation is [η] = MLT⁻¹ similar to that of momentum. The first two constants naturally depend on the concentration or the volume fraction ɸ of the nanoparticles in the basic lubricant. It is then necessary to use laws of variation of these constants with concentration such as the law of Krieger-Dougherty [32,33] currently used in literature for the calculation of the effective dynamic viscosity of the nanolubricant.

In the parametric study, we introduce a novel parameter denoted $l=\sqrt{\frac{\eta }{\mu }}$ which has a dimension of length, i.e. [l] = L . This parameter accounts for the couple stress property and can be regarded as the characteristic diameter of the fullerene NPs or the agglomerates of fullerene NPs.

Basically, the present research is concerned with the study of the combined effects of the characteristic size and concentration of inorganic fullerene-like tungsten disulphide nanoparticles (IF-WS₂ NPs) or molybdenum disulphide nanoparticles (IF-MoS₂ NPs) on the nonlinear dynamic behaviour of a gasoline engine crankshaft bearing subject to an arbitrary force torsor. It is assumed that the journal (crankshaft) currently made of a forged steel is infinitely rigid and the main bearing consists of a thin poroelastic liner made of low elasticity modulus materials like Babbitt metals fixed in a stiff housing as defined by ASTM B23-00. The simplified thin elastic liner model (TELM), derived from a semi-analytical approach based on the use of complex Kolosov and Muskhelishvili potentials [17] considering asymptotic cases for plane strain hypothesis, is retaken to determine the pressure induced deformation at the bearing liner-fluid film interface. Thus, this elasticity model avoids us the use of the 3-D finite element analysis of the bushing structure.

Reynolds equation is derived in transient conditions and modified to account for the size of NPs or the aggregates of NPs dispersed in the base lubricant, and the bearing-liner permeability property. This material property is introduced by means of the Morgan-Cameron approximation [34].

In the proposed PEHD model, the normalized Reynolds equation is first written in residual form and discretized by the finite difference method using two different computational molecules (stencils).

For an arbitrary force torsor, the hydrodynamic pressure and film thickness distributions, the squeeze film velocities, and the misalignment angular velocities are determined simultaneously by solving the discretized Reynolds equation and the equilibrium equations with the damped Newton-Raphson iterative method at each crank angle step. The crankshaft center trajectories in three sections of the main journal axis, as well as the misalignment angles, are deduced from the squeeze film velocities and the misalignment angular velocities by means of a Runge-Kutta scheme.

Note that the dynamic misalignment considered in the present investigation is due to the non-central application of the dynamic load W which therefore causes a misalignment torque T over the entire cycle, i.e. $\left\{\begin{array}{c}\hfill T_X\left(t\right)\hfill \\ \hfill T_Y\left(t\right)\hfill \end{array}\right\}={\ell }_w\left\{\begin{array}{c}\hfill W_Y\left(t\right)\hfill \\ \hfill W_X\left(t\right)\hfill \end{array}\right\}$ where T_X (t) andT_Y (t) arethetorque components, ℓ_w is the lever arm, and W_X and W_Y are the dynamic load components expressed in the frame related to the engine cylinder.

Fig. 1 Representation of a rough contact lubricated by a fluid with fullerene-like NPs.

2 Governing equations

Figure 2 shows the details of the geometric configuration of a layered engine crankshaft bearing subject to an arbitrary force torsor. A coordinate system (X, Y, Z) is fixed to the non-rotating bushing with its origin located at the bearing center. The cylindrical main journal of radius R, made of a hard steel, is assumed to be rigid and rotates at a constant angular velocity $\omega ={\theta }_c^{\text{'}}=\frac{d{\theta }_c}{dt}$ about the Z-axis.

Fig. 2 Geometry of a five degree of freedom main journal bearing, and the system of coordinates.

2.1 Modified Reynolds equation and pressure boundary conditions

Using the usual assumptions of thin film lubrication theory, and considering the nanolubricant as an incompressible non-Newtonian couple-stress fluid [12–14], the transient modified Reynolds equation for laminar flow in the connecting-rod big-end bearing with a poroelastic layer may be derived in its more general form as [18]

$$\frac{1}{R2}\frac{\partial }{\partial \theta }\left[G\frac{\partial p}{\partial \theta }\right]+\frac{\partial }{\partial z}\left[G\frac{\partial p}{\partial z}\right]=12{\mu }_{nf}\left[\frac{\omega }{2}\left(1-\frac{\frac{{\ell }_2}{{\ell }_3}cos\left({\theta }_c\right)}{\sqrt{1-{\left(\frac{{\ell }_2}{{\ell }_3}\right)}^{2}si{n}^2\left({\theta }_c\right)}}\right)\frac{\partial h}{\partial \theta }+\frac{\partial h}{\partial t}\right].$$(1a)

For the main crankshaft bearing, the modified Reynolds equation (1a) reduces to

$$\frac{1}{R2}\frac{\partial }{\partial \theta }\left[G\frac{\partial p}{\partial \theta }\right]+\frac{\partial }{\partial z}\left[G\frac{\partial p}{\partial z}\right]=12{\mu }_{nf}\left[\frac{\omega }{2}\frac{\partial h}{\partial \theta }+\frac{\partial h}{\partial t}\right]$$(1b)

where

$$G\left(h,l,k\right)=h^3-12l^{2}h+24l^{3}Tanh\left(\frac{h}{2l}\right)+12\frac{k{t}_l}{1-\beta }$$(2)

and l is the characteristic size of nanoparticle added to the base lubricant.

The last term appearing in equation (2) is incorporated by using the Morgan-Cameron approximation [24] to consider the bearing-liner permeability property k.

In equations (1) and (2), h is the film thickness expressed by equation (3a) considering the initial geometrical clearance, misalignment angles of the main journal, and the elastic deformation [35]:

$$h\left(\theta ,z,t\right)=C-\left(e_X\left(t\right)+z{\alpha }_Y\left(t\right)\right)\cos \theta -\left(e_Y\left(t\right)-z{\alpha }_X\left(t\right)\right)\sin \theta +ℒp\left(\theta ,z,t\right).$$(3a)

Its rate of change is given by

$$\frac{\partial h}{\partial t}=-{e^{\prime }}_X\cos \theta -{e^{\prime }}_Y\sin \theta -\text{z}{{\alpha }^{\prime }}_Y\cos \theta +\text{z}{{\alpha }^{\prime }}_X\sin \theta +ℒ\frac{\partial p}{\partial t},$$(3b)

ℒ being the compliance operator and $\left(\right)’=\frac{d\left(\right)}{dt}$ which is the total temporal derivative.

In equation (3a),e_X and e_Y are the eccentricity components of the main journal axis in the mid-plane section. The terms ${e^{\prime }}_X$, ${e^{\prime }}_Y$, ${{\alpha }^{\prime }}_X,\text{and}{{\alpha }^{\prime }}_Y$ appearing in equation (3b) stand for squeeze film velocities in X and Y-direction, and misalignment angular velocities about X and Y axis, respectively.

Pressure within the lubricating film must satisfy the modified Reynolds equation (1) and the following boundary conditions:

$$p\left(\theta ,z=\frac{L}{2},t\right)=0,$$(4)

$$p\left(\theta =2\pi ,z,t\right)=p\left(\theta =0,z,t\right).$$(5)

Equation (4) represents conditions at the bearing ends, while (5) expresses the periodicity of pressure in the circumferential direction.

Additionally, the Swift-Stieber conditions or simply the Reynolds conditions are needed to take into account the cavitation phenomenon (fluid-film rupture) occurring in the diverging space of the bearing clearance filled with lubricant, i.e.,

$$p\left(\theta ={\theta }_{cav},z,t\right)=0\text{and}{\frac{\partial p}{\partial \theta }}_{\left(\theta ={\theta }_{cav},z,t\right)}={\frac{\partial p}{\partial z}}_{\left(\theta ={\theta }_{cav},z,t\right)}=0$$(6)

where θ_cav is the cavitation angle at which the film rupture occurs.

If no cavitation condition is prescribed, the pressure obtained from the Reynolds equation corresponds to the Sommerfeld solution. This solution is physically unrealistic because the pressure becomes negative in the diverging part of the clearance.

Once the pressure distribution is obtained, resultant film force and moment components acting on the main journal surface are evaluated by the following relations

$$F_X={\displaystyle \underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}}{\displaystyle \underset{{\theta }_1\left(t\right)}{\overset{{\theta }_2\left(t\right)}{\int }}}p\left(-cos\theta \right)Rd\theta dz;F_Y={\displaystyle \underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}}{\displaystyle \underset{{\theta }_1\left(t\right)}{\overset{{\theta }_2\left(t\right)}{\int }}}p\left(-sin\theta \right)Rd\theta dz,$$(7)

$$M_X={\displaystyle \underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}}{\displaystyle \underset{{\theta }_1\left(t\right)}{\overset{{\theta }_2\left(t\right)}{\int }}}zp\left(-sin\theta \right)Rd\theta dz;M_Y={\displaystyle \underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}}{\displaystyle \underset{{\theta }_1\left(t\right)}{\overset{{\theta }_2\left(t\right)}{\int }}}zp\left(-cos\theta \right)Rd\theta dz.$$

2.2 Viscosity model of nanolubricants

Nanolubricants are fluids containing nanoparticles such as fullerene NPs. These nontraditional fluids, used as efficient liquid lubricants in machinery, are suspensions of particles in a base lubricant.

The effective viscosity of the nanolubricant, μ_nf, is calculated using the Krieger-Dougherty viscosity model [23].

$${\mu }_{nf}={\mu }_0{\left(1-\frac{\phi }{{\phi }_m}\right)}^{-\left[{\displaystyle \overline{\eta }}\right]{\phi }_m}$$(8)

where μ₀ is the viscosity of base fluid, φ is the volume fraction of spherical fullerene NPs in base fluid. The symbol φ_m is the maximum particle packing fraction, which varies from 0.495 to 0.54 under quiescent condition, and is approximately 0.605 at high shear rates [28], and [${\displaystyle \overline{\eta }}]$ is the intrinsic viscosity whose typical value for mono disperse suspensies of hard spheres is 2.5 [36].

Typical variation of the relative viscosity, $\frac{{\mu }_{nf}}{{\mu }_0}$, with the volume fraction φ, is presented in Figure 3 for φ_m = 0.605.

Fig. 3 Graphical representation of the Krieger-Dougherty viscosity model for a range of volume fraction and φm = 0.605.

2.3 Force balance and journal center orbits

For a misaligned journal bearing, subject to an arbitrary force torsor (Fig. 2), the equations of motion are given as follows [37,38]:

$$W_X\left(t\right)+F_X\left(e_X\left(t\right),e_Y\left(t\right),e{\text{'}}_X\left(t\right),e{\text{'}}_Y\left(t\right),{\alpha }_X\left(t\right),{\alpha }_Y\left(t\right),\alpha {\text{'}}_X\left(t\right),\alpha {\text{'}}_Y\left(t\right)\right)=m_j\left(\frac{e_X^{\text{'}}\left(t\right)-e_X^{\text{'}}\left(t-\Delta t\right)}{\Delta t}\right),$$(9a)

$$W_Y\left(t\right)+F_Y\left(e_X\left(t\right),e_Y\left(t\right),e{\text{'}}_X\left(t\right),e{\text{'}}_Y\left(t\right),{\alpha }_X\left(t\right),{\alpha }_Y\left(t\right),\alpha {\text{'}}_X\left(t\right),\alpha {\text{'}}_Y\left(t\right)\right)=m_j\left(\frac{e_Y^{\text{'}}\left(t\right)-e_Y^{\text{'}}\left(t-\Delta t\right)}{\Delta t}\right),$$(9b)

$$T_X\left(t\right)+M_X\left(e_X\left(t\right),e_Y\left(t\right),e{\text{'}}_X\left(t\right),e{\text{'}}_Y\left(t\right),{\alpha }_X\left(t\right),{\alpha }_Y\left(t\right),\alpha {\text{'}}_X\left(t\right),\alpha {\text{'}}_Y\left(t\right)\right)=I_{0X}{{\alpha }^{″}}_X+I_{0z}\omega \alpha {\text{'}}_Y,$$(9c)

$$T_Y\left(t\right)+M_Y\left(e_X\left(t\right),e_Y\left(t\right),e{\text{'}}_X\left(t\right),e{\text{'}}_Y\left(t\right),{\alpha }_X\left(t\right),{\alpha }_Y\left(t\right),\alpha {\text{'}}_X\left(t\right),\alpha {\text{'}}_Y\left(t\right)\right)=I_{0Y}\alpha \text{'}{\text{'}}_Y+I_{0z}\omega \alpha {\text{'}}_X,$$(9d)

where $\alpha \text{'}{\text{'}}_X=\left(\frac{{\alpha }_X^{\text{'}}\left(t\right)-{\alpha }_X^{\text{'}}\left(t-\Delta t\right)}{\Delta t}\right)$ $\alpha \text{'}{\text{'}}_Y=\left(\frac{{\alpha }_Y^{\text{'}}\left(t\right)-{\alpha }_Y^{\text{'}}\left(t-\Delta t\right)}{\Delta t}\right)$ are the angular acceleration components about X and Y-axis We denote T_X = ℓ _wW_Y, T_Y = ℓ _wW_X as the applied couples,m_j is the journal mass, and I_0X, I_0Y, and I_0Z are components of the inertia tensor:

$$\left[\begin{array}{ccc}\hfill I_{0X}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill I_{0Y}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill I_{0Z}\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill \frac{m_j}{12}\left(3R^2+L^2\right)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{m_j}{12}\left(3R^2+L^2\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{m_j{R}^2}{2}\hfill \end{array}\right].$$

The four equations of motion (9) have eight unknowns, namely e_X (t) , e_Y (t) , e '_X (t) , e '_Y (t) , α_X (t) , α_Y (t) , α '_X (t), and α '_Y (t) . Therefore, four more equations are required and are given below [39]

$$e_X\left(t\right)=e_X\left(t-\text{Δ}t\right)+e{\text{'}}_X\left(t\right)\Delta t,$$(10a)

$$e_Y\left(t\right)=e_Y\left(t-\Delta t\right)+e{\text{'}}_Y\left(t\right)\Delta t,$$(10b)

$${\alpha }_X\left(t\right)={\alpha }_X\left(t-\Delta t\right)+\alpha {\text{'}}_X\left(t\right)\Delta t,$$(10c)

$${\alpha }_Y\left(t\right)={\alpha }_Y\left(t-\Delta t\right)+\alpha {\text{'}}_Y\left(t\right)\Delta t.$$(10d)

Taking into account equations (10), the equations of motion (9) become

See equation (11) below.

If the journal bearing is aligned and the inertia forces are neglected, equations (11) reduce to

$$\{\begin{array}{c}\hfill W_X-F_X=0\hfill \\ \hfill W_Y-F_Y=0\hfill \end{array}.$$(12)

The coupling of equations (1) and (11) results in a nonlinear system whose unknowns are the fluid film pressure, the instantaneous translational journal velocities (e ' _X, e ' _Y) also called squeeze velocities, and the instantaneous misalignment angular velocities (α '_X, α '_Y).

The instantaneous journal center position at the mid-plane section (e_X(t), e_Y(t)) of bearing as well as the instantaneous misalignment angles (α_X (t) , α_Y (t)) are then calculated from equations (10) using the modified Euler scheme or the Runge-Kutta method of second order.

The journal eccentricity components (e_XF, e_YF) at the front bearing end (z = +L/2) and the journal eccentricity components(e_XR, e_YR)at the rear bearing end (z = −L/2) (Fig. 2) are determined by the following relation:

$$\{\begin{array}{c}\hfill e_X{\mathrm{\_\_F}}_{}\left(t\right)=e_X\left(t\right)+\frac{L}{2}{\alpha }_Y\left(t\right)\hfill \\ \hfill e_Y{F}_{}\left(t\right)=e_Y\left(t\right)-\frac{L}{2}{\alpha }_X\left(t\right)\hfill \\ \hfill e_X{R}_{}\left(t\right)=e_X\left(t\right)-\frac{L}{2}{\alpha }_Y\left(t\right)\hfill \\ \hfill e_Y{R}_{}\left(t\right)=e_Y\left(t\right)+\frac{L}{2}{\alpha }_X\left(t\right)\hfill \end{array}.$$(13)

3 Numerical treatment of the governing equations

There is no analytical solution of the transient nonlinear inverse problem governed by equations (1) and (11), thus numerical calculation using the finite difference method is used to determine the pressure distribution within the oil film. Detailed description of the procedure is given in Annex A.

4 Elasticity model: expression of the compliance factor ℒ

In the film thickness equation (3), the radial component of the elastic displacement field U = ℒ p at the fluid film-bearing liner interface is obtained from elastic analysis of the bearing-liner considered as a cylindrical shell made of low elasticity modulus materials like white metals or PEEK (Poly Ether Ether Ketone, ASTM B23-00) rigidly fixed in a stiff housing (Cast steel).

Hydrodynamic pressures developed in the fluid film to endure external forces act as normal stresses on the fluid film-bearing liner interface and produce elastic deformations of the bearing shell structure. These deformations occur in the three dimensions of space; i.e. radial, axial, and circumferential.

It has been demonstrated [17] that if the relative thickness of the bearing-liner is very small; i.e. $\frac{t_l}{R}$ ≪1, the axial and circumferential components of the deformations are negligible compared to the radial. In a such case, the radial deformation at any point U (θ, z, t) belonging to the bearing-liner interface can be expressed in terms of the pressure at that point by the relation (14):

$$U\left(\theta ,z,t\right)=\frac{\left(1+\sigma \right)\left(1-2\sigma \right)}{1-\sigma }\frac{t_l}{E}p\left(\theta ,z,t\right).$$(14)

This relationship is derived assuming plane strain hypothesis [40].

From equation (14), the compliance factor is then a scalar defined by the following expression

$$ℒ=\frac{\left(1+\sigma \right)\left(1-2\sigma \right)}{1-\sigma }\frac{t_l}{E},$$(15)

where E and σ are Young's elasticity modulus and Poisson's ratio of the homogeneous and isotropic bearing-liner, respectively. The compliance matrix [L] only contains non-zero coefficients in the principal diagonal, i.e. L_ij = ℒ δ_ij where δ_ij is the Kronecker symbol or the unit second order tensor.

The scalar compliance factor ℒ is also defined as the deformation coefficient and can be regarded as a measure of flexibility and compressibility of the bearing-liner since it involves Young's modulus and Poisson's ratio.

The simplified elasticity model derived for a thin cylindrical shell with finite dimensions and expressed by equation (14) is commonly known as the “Thin Elastic Liner Model” and has been successfully integrated by several authors [41–43] in EHD simulations of single-layered cylindrical journal bearings.

The reason for using the TELM is the easy implementation and very low CPU time requirement since very little time is needed to find the values of the radial displacement U .

However, this model is only approximate and is expected to be less accurate when a journal bearing with large shell thickness particularly made of incompressible or almost-incompressible materials σ > 0.4 is to be investigated. In such cases, the three-dimensional FEM is more suitable to generate the compliance matrix [L] of the bearing-liner structure. Hence, the radial displacement at node k should be calculated by equation (16).

$$U_k={\displaystyle \sum _{m=1}^K}{L}_{km}{p}_m{}_{}\forall k\in \left\{1,K\right\}$$(16)

where k = (j − 1) (N − 1) + i, and K = (N − 1) M. If L_km = ℒ δ_km, equation (39) becomes U_k = ℒ p_k.

5 Hydrodynamic characteristics

5.1 Side leakage flow

The side leakage flow at bearing edges $\left(z=\pm \frac{L}{2}\right)$ is obtained from the following relation:

See equation (17) below.

where θ₁(t) and θ₂(t) are the angles delimiting the active zone of the bearing.

In the above equations, h is the film thickness and ${\displaystyle \overline{w}}=\frac{1}{h}{\displaystyle \underset{0}{\overset{h}{\int }}}wdy$ is the axial mean flow velocity:

$${\displaystyle \overline{w}}=-\frac{h^3-12l^{2}h+24l^3\text{Tan h}\left(\frac{h}{2l}\right)}{12h{\mu }_{nf}}\frac{\partial p}{\partial z}.$$(18)

The axial flow velocity w as well as the circumferential velocity u are determined from integration of equations governing the flow of the nano lubricant assumed as a non-Newtonian couple-stress fluid in x and z-directions [17]

$$l^2\frac{{\partial }^{4}u}{\partial y^4}-\frac{{\partial }^{2}u}{\partial y^2}=-\frac{1}{{\mu }_{nf}}\frac{\partial p}{\partial x},$$(19a)

$$l^2\frac{{\partial }^{4}w}{\partial y^4}-\frac{{\partial }^{2}w}{\partial y^2}=-\frac{1}{{\mu }_{nf}}\frac{\partial p}{\partial z}.$$(19b)

The general solutions have the following form:

$$u\left(x,y,z\right)=A_0+A_{1}y+A_2\cosh \left(\frac{y}{l}\right)+A_3\sinh \left(\frac{y}{l}\right)+\frac{1}{2{\mu }_{nf}}\frac{\partial p}{\partial x}{y}^2,$$(20a)

$$w\left(x,y,z\right)=B_0+B_{1}y+B_2\cosh \left(\frac{y}{l}\right)+B_3\sinh \left(\frac{y}{l}\right)+\frac{1}{2{\mu }_{nf}}\frac{\partial p}{\partial z}{y}^2.$$(20b)

Using the boundary conditions:

$$u\left(x,0,z\right)=0,\frac{{\partial }^{2}u}{\partial y^2}\left(x,0,z\right)=0,u\left(x,h,z\right)=U_j=\omega R,\text{and}\frac{{\partial }^{2}u}{\partial y^2}\left(x,h,z\right)=0,$$(21a)

$$w\left(x,0,z\right)=0,\frac{{\partial }^{2}w}{\partial y^2}\left(x,0,z\right)=0,w\left(x,h,z\right)=0,\text{and}\frac{{\partial }^{2}w}{\partial y^2}\left(x,h,z\right)=0.$$(21b)

we obtain

See equation (22a),(22b) below.

5.2 Power loss

The total power loss is evaluated on the active zone of the bearing from

$$P=\left|P_1+P_2\right|$$(23)

where

In equation (24), Φ₁ and Φ₂ are the dissipation functions due to the shear stress and the couple stress effects, respectively. These two functions can be defined in hydrodynamic lubrication theory as [44]

$${\text{Φ}}_1={\mu }_{nf}\left(\frac{\partial v_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}+\frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i}\right)\approx {\mu }_{nf}\left[{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial y}\right)}^2\right],$$(25)

$${\Phi }_2=l^2{\mu }_{nf}\frac{{\partial }^2{v}_i}{\partial x_j\partial x_j}\frac{{\partial }^2{v}_i}{\partial x_k\partial x_k}\approx l^2{\mu }_{nf}\left[{\left(\frac{{\partial }^{2}u}{\partial y^2}\right)}^2+{\left(\frac{{\partial }^{2}w}{\partial y^2}\right)}^2\right].$$(26)

After integration with respect to y, we get

where $F\left(h,l\right)=h^3-12l^{2}h+24l^3\text{Tan h}\left(\frac{h}{2l}\right)$

and

$$P_2=F_X·e_{\mathrm{X\prime }}+F_Y·e_{\mathrm{Y\prime }}+M_X·{\alpha }_{\mathrm{X\prime }}+M_Y·{\alpha }_{\mathrm{Y\prime }}$$(28)

6 Solution procedure of the inverse and transient EHD lubrication problem for dynamically loaded journal bearings operating with dynamic misalignment

After selecting the input parameters of the problem, for each time step t or crankshaft rotation angle θ_c ranging from 0 to 1440 degrees (two engine load cycles), the procedure for obtaining a numerical solution to the inverse EHD problem governed by equations (1) to (6), as well as equation (9) is outlined as follows:

Step 1

Initialize the iteration counter n to zero, choose initial values of ${ϵ}_X^{\left(n\right)},{ϵ}_Y^{\left(n\right)},{\alpha }_X^{\left(n\right)},{\alpha }_Y^{\left(n\right)},$ $ϵ{\text{'}}_X^{\left(n\right)},ϵ{\text{'}}_Y^{\left(n\right)},\alpha {\text{'}}_X^{\left(n\right)},\alpha {\text{'}}_Y^{\left(n\right)}$, and the initial dimensionless pressure ${{\displaystyle \tilde{p}}}_k^{\left(n\right)}$ at each interior grid point.

Step 2

Construct the initial solution $⟨X^{\left(n\right)}⟩={\left\{X^{\left(n\right)}\right\}}^T=⟨X_1^{\left(n\right)}⋮X_2^{\left(n\right)}⋮\dots ⋮X_K^{\left(n\right)}⋮\dots ⋮X_{K+1}^{\left(n\right)}⋮\dots ⋮X_{K+4}^{\left(n\right)}⟩$ such that $X_j^{\left(n\right)}={\displaystyle \widetilde{p_j}}$ for j = 1, K, and

$$X_{K+1}^{\left(n\right)}=ϵ{\text{'}}_X^{\left(n\right)};X_{K+2}^{\left(n\right)}=ϵ{\text{'}}_Y^{\left(n\right)};X_{K+3}^{\left(n\right)}=\alpha {\text{'}}_X^{\left(n\right)},\text{and}{X}_{K+4}^{\left(n\right)}=\alpha {\text{'}}_Y^{\left(n\right)}.$$

Step 3

Calculate the dimensionless film thickness ${{\displaystyle \tilde{h}}}_k$ and its rate of change at each grid point, and effective viscosity of nanolubricant ${\displaystyle \tilde{\mu }}\left(\phi \right)=\frac{{\mu }_{nf}}{{\mu }_0}$.

Step 4

Calculate the residual functions $f_i^{\left(J\right)}$ and f_i for i = 1, K + 4.

Step 5

Generate the Jacobian matrix $E_{ij}=\partial f_i^{\left(J\right)}/\partial X_j$ and solve the linear algebraic system E_ij δX_j = − f_i for δX_j such that i, j = 1, K + 4.

Step 6

Update the solution according to the Newton-Raphson recurrence formula

$$X_i^{n+1}=X_i^n+{\omega }_{NR}\delta X_i^n$$

where ω_NR is an under-relaxation factor which ranges from 0 to 1 (i.e., 0 < ω_NR ≤ 1), and n = 0, 1, 2, 3 … , n_max.

Step 7

Extract the nodal film pressures ${{\displaystyle \tilde{p}}}_k$ from the solution ⟨Xⁿ⁺¹⟩, set all negative values of pressure to zero and calculate the corresponding radial deformations ${{\displaystyle \tilde{U}}}_k=L_{km}{{\displaystyle \tilde{p}}}_m$.

Step 8

Check convergence of the Newton Raphson algorithm: $\sqrt{{\displaystyle \sum _{i=1}^{K+4}}{f}_i^2}$ <ϵ_tol where ϵ_tol = 10⁻⁶. If the convergence test is not satisfied, increment the iteration counter n (n ← n + 1), and return to step 3 to perform another iteration while n < n_max.

Step 9

Extract the velocities of the journal center and the angular velocities of misalignment ε^′_X, ε^′_Y, α^′_X, and α^′_Y from ⟨Xⁿ⁺¹⟩ and calculate the new values of the eccentricity components and misalignment angles using the modified Euler's scheme or the Runge-Kutta method of second order, i.e.,

$${ϵ}_X\left(t+\Delta t\right)={ϵ}_X\left(t\right)+\Delta t\left(\frac{{ϵ}_X^{\text{'}}\left(t\right)+{ϵ}_X^{\text{'}}\left(t-\Delta t\right)}{2}\right)$$

and similar relations for ϵ_Y (t + Δt) , α_X (t + Δt) , andα_Y (t + Δt) where $\Delta t=\left|\frac{\Delta {\theta }_c}{6n_j}\right|$, Δθ_c being the crank angle step [degrees], and n_j is the engine rotational speed [rpm].

Step 10

Calculate the instantaneous hydrodynamic performance characteristics of the journal bearing and store the obtained numerical results for post-processing.

7 Results and discussion

7.1 Validation of the thin elastic liner model (TELM)

To confirm the validity of the TELM, a separate 3D-FE computer code using eight-node hexahedral elements has been developed to generate the compliance matrix [L] for a cylindrical shell structure. The FE formulation is obtained using the principle of virtual work which corresponds to the Galerkin weak formulation of the three-dimensional elasticity problem [45].

The bearing-liner is modeled as a cylindrical shell structure made of a homogenous and isotropic material having Young's elasticity modulus of the lead-base white metal (E = 29 GPa) [46].

Two values of Poisson's ratio have been selected to study the compressibility effect on the accuracy of the TELM.

The comparison between the more accurate 3D-FE results and results obtained from the analytical TELM (Eq. (14)) for a given pressure field (Fig. 4) confirms the accuracy of the TELM, when the shell thickness t_l is much smaller than the bearing radius R and the shell material is compressible (σ < 0.40) as depicted in Figure 5.

Fig. 4 Circumferential distribution of the hydrodynamic pressure field exerted on the Bearing liner—fluid film interface.

Fig. 5 Comparison of the radial displacement field calculated by the simplified thin elastic liner model at the bearing-liner—fluid film interface with that obtained by the 3D finite element analysis.

7.2 Validation of the computer code

There are no benchmark results available for crankshaft bearing operating with non-Newtonian couple-stress nanolubricant. Therefore, the proposed PEHD model is first compared with the published results of Paranjpe et al. [47], and Hirani et al. [48] for a gasoline engine crankshaft main journal bearing using a Newtonian fluid as lubricant. The bearing and lubricant data are reported in Table 1. Figure 6 shows the load components acting on the crankshaft main bearing at 5000 rpm obtained by interpolation with cubic splines for each one-degree crank rotation angle from data given in Table 2 to ensure convergence of the Euler's method and to faithfully represent the peak loads. The two peaks in the Y-component of the load at crank angles, θ_c, at approximately 135 degrees and 255 degrees crank angle are caused by firing of the neighboring cylinder. In Table 3, we compare the predicted values of both the greatest peak pressure and the smallest film thickness with those published by Paranjpe et al. [47] and Hirani et al. [48]. It may be deduced that the computed results lie within the scatter of the results obtained by using different approaches such that the powerful technique of FEM (FEJOB) which uses the complementarity principle for implementing Swift-Stieber or Reynolds condition at both film rupture and reformation boundaries. However, there are slight discrepancies between the results which may be attributed to the feeding conditions of the bearing that were not considered in the present calculations. The journal bearing system is supplied with lubricant through a half-circumferential groove which is machined in the mid-plane section of the bearing with a supply pressure p_s = 2.758 × 10⁵ Pa.

To further validate the correctness and accuracy of the algorithm and the computer code, the Diesel Ruston and Hornsby 6 VEB-X Mk III connecting rod big-end bearing with circumferential groove was analyzed because the data are widely available in technical literature [49,50] related to dynamically loaded bearings. The engine and bearing data are given below: crank length ℓ₂ = 0.184 m, connecting rod length ℓ₃ = 0.782 m, angular velocity of the crankshaft ω_c = −20π rad/s, length of the half-bearing L = 0.057 m, bearing diameter = 0.2032 m, bearing radial clearance C = 82.55 × 10⁻⁶ m, and dynamic viscosity of the lubricant μ = 0.015 Pa s.The values of the greatest peak pressure obtained by the proposed methodology and those obtained by various other techniques are reported in Table 4. The results show that the value of the greatest peak pressure is comparable to that obtained by the FEA.

Fig. 6 Cartesian representation of dynamic loading on the Paranjpe's main bearing [47].

7.3 Parametric study

After the successful validation, the code has been used to study the combined effects of dynamic misalignment, size and concentration of NPs; and poroelasticity on the crankshaft axis stationary orbits, the minimum film thickness calculated in rear, middle, and front sections of the main bearing. The side leakage flow, and the power loss for an entire engine cycle are also computed. All calculations have been carried out for the Paranjpe's main bearing [47] consisting of an elastic layer made of lead base white metal for which the geometric and elastic characteristics are t_l = 1 mm, E = 29 GPa, σ = 0.33. The obtained results correspond to the same dynamic loading as depicted in Figure 6.

7.3.1 Effects of dynamic misalignment

Various numerical simulations have been performed considering a misaligned case. The dynamic misalignment is due to the non-central application of the dynamic load that therefore causes a misalignment torque over the entire cycle. The graphical results are presented in Figures 7–14  for Newtonian $\left(\phi =0.,\raisebox{1ex}{$l$}\!\left/ \!\raisebox{-1ex}{$C$}\right.=0.\right)$ and non-Newtonian fluids, i. e. $\left(\phi =0.3,\raisebox{1ex}{$l$}\!\left/ \!\raisebox{-1ex}{$C$}\right.=0.3\right)$, and two values of the relative lever arm $L_W=\frac{2{\ell }_W}{L}$. The aligned case correponds to L_W = 0.

Figures 7 and 8 show the comparison of the stationary main journal bearing axis trajectories for both Newtonian and non-Newtonian couple stress cases. At the ends of the layered bearing (rear and front sections), the influence of the misalignment is clearly significant. However, at the mid-plane section of the bearing, the same shape of trajectories is obtained in both aligned and misaligned cases. On the other hand, the EHD simulations performed in Newtonian case (Fig. 7) for L_W = 10% showed that the maximum eccentricity occurs at the front section plane (Fig. 7c) which tends perilously towards the bearing radial clearance. In contrast, the EHD calculations made for the same bearing using a nanofluid as lubricant predicted that the presence of NPs leads to a contraction of the orbits at the three sections for both aligned and misaligned cases as depicted in Figure 8. Therefore, the presence of NPs or aggregates of NPs moves the crankshaft center orbit towards the bushing center. This allows the system to operate safely.

In both the Newtonian and non-Newtonian couple stress cases, the minimum film thickness occurs at one of the bearing edges, i.e. the front section. Figures 9 and 10 show the variation of the minimum film thickness as a function of crank angle for both Newtonian and non-Newtonian cases. In almost the entire cycle, the values of the minimum film thickness calcultated in the Newtonian case are smaller than those obtained in the non-Newtonian couple stress case.

At the front section of the bearing (Figs. 9c and 10c), the value of minimum film thickness is calculated as 1.56 [µm] at θ_c = 231° which is about seven times smaller than the value 10.73 [µm] at θ_c = 238° obtained in the non-Newtonian couple stress case for the misaligned case, i. e. L_W = 10%.

It can be concluded that the presence of NPs which are responsible for the couple stress in the lubricating fluid can prevent the occurrence of metal-to-metal contact at the bearing edges.

Figures 11–14 show that there is no significant influence of the dynamic misalignment on side leakage flow and power loss for the entire engine cycle in comparison to the aligned case. Furthermore, the average side leakage flow calculated in the non-Newtonian case when there is a presence of NPs is smaller than that predicted in Newtonian case as depicted in Figures 11 and 13, i.e., when the bearing is lubricated with a pure base fluid. However, the opposite trend may be observed for the average power loss but only slightly Figures 12 and 14.

Fig. 7 Effects of the dynamic misalignment parameter LW on the stationary crankshaft center orbits for the Newtonian case.

Fig. 8 Effects of the dynamic misalignment parameter LW on the stationary Journal center orbits for the non-Newtonian case.

Fig. 9 Effects of the dynamic misalignment parameter LW  on the minimum film thickness over an engine cycle for the Newtonian case.

Fig. 10 Effects of the dynamic misalignment parameter LW on the minimum film thickness over an engine cycle for the non-Newtonian case.

Fig. 11 Effects of dynamic misalignment on side leakage flow over an entire cycle.

Fig. 12 Effects of dynamic misalignment on power loss over an entire cycle.

Fig. 13 Effects of dynamic misalignment on side leakage flow over an entire cycle.

Fig. 14 Effects of dynamic misalignment on power loss over an entire cycle.

7.3.2 Effects of NPs size

Figure 15 shows the influence of the NPs size parameter, l/C, on the main journal center orbits for both Newtonian and non-Newtonian cases. We note that the size of the crankshaft axis trajectories decreases with the increasing of l/C leading to an increase of the minimum film thickness as shown in Figure 16. The orbits calculated for the various values of the size parameter have the same shape.

The effect of the size of NPs gives a slight reduction of the side leakage flow as shown in Figure 17. However, the power loss is significantly affected by the NP size over the entire engine cycle as depicted in Figure 18. For the same crank angle, the power loss calculated in Newtonian case ($\frac{l}{C}=0$) is more pronounced than that obtained in non-Newtonian case especially for the higher values of the NPs size parameter ($\frac{l}{C}=0.3$).

Fig. 15 Effects of the particle size parameter, l/C, on the stationary Journal center orbits for the deformable case.

Fig. 16 Effects of the particle size parameter, l/C, on the minimum film thickness over an engine cycle for the deformable case.

Fig. 17 Effects of the particle size parameter on side leakage flow over an engine cycle.

Fig. 18 Effects of the particle size parameter on power loss over an engine cycle.

7.3.3 Effects of NPs concentration

Figure 19 depicts for the case of aligned and layered bearing, the influence of the NPs concentration parameter φ on the stationary main journal center orbits for both Newtonian and non-Newtonian cases. One can observe that the increase of this parameter decreases the size of the crankshaft axis trajectories and increases the minimum film thickness in both cases Figure 20. This moves the shaft center orbits towards the bushing center due to the increase of the nanolubricant viscosity. The orbits calculated for the various values of the NPs concentration have the same shapes.

As expected, the increase of the dynamic viscosity due to the presence of NPs leads to an increase of the power loss and a decrease of the side leakage flow as depicted in Figures 21 and 22.

Fig. 19 Effects of the NPs concentration parameter, φ, on the stationary journal center orbits for the deformable case.

Fig. 20 Effects of the nanoparticles concentration parameter, φ, on the minimum film thickness over an engine cycle for the deformable case.

Fig. 21 Effects of the nanoparticles concentration parameter on side leakage flow over an engine cycle.

Fig. 22 Effects of the nanoparticles concentration parameter on power loss over an engine cycle.

7.3.4 Effects of the poroelasticity of the bearing liner

Figure 23 depicts the effect of the liner permeability on the shaft trajectory in the case of an aligned bearing. When the shaft is far from the bearing, the effect of permeability is not important. The pressure is not sufficient enough to allow the fluid to flow into the porous layer. As the shaft approaches the bearing the porosity of the material begins to play its role as a flow regulator. The compressed fluid can penetrate the porous layer and flow inside it. The trajectory compared to the impermeable case is larger.

Figure 24 represents the minimum film thickness evolution over an engine cycle. It confirms that for the porous liner case the minimum film thickness is lower when compared to the impermeable case. The fluid is incompressible and offers greater resistance to squeeze when the material is not porous.

The effect of a permeable layer on the leakage rate is not significant as depicted in Figure 25. This is due to the low thickness of the thin layer (1 mm) and the moderate pressure level in this layer.

As the film thicknesses are lower in the case of a permeable thin layer, the shear rates are higher, and a higher power dissipation is thus observed (Fig. 26).

Fig. 23 Effects of bearing-liner permeability parameter on the stationary Journal axis orbits.

Fig. 24 Effects of bearing-liner permeability parameter on the minimum film thickness over an engine cycle.

Fig. 25 Effects of bearing-liner permeability parameter on side leakage flow over entire cycle.

Fig. 26 Effects of bearing-liner permeability parameter on power loss over entire cycle.

8 Conclusions

On the basis of the V. K. Stokes micro-continuum theory for describing the flow of nanolubricants, the combined effects of size and concentration of NPs or aggregates of NPs on the dynamic behavior of a layered main crankshaft bearing subject to an arbitrary force torsor were investigated.

The analytical thin elastic liner model (TELM) was successfully used to determine the induced pressure and the radial deformation at the bearing liner-fluid film interface.

The Krieger-Dougherty law was included in the proposed PEHD model to account for the viscosity variation with respect to the fraction volume of NPs dispersed in the base lubricant.

A transient modified Reynolds equation was derived to take into consideration the effects of couple stresses resulting from the presence of NPs in the base oil as well as the porosity of elastic thin layer. This latter was incorporated in the pressure equation by means of the Morgan-Cameron approximation.

The damped Newton-Raphson iterative method was applied to predict the hydrodynamic pressure distribution, the squeeze film velocities, and the misalignment angular velocities. The Runge-Kutta of second order scheme was used to calculate the crankshaft center trajectories and misalignment angles.

From this investigation, the following conclusions may be drawn:

• The presence of the NPs in base oil produces higher oil-film thickness, more contracted trajectories and decreases the average side leakage flow and power loss.

• The dynamic misalignment due to the non-central application of the dynamic load produces smaller oil-film thickness at the ends of bearing (front section) and more divergent trajectories for both Newtonian and non-Newtonian cases but have no influence on the leakage flow rate and power loss.

• The minimum film thickness increases with the increase of NPs size parameter.

• In Newtonian cases, the porosity of the thin layer decreases the minimum film thickness and increases the trajectories of the main crankshaft center.

The current findings will be beneficial for automotive engine designers to improve the performance characteristics of crankshaft bearings.

Nomenclature

Symbols: Definitions, SI Units

C: Bearing radial clearance, m

E: Young's elasticity modulus of the bearing-liner, Pa

e: Operating eccentricity, $=\sqrt{e_X^2+e_Y^2}$, m

F_X (t) F_Y (t): Applied load components, N

T_X (t) , T_Y (t): Torque components, N.m

ℓ_w: Lever arm, m

m_j: Crankshaft mass per main bearing, kg

W_X (t) , W_Y (t): Carrying load components, N

n_j: Engine rotational speed, rpm

k: Permeability of the porous matrix, m²

h: Oil-film thickness, m

L: Length of bearing, m

l: Size of NP (or agglomerates of NPs), m

p: Oil-film pressure, Pa

Q_z: Side leakage flow, m³/s

R: Journal radius, m

t: Time, s

Δt: Time increment, s

t_l: Thickness of bearing-liner, m

ℓ₂: Crankshaft-arm length, m

ℓ₃: Connecting-rod length, m

e_X, e_Y: Eccentricity components of the main journal axis in the mid-plane section., m

e '_X, e '_Y: Instantaneous translational journal velocities, m/s

α '_X, α '_Y: Misalignment angular velocities, rad/s

e_XF, e_YF: Journal eccentricity components at the front bearing end, m

e_XR, e_YR: Journal eccentricity componentsat the rear bearing end, m

α '_X, α '_Y: Misalignment angular acceleration components, rad/s²

M_X, M_Y: Misalignment moment components, N.m

I_0X, I_0Y, I_0Z: components of the inertia tensor, kg.m²

U (θ, z, t): Radialdeformation of the oil film-bearing liner interface, m

φ: Volume fraction of spherical fullerenes NPs in base fluid, −

β: Percolation coefficient, −

ψ: Non-dimensional permeability, −

ω_NR: Under-relaxation factor for the Newton Raphson algorithm, −

μ₀: Dynamic viscosity of base fluid, Pa.s

${\displaystyle \tilde{\mu }}$ : Relative viscosity of nanolubricant, −

P: Power loss, W

ϵ: Eccentricity ratio, ϵ = e/C, −

η: Material constant responsible for couple-stresses, N.s

σ: Poisson's ratio of the bearing-liner, −

θ: Bearing angle, rad

△θ_c: Crank angle step, rad

θ_c: Crank rotation angle, rad

ℒ: Compliance operator, m/Pa

ρ: Mass density of lubricant, kg/m³

Abbreviations

ASTM: American Society for Testing and Materials

CFD: Computational Fluid Dynamics

CPU: Central Processing Unit

FEM: Finite element method

TELM: Thin Elastic Liner Model

IF: Inorganic fullerene

NPs: Nanoparticles

GS: Gauss-Seidel

rpm: Revolution per minute

EHD: Elasto-Hydro-Dynamic

PEHD: Poro-Elasto-Hydro-Dynamic

TEHD: Thermal Elasto-Hydro-Dynamic

PEEK: Poly-Ether-Ether-Ketone

SAE: Society of automotive Engineers

SEM: Scanning Electron Microscopy

VG: Viscosity Grade of oil

Frames

(O, X, Y, Z): Frame related to cylinder of engine

(x, y, z): Local frame for writing of Reynolds equation

Annex A: Numerical treatment of the governing equations

There is no analytical solution of the transient nonlinear inverse problem governed by equations (1) and (13), thus numerical calculation using the finite difference method is used to determine the pressure distribution within the oil film.

Figure A1 shows the (2-D) finite difference grid of fluid medium showing both double-indexed-node numbers and renumbered single-indexed nodes. As illustrated in this figure, the unique number or the global index, k, assigned to each grid-node can actually be derived by using the node's i and j indices, as follows:

$$k=\left(j-1\right)\left(N-1\right)+i;i=1,N-1;\text{and}j=1,M$$(A1)

where k ∈ {1, K}.

Fig. A1 2-D orthogonal Cartesian finite difference grid used for spatial discretization of modified Reynolds equation.

It is to be noted that equation (A1) represents just one way of converting the two indices to a single index and is, by no means, the only way. Based on this numbering pattern for a periodic bearing, the five nodes in the stencil (computational molecule) now have the indices k, k − 1, k + 1, k − N + 1, and k + N − 1 instead of the central node (i, j), the western node (i − 1, j), the eastern node (i + 1, j), the southern node (i, j − 1), and the northern node (i, j + 1), respectively.

For a non-periodic bearing, equation (A1) becomes:

$$k=\left(j-1\right)N+i;i=1,N\text{and}j=1,M$$(A2)

where N stands for the number of nodes in the circumferential direction.

The advantage of using an automatic formula, such as equation (A1) in the case of periodic bearing, for generation of the unique node number is that it allows numerical generation of the Jacobian matrix $E_{km}=\frac{\partial f_k^{\left(J\right)}}{\partial p_m}$ more easily, as will be seen later.

For simplicity, we will assume that the nodes are uniformly spaced both in the θ and z- direction such that $\Delta \theta ={\theta }_{i+1}-{\theta }_i=\frac{2\pi }{N-1}$ and $\Delta z=z_{j+1}-z_j=\frac{L}{M-1}$ where N and M are the number of nodes in the θ and z-direction, respectively.

To solve the nonlinear problem, we employ two finite difference schemes in the approximation of the modified Reynolds equation. Equation (1b) written for the main crankshaft bearing is first normalized and expressed in the following residual form:

$$f=\frac{\partial }{\partial \theta }\left({\displaystyle \tilde{G}}\frac{\partial {\displaystyle \tilde{p}}}{\partial \theta }\right)+{\left(\frac{R}{L}\right)}^2\frac{\partial }{\partial {\displaystyle \tilde{z}}}\left({\displaystyle \tilde{G}}\frac{\partial {\displaystyle \tilde{p}}}{\partial {\displaystyle \tilde{z}}}\right)-12{\displaystyle \tilde{\mu }}\left(\phi \right)\left[\frac{1}{2}\frac{\partial {\displaystyle \tilde{h}}}{\partial \theta }+\frac{1}{\omega }\frac{\partial {\displaystyle \tilde{h}}}{\partial t}\right],$$(A3)

where

$${\displaystyle \tilde{G}}\left({\displaystyle \tilde{h}},{\displaystyle \tilde{l}},\psi \right)={{\displaystyle \tilde{h}}}^3-12{{\displaystyle \tilde{l}}}^2{\displaystyle \tilde{h}}+24{{\displaystyle \tilde{l}}}^{3}Tanh\left(\frac{{\displaystyle \tilde{h}}}{2{\displaystyle \tilde{l}}}\right)+\frac{12\psi }{1-\beta },$$(A4)

$${\displaystyle \tilde{h}}=1-{ϵ}_X\cos \theta -{ϵ}_{Y}sin\theta +\frac{L}{C}{\displaystyle \tilde{z}}\left({\alpha }_{X}sin\theta -{\alpha }_{Y}cos\theta \right)+{\displaystyle \widetilde{ℒ}}{\displaystyle \tilde{p}}$$(A5)

$$\frac{\partial {\displaystyle \tilde{h}}}{\partial t}={\displaystyle \widetilde{h\text{'}}}=-{ϵ}_X^{\text{'}}\cos \theta -{ϵ}_Y^{\text{'}}sin\theta +\frac{L}{C}{\displaystyle \tilde{z}}\left(\alpha {\text{'}}_{X}sin\theta -\alpha {\text{'}}_{Y}cos\theta \right)+{\displaystyle \widetilde{ℒ}}\frac{\partial {\displaystyle \tilde{p}}}{\partial t},$$(A6)

and ${\displaystyle \tilde{\mu }}=\frac{{\mu }_{nf}}{{\mu }_0}={\left(1-\frac{\phi }{{\phi }_m}\right)}^{-\left[{\displaystyle \overline{\eta }}\right]{\phi }_m}$, $\frac{\partial {\displaystyle \tilde{h}}}{\partial {\theta }_{}}={ϵ}_X\sin \theta -{ϵ}_{Y}cos\theta +\frac{L}{C}{\displaystyle \tilde{z}}\left({\alpha }_{X}cos\theta +{\alpha }_{Y}sin\theta \right)+{\displaystyle \widetilde{ℒ}}\frac{\partial {\displaystyle \tilde{p}}}{\partial {\theta }_{}}.$

The other parameters such as β and ψ are defined in the Nomenclature.

When the correct pressure distribution is found, the residual function f is identically zero. Thus, to obtain an accurate numerical solution to the inverse isothermal PEHD problem, the residual function should be accurately discretized.

To this end, the standard second-order central difference formulas are used to approximate the residual function. For ease of discretization, we first rewrite equation (A3) as:

$$f={\displaystyle \tilde{G}}\frac{\partial 2{\displaystyle \tilde{p}}}{\partial \theta 2}+\frac{\partial {\displaystyle \tilde{G}}}{\partial \theta }\frac{\partial {\displaystyle \tilde{p}}}{\partial \theta }+{\left(\frac{R}{L}\right)}^2\left[{\displaystyle \tilde{G}}\frac{\partial 2{\displaystyle \tilde{p}}}{\partial {\displaystyle \widetilde{z2}}}+\frac{\partial {\displaystyle \tilde{G}}}{\partial {\displaystyle \tilde{z}}}\frac{\partial {\displaystyle \tilde{p}}}{\partial {\displaystyle \tilde{z}}}\right]-12{\displaystyle \tilde{\mu }}\left(\phi \right)\left[\frac{1}{2}\frac{\partial {\displaystyle \tilde{h}}}{\partial \theta }+\frac{1}{\omega }\frac{\partial {\displaystyle \tilde{h}}}{\partial t}\right].$$(A7)

Now applying the central difference approximation, for the interior nodes (i.e. 1≤i ≤ N − 1, 2 ≤ j ≤ M − 1) we get

See equation (A8) below.

where

$${{\displaystyle \tilde{G}}}_{ij}={{\displaystyle \tilde{h}}}_{ij}{3}^{}-12{{\displaystyle \tilde{l}}}^2{{\displaystyle \tilde{h}}}_{ij}+24{{\displaystyle \tilde{l}}}^3\text{Tan h}\left(\frac{{{\displaystyle \tilde{h}}}_{ij}}{2{\displaystyle \tilde{l}}}\right)+\frac{12\psi }{1-\beta },$$(A9)

$${{\displaystyle \tilde{h}}}_{ij}=1-{ϵ}_X\cos {\theta }_i-{ϵ}_Y\sin {\theta }_i+\frac{L}{C}{{\displaystyle \tilde{z}}}_j\left({\alpha }_X\sin {\theta }_i-{\alpha }_Y\cos {\theta }_i\right)+{\displaystyle \widetilde{ℒ}}{{\displaystyle \tilde{p}}}_{ij},$$(A10)

$${\frac{\partial {\displaystyle \tilde{h}}}{\partial \theta }|}_{i,j}\approx {ϵ}_X\sin {\theta }_i-{ϵ}_Y\cos {\theta }_i+\frac{L}{C}{{\displaystyle \tilde{z}}}_j\left({\alpha }_X\cos {\theta }_i+{\alpha }_Y\sin {\theta }_i\right)+{\displaystyle \widetilde{ℒ}}\left(\frac{{{\displaystyle \tilde{p}}}_{i+1,j}-{{\displaystyle \tilde{p}}}_{i-1,j}}{2\Delta \theta }\right),$$(A11)

$${\frac{\partial {\displaystyle \tilde{h}}}{\partial t}|}_{i,j}\approx -ϵ{\text{'}}_X\cos {\theta }_i-ϵ{\text{'}}_Y\sin {\theta }_i+\frac{L}{C}{{\displaystyle \tilde{z}}}_j\left(\alpha {\text{'}}_X\sin {\theta }_i-\alpha {\text{'}}_Y\cos {\theta }_i\right)+{\displaystyle \widetilde{ℒ}}\left(\frac{{{\displaystyle \tilde{p}}}_{i,j}^{\left(t\right)}-{{\displaystyle \tilde{p}}}_{i,j}^{\left(t-\Delta t\right)}}{\Delta t}\right).$$(A12)

Next, we rewrite the discrete equation (A8) in terms of the unique global index k to yield

$$f_k={{\displaystyle \tilde{G}}}_k\left[\frac{{{\displaystyle \tilde{p}}}_{k+1}-2{{\displaystyle \tilde{p}}}_k+{{\displaystyle \tilde{p}}}_{k-1}}{\left(\Delta \theta \right)2}\right]+\left(\frac{G_{k+1}-{{\displaystyle \tilde{G}}}_{k-1}}{2\Delta \theta }\right)\left(\frac{{{\displaystyle \tilde{p}}}_{k+1}-{{\displaystyle \tilde{p}}}_{k-1}}{2\Delta \theta }\right)+{\left(\frac{R}{L}\right)}^2\left\{{{\displaystyle \tilde{G}}}_k\left(\frac{{{\displaystyle \tilde{p}}}_{k+N-1}-2{{\displaystyle \tilde{p}}}_k+{{\displaystyle \tilde{p}}}_{k-N+1}}{\left(\Delta {\displaystyle \tilde{z}}\right)2}\right)+\left(\frac{{{\displaystyle \tilde{G}}}_{k+N-1}-{{\displaystyle \tilde{G}}}_{k-N+1}}{2\Delta {\displaystyle \tilde{z}}}\right)(\frac{{{\displaystyle \tilde{p}}}_{k+N-1}-{{\displaystyle \tilde{p}}}_{k-N+1}}{2\Delta {\displaystyle \tilde{z}}}\}\}\right\}-12{\displaystyle \tilde{\mu }}\left(\phi \right)\left[\frac{1}{2}{\frac{\partial {\displaystyle \tilde{h}}}{\partial \theta }|}_k+{\frac{1}{\omega }\frac{\partial {\displaystyle \tilde{h}}}{\partial t}|}_k\right],\forall k\in \text{interior}$$(A13),

$f_k={{\displaystyle \tilde{p}}}_k-{{\displaystyle \tilde{p}}}_{amb},\forall k\in \text{free boundary (bearing edges)}$(A14)

where

$${{\displaystyle \tilde{G}}}_k={{\displaystyle \tilde{h}}}_k^3-12{{\displaystyle \tilde{l}}}^2{{\displaystyle \tilde{h}}}_k+24{{\displaystyle \tilde{l}}}^3\text{Tanh}\left(\frac{{{\displaystyle \tilde{h}}}_k}{2{\displaystyle \tilde{l}}}\right)+\frac{12\psi }{1-\beta },$$