Elasto-hydrodynamic lubrication analysis of a porous misaligned crankshaft bearing operating with nanolubricants

. In this paper, the combined effects of the characteristic size and concentration of inorganic fullerene-like tungsten disulphide nanoparticles (IF-WS 2 NPs) or molybdenum disulphide nanoparticles (IF-MoS 2 NPs) on the nonlinear dynamic behaviour of a gasoline engine crankshaft bearing subject to an arbitrary force torsor (effective applied force and moment vector) are theoretically and numerically investigated using the V. K. Stokes micro-continuum theory. These NPs arethe most common additives for lubrication purposes dueto their excellent tribological characteristics along with their effect on reducing friction and wear. It is assumed that the journal (crankshaft)currentlymadeofaforgedsteelisrigidandthemainbearingconsistsofathinporoelasticlinermadeof lowelasticmodulusmaterialslikeBabbittmetals ﬁ xed in a stiff housing as de ﬁ ned by ASTM B23-00. The Krieger-Dougherty law is included in the proposed EHD model to account for the viscosity variation with respect to the volume fraction of nanoparticles dispersed in the base lubricant. On the other hand, the characteristic size of nanomaterials is introduced by a new material entity, denoted l , which is responsible for a couple-stress property. The Reynolds equation is derived in transient conditions and modi ﬁ ed to account for the size of nanoparticles and the bearing-liner permeabilityproperty. For anarbitraryforce torsor, the hydrodynamicpressuredistribution, the squeeze ﬁ lm 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 misalignmentanglesarededucedfromthesqueeze ﬁ lmvelocitiesandthemisalignmentangularvelocitiesbymeans of a Runge-Kutta scheme. According to the obtained results, the combined effects of the size and concentration of fullerene-like nanoparticles on the dynamic behavior of a compliant dynamically loaded crankshaft bearing operating with dynamic misalignment are signi ﬁ cant and cannot be overlooked.


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 antipollution 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][4][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 nanoobjects 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 2 and IF-WS 2 , 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 2 and the WS 2 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][13][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][16][17][18][19][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 2 , and Al 2 O 3 ). 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 2 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 2 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 2 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 2 ), and aluminum oxide (Al 2 O 3 ). 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 2 based nanolubricant considering a variable viscosity model. The predicted shear viscosities of TiO 2 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 2 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 2 nanoparticle additives at a concentration of 0.01 volume fraction and couple stress parameter,l ¼ l C , of 0.03108, which corresponds to an average TiO 2 nanoparticle aggregate size of 777 nm. The analysis also revealed a 40% increase in friction force for journal bearings using TiO 2 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 2 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 2 NP additive. In addition, the SEM images showed that the high aggregation tendency of WS 2 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 2 O 3 nanolubricant with different volume fractions by using fluidstructure interaction with CFD technique. The relevant results obtained show that the use of Al 2 O 3 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 2 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 2 , CuO, and Al 2 O 3 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 2 ) 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 2 O 3 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 r nf , the dynamic or absolute viscosity m nf , and a constant responsible for the presence of stress couples in the lubricant denoted h whose the dimensional equation is [h] = MLT À1 similar to that of momentum. The first two constants naturally depend on the concentration or the volume fraction f 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 ¼ 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 2 NPs) or molybdenum disulphide nanoparticles (IF-MoS 2 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.
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.
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 v ¼ u 0 c ¼ du c dt about the Z-axis.

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][13][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] 1 For the main crankshaft bearing, the modified Reynolds equation (1a) reduces to 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]: Its rate of change is given by L being the compliance operator and ðÞ' ¼ dðÞ 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 0 X , e 0 Y , a 0 X ; and a 0 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: 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., where u 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

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, m nf , is calculated using the Krieger-Dougherty viscosity model [23].
where m 0 is the viscosity of base fluid, f is the volume fraction of spherical fullerene NPs in base fluid. The symbol f 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 [h is the intrinsic viscosity whose typical value for mono disperse suspensies of hard spheres is 2.5 [36]. Typical variation of the relative viscosity, m nf m 0 , with the volume fraction f, is presented in Figure 3 for f m = 0.605.

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]: ð9aÞ ð9bÞ are the angular acceleration components about X and Y-axis We denote T X = ℓ w W Y , T Y = ℓ w W 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: The four equations of motion (9) have eight unknowns, namely e X (t) , e Y (t) , e ' X (t) , e ' Y (t) , a X (t) , a Y (t) , a ' X (t), and a ' Y (t) . Therefore, four more equations are required and are given below [39] 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 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 (a ' X , a ' Y ).
The instantaneous journal center position at the midplane section (e X (t), e Y (t)) of bearing as well as the instantaneous misalignment angles (a X (t) , a 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: 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.

Elasticity model: expression of the compliance factor L
In the film thickness equation (3), the radial component of the elastic displacement field U = L p at the fluid filmbearing 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 filmbearing 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. 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 (u, z, t) belonging to the bearing-liner interface can be expressed in terms of the pressure at that point by the relation (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 where E and s are Young's elasticity modulus and Poisson's ratio of the homogeneous and isotropic bearing-liner, respectively. The compliance matrix [L] only contains nonzero coefficients in the principal diagonal, i.e. L ij = L d ij where d ij is the Kronecker symbol or the unit second order tensor. The scalar compliance factor L 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][42][43] in EHD simulations of singlelayered 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 s > 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).
where u 1 (t) and u 2 (t) are the angles delimiting the active zone of the bearing.
In the above equations, h is the film thickness and w ¼ 1 h ∫ h 0 wdy is the axial mean flow velocity: 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] The general solutions have the following form: we obtain See equation (22a),(22b) below.

Power loss
The total power loss is evaluated on the active zone of the bearing from where In equation (24), F 1 and F 2 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] After integration with respect to y, we get 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 u 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: : (24) (22a) Step 1 Initialize the iteration counter n to zero, choose initial values of e n ð Þ X ; e n ð Þ Y ; a n ð Þ X ; a n ð Þ Y ; e 0 X n ð Þ; e 0 Y n ð Þ; a 0 X n ð Þ; a 0 Y n ð Þ, and the initial dimensionless pressurep n ð Þ k at each interior grid point.

Step 2
Construct the initial solution ⟨X n ð Þ ⟩ ¼ X n Step 3 Calculate the dimensionless film thicknessh k and its rate of change at each grid point, and effective viscosity of Calculate the residual functions f J ð Þ i and f i for i = 1, K + 4.
Step 5 Generate the Jacobian matrix E ij ¼ ∂f J ð Þ i =∂X j and solve the linear algebraic system E ij dX j = À f i for dX j such that i, j = 1, K + 4.

Step 6
Update the solution according to the Newton-Raphson recurrence formula where v NR is an under-relaxation factor which ranges from 0 to 1 (i.e., 0 < v NR 1), and n = 0, 1, 2, 3 … , n max .
Step 7 Extract the nodal film pressuresp k from the solution ⟨X n+1 ⟩, set all negative values of pressure to zero and calculate the corresponding radial deformations U k ¼ L kmpm .
Step 8 Check convergence of the Newton Raphson algorithm: <e tol where e tol = 10 À6 . 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 e 0 X , e 0 Y , a 0 X , and α 0 Y from ⟨X n+1 ⟩ 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., and similar relations for e Y (t + Dt) , a X (t + Dt) , and a Y (t + Dt) where Dt ¼ j Du c 6n j j, Du 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.

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 (s < 0.40) as depicted in Figure 5.

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, u 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   Table 2. Numerical values of the dynamic loading on the crankshaft main bearing for one engine cycle at 5000 rpm [47,48].  Hirani et al. [48] 34.66 Booker's short bearing approximation [51] 35.84 Goenka's new curve fits [52] 34.57 Finite element method [52] 34.40 Present study 35.82 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 5 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] Table 4. The results show that the value of the greatest peak pressure is comparable to that obtained by the FEA.

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

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 f ¼ 0:; l C ¼ 0: = Þ À and non-Newtonian fluids, i. e. f ¼ 0:3; l C ¼ 0:3 = Þ À , and two values of the relative lever arm L W ¼ 2ℓ 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 [mm] at u c = 231°which is about seven times smaller than the value 10.73 [mm] at u 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. 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.

Effects of NPs size
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 ( l C ¼ 0) is more pronounced than that obtained in non-Newtonian case especially for the higher values of the NPs size parameter ( l C ¼ 0:3).    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.

Effects of the poroelasticity of the bearing liner
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).

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.

Symbols
Definitions, SI Units C Bearing radial clearance, m E Young's elasticity modulus of the bearing-liner, Pa e Operating eccentricity, ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  (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: where k ∈ {1, K}. 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: 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 ¼ ∂f J ð Þ k ∂p m more easily, as will be seen later.
For simplicity, we will assume that the nodes are uniformly spaced both in the u and z-direction such that Du ¼ u iþ1 À u i ¼ 2p NÀ1 and Dz ¼ z jþ1 À z j ¼ L MÀ1 where N and M are the number of nodes in the u 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: whereGh ;l; c À Á ¼h 3 À 12l 2h þ 24l 3 Tan hh 2l andm ¼ m nf The other parameters such as b and c 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: 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.
wherẽ G ij ¼h ij 3 À 12l 2h ij þ 24l 3 Tan hh ij 2 l þ 12c 1 À b ; ðA9Þ (A8) ∂h ∂u j i;j ≈ e X sin u i À e Y cos u i þ L Cz j a X cos u i þ a Y sin u i ð Þ þLp iþ1;j Àp iÀ1;j 2Du ; ðA11Þ Next, we rewrite the discrete equation (A8) in terms of the unique global index k to yield f k ¼p k Àp amb , ∀k ∈ free boundary (bearing edges) (A14) whereG k ¼h 3 k À 12l 2h k þ 24l 3 Tanhh k 2 l þ 12c 1 À b ; h k ¼ 1 À e X cos u k À e Y sin u k þ L Cz a X sin u k À a Y cos u k ð Þ þ Lp k ; ∂h ∂u j k ≈ e X sin u k À e Y cos u k þ L Cz k a X cos u k þ a Y sin u k ð Þ þLp kþ1 Àp kÀ1 2Du ; In all, we have to solveK nonlinear equations where = (N À 1) M . To these K equations, we must add the following equations also written in residual form: Recall that k = (j À 1) (N À 1) + i such that i = 1, N À 1 andj = 2, M À 1.
The computational molecule shown in Figure A2 is used to approximate the residual function f, from which the Jacobian matrix is generated. Thus, the discrete form of equation (A3) reads as