A CFD-Based Frequency Response Method Applied in the Determination of Dynamic Coefficients of Hydrodynamic Bearings . Part 1 : Theory

A general, CFD-based frequency response method for obtaining the dynamic coefficients of hydrodynamic bearings is presented. The method is grounded in experimental parameter identification methods and is verified for an extremely long, slider bearing geometry as well as short and long journal bearing geometries. The influence of temporal inertia on the dynamic response of the bearings is discussed and quantified through the inclusion of added mass coefficients within the mechanical models of the hydrodynamic bearing films. Methods to separate the dynamic stiffness into static stiffness and added mass contributions are presented and their results compared. Harmonic perturbations are applied to the bearings at varying frequencies to determine the frequency dependence of the dynamic coefficients and to facilitate the decomposition of the dynamic stiffness into its constituents. Added mass effects are shown to be significant for the extremely long slider bearing geometry and negligible for the short and long journal bearing geometries under operating conditions motivated by those typical of marine bearings.


Introduction
In their review paper, Dimond et al. [1] presented methods for the identification of bearing dynamic coefficients (stiffness, damping, and added mass) based on experimentally measured data.They discussed bearing coefficient identification employing force excitation using incremental loading, impulse, pseudo-random, and unbalance excitation applied to a moving shaft and fixed bearing housing, or viceversa.Goodwin [2] also reviewed experimental techniques used for the measurement of bearing coefficients.He reported that typically there was a 20% discrepancy between the measured and theoretical results.In his review, Goodwin covered static loading, shakers, unbalance forces, and transient methods.Tiwari [3] offered a survey of methods applied in the determination of bearing and seal dynamic coefficients that covered fixed geometry and tilting pad bearings as well as rolling element bearings, active magnetic bearings, and hydrodynamic seals.Someya [4] assembled an exhaustive quantity of hydrodynamic journal bearing data (Sommerfeld number and dynamic coefficients with respect to eccentricity ratio) in addition to results obtained numerically through solution of the static and perturbed forms of Reynolds equation.
Kanki et al. [5] used a special testing system to determine the dynamic characteristics of water lubricated bearings by applying sinusoidal excitation forces in the x-and y-directions (perpendicular to the rotor axis) by means of hydraulic actuators and measuring the journal displacement to determine the transfer function.They used two methods: (a) the first provided for the application of excitation force with phase difference of 90 deg in the x and y directions and (b) by applying the excitation forces at two different frequencies in the x and y directions simultaneously.The twelve dynamic coefficients were determined by solving the equation of motion for several different frequencies and by curve fitting the obtained transfer function.Diana et al. [6] also obtained dynamic coefficients using two sinusoidal forces applied at dissimilar excitation frequencies.The vertical and horizontal forces generated in the film were measured, and the bearing dynamic coefficients were identified through a frequency-domain approach.
Xu [7] studied a hybrid hydrostatic bearing both from the point of view of pressure performance and the dynamic coefficients.An experimental rig was outfitted to allow for axial misalignment and thus verification of the numerical model for a multitude of circumstances.The dynamic bearing coefficients were determined by means of a frequency response function with excitations at two frequencies to yield enough equations to determine the 8 dynamic coefficients.
Ha and Yang [8] performed experiments using a tilting pad bearing to assess the effect of excitation force frequency on the linear stiffness and damping coefficients of this type of bearing.The excitation frequency was controlled independently for each of the hydraulic shakers mounted at 90 deg from each other.The acceleration of the bearing bushing as well as the relative motion between the bushing and the journal were measured with respect to the excitation frequency.The authors found that the stiffness and damping coefficients only changed slightly with the excitation frequency, which was varied using slightly subsynchronous values (down to 0.5× rotational speed).
Ono [9] studied analytically/numerically the dynamic characteristics and response of an air lubricated slider bearing for use in non-contact magnetic recording storage.The dynamic pressure distribution of the air film was calculated in the frequency domain starting with, and using a non-dimensional Reynolds equation for compressible fluid of the form: where σ = 12ωL 2 µ p re f h 2 0 is the squeeze, or frequency number.The frequency domain analysis was performed by means of a perturbation method applied both to the equations of motion of the slider and the compressible Reynolds equation.The frequency number played a vital role in the subsequent dynamic (identification) analysis.An experimental study accompanied the theoretical/numerical development.
The author concluded that the stiffness and damping in both translation and rotation were frequency dependent, but became almost constant for perturbation frequencies below u/2L, with u and L represent the sliding velocity and the slider length, respectively.Parsell et al. [10] showed that for a tilting-pad bearing the frequency of the excitation force can be an important parameter with regard to stiffness and damping characteristics for high-speed bearings.The authors found stiffness and damping coefficients to be insensitive to excitation frequency except in the absence of bearing preload, or for large values (>10) of Sommerfeld number, S = µN/P b (R/C) 2 , where N is rotations per seconds and P b is the projected bearing pressure.When very lightly loaded and/or running at high speeds, the dynamic coefficients exhibited strong dependence on the excitation frequency, which was varied over a range of subsynchronous values.The authors' results were presented in non-dimensional form in absence of the lubricant type.
Franchek and Childs [11] presented experimental results for four different high speed, high pressure, orifice compensated hybrid bearing designs lubricated with water.Due to the high speed of rotation, added mass terms due to temporal and advective inertia were considered in the determination of the dynamic coefficients.The authors used hydraulic shakers to alternately deliver, in the xand y-directions respectively, priorly defined, pseudo-random excitations (Bendat and Piersol [12]).The excitation frequency spanned a range of 40-440 Hz.The input force and output bearing acceleration were cast in the frequency domain using a Power-Spectral-Density technique.The real and imaginary components of the resulting frequency response functions (FRFs) were curve fit to obtain the stiffness, damping, and added mass coefficients [13].
Al-Ghasem and Childs [14] investigated the dynamic properties of a high speed tilting pad bearing lubricated with ISO VG 32 turbine oil.The authors found that the direct dynamic-stiffness coefficients were frequency dependent as a result of pad and fluid inertia.Fluid inertia effects were accounted for through the introduction of added mass coefficients within the dynamic model of the bearing.Experimentally determined, and bulk-flow predicted, added mass coefficients were found to be comparable in magnitude.
Delgado et al. [15] presented the identification of rotordynamic coefficients for four-and five-pad tilting pad bearings lubricated with oil.The dynamic coefficients were identified by applying multiple frequency excitations varying between 20 and 300 Hz, and exerted independently through two hydraulic shakers arranged at 90 deg from each other.The direct dynamic coefficients were identified from curve fits of the complex dynamic stiffness.The direct stiffness and damping coefficients were found to be frequency independent if an added mass term, smaller than the test device modal mass, was introduced.Both the parameter identification procedure and the experimental technique were based on the earlier work of Rouvas and Childs [16] with water-lubricated bearings.Numerical results were found to be within 50% for direct stiffness and 30% for direct damping, respectively, of the experimental results.
Zhang et al. [17] numerically evaluated the dynamic coefficients of a water-lubricated journal bearing based on the steady-state solution of the Navier-Stokes (NS) equations.The CFD model employed ANSYS Fluent to calculate the pressure profiles while taking into consideration both fluid inertia and cavitation through the Zwart cavitation model.The authors determined the bearing stiffness coefficients by curve-fitting the steady-state forces over a range of prescribed static journal eccentricities.The stiffnesses obtained using CFD were compared with Constantinescu's [18] with good coincidence.
Certainly there is a wealth of experimental research that concerns the identification of the dynamic coefficients and a measured amount of theoretical studies that propose and identify methodologies for the numerical prediction of these coefficients.This brief introduction cannot, unfortunately, do justice by covering all such papers.However, for the reader further interested in methods for identification of the dynamic coefficients that use (i) incremental loading, (ii) sinusoidal excitation, (iii) transient excitation or (iv) pseudo-random excitation, the review papers previously mentioned by Dimond et al. [1] and Tiwari et al. [3], and Goodwin [2] are strongly recommended.Tiwari et al. [19] also presents in its Table 1 an exhaustive list of the experimental investigations, their authors, and the methods to identify the dynamic coefficients of bearings and seals.

Scope of Work
Snyder and Braun [20] previously applied a frequency domain method to CFD-predicted hydrodynamic forces to determine the dynamic coefficients of bearings.The method was shown to be a viable alternative to solving the perturbed Reynolds equation for bearing dynamic coefficients.Unlike the perturbed Reynolds equation, CFD more easily accommodates complex flow physics such as inertia, turbulence, thermal effects, and multiphase flows and with less empiricism and model tuning.Moreover, the frequency domain approach, borrowed from experimental methods to determine dynamic coefficients of bearings, can be applied in Fluid Structure Interaction (FSI) problems in which the hydrodynamic fluid film interacts with a solid, such as the deforming rubber liner in a marine bearing (E = 1 − 10 [MPa], small strain elastic modulus approximated from indentation hardness [21]).It is noteworthy, however, that in "traditional" bearings, solid deformation primarily results from thermal strain rather than mechanical strain induced from the hydrodynamic pressure [22].The term "traditional" is used here in reference to bearing with conformal geometry and in which the elastic modulus of the bushing is large compared with rubber, e.g., E = 115 [GPa] for bronze (leaded, commercial grade [23]).
As mentioned, the CFD-based frequency response method follows from experimental methods which perturb bearings or seals with electromechanical shakers, Active Magnetic Bearings (AMBs), or rotating unbalance.The notable difference of the proposed method is that the bearing forces are resolved numerically using CFD rather than obtained from experiment.
The present work is a follow-up to Snyder and Braun [20], and extends the authors' previous theoretical work in two main ways: 1.
Alternative methods to separate the dynamic stiffness into added mass and static stiffness effects are investigated.The authors [20] previously demonstrated that temporal inertia effects as embodied in added mass coefficients are inherent within transient CFD simulations.It then becomes necessary to perturb the bearing at multiple frequencies in order to fit the static stiffness and added mass parameters.However, different techniques for fitting the dynamic stiffness will provide different estimates of the coefficients depending on assumed behavior of the coefficients.Details regarding the implementation and robustness of these techniques are explored herein.

2.
The proposed methodology is extended to hydrodynamic journal bearing geometries.These bearings are generally of more practical interest than simple slider bearings and an extension of the CFD-based frequency response method to determine the dynamic coefficients is a non-trivial task, warranting a thorough explanation.The dynamic coefficients predicted for short L/D = 0.25 and long L/D = 3.0 journal bearing geometries are compared directly with those obtained from the perturbed Reynolds equation.

Linear Slider Bearing
Figure 1 depicts the extremely long linear slider bearing geometry (length L in the z-direction is considered infinite) and a mechanical model of the hydrodynamic fluid film.The geometric and operational parameter values are summarized in Table 1 and are motivated by the physical characteristics of rubber-lined, marine bearings where individual staves create a local geometry similar to that of a series of extremely long slider bearings assembled in the circumferential direction.The lubricant considered is water and the minimum film thickness, h 2 = 266.7 [µm], is based on the radial clearance of a marine journal bearing with a diameter of 50.8 [mm].For a traditional (oil-lubricated, plain) hydrodynamic journal bearing, the ratio of the radial clearance and the bearing radius, C R , is generally between 1/500 and 1/1000.For marine bearings, C R is slightly larger (3.4/500) [24].The nominal bearing width b = 39.22 [mm] satisfies h 2 b = 3.4 500 motivated by marine bearing dimensions.Nominally, a bearing taper ratio r = 2.0 is considered, where r = h 1 /h 2 .
The bearing, inflow, and outflow surfaces s i (where i = 1, 2, 3, 4) will be subsequently used for identification of CFD boundary conditions.

Parameter
Value Unit The sliding velocity u = 1.88 [m s −1 ] results from fixing the sliding Reynolds number Re u = 500 (Re = uh Perturbations are applied to bearing surface s 1 with an amplitude A and and frequency ω, which is varied.The oscillation amplitude is taken to be 10% of the minimum film height h 2 .

Journal Bearing
The geometry and mechanical model of the hydrodynamic fluid film of a plain journal bearing is depicted in Figure 2. The journal diameter is D and the radial clearance of the bearing is C. The bearing length is L in the z-direction (not shown).The rotational speed of the bearing is Ω.The film thickness of the bearing h is defined through the angular coordinate θ which in-turn is defined by a line extending along the bearing center O b and the journal center O j (line-of-centers).The reference frame (X, Y) is fixed (inertial) with its origin located at the bearing center.The reference frame (R, T) is also fixed (inertial) with its origin located at O j (O j is translated a distance e from O b ) and rotated through an angle φ from (X, Y).The position and orientation of reference frame (R, T) is defined in terms of the bearing's static equilibrium position, which may be expressed as (x 0 , y 0 ) or (e cos φ, e sin φ) in reference frame (X, Y).
The geometric and parameter values of the journal bearing considered in the present work are summarized in Table 2.

Parameter
Value Unit The bearing clearance C was determined from the nominal value of C/R equal to 0.0068 for a marine bearing and a bearing radius of 25.4 [mm] (i.e., 1 inch).The bearing Reynolds number Re = ΩRC/ν is equal to 323.9.As the linear slider bearing, the flow may be considered laminar, considering Re = 2000 as the onset for turbulence [25].Both short and long bearing geometries were investigated, the latter having a L/D = 3.0, which is typical of marine bearings.

Linear Slider Bearing
A linear mechanical model of the hydrodynamic fluid film dynamic response can be seen in Figure 1.The hydrodynamic film is modeled through a linear spring, damper, and added mass.Applying Newton's second law to the upper surface s 1 considering a small displacement y from the static equilibrium position results in Equation (1) for the dynamic bearing load w y,d (t).
In Equation ( 1), the mass of the slider has been neglected and slider is constrained to move only in the y-direction and for which the initial position of the upper surface of the bearing s 1 is considered zero, y = 0.A notable feature of Equation ( 1) is retention of an added mass coefficient associated with inertia of the fluid.Snyder et al. [20] previously demonstrated that lubricant inertia can cause an appreciable difference between the dynamic bearing stiffness (k yy + m yy ω 2 ) predicted using CFD (lubricant inertia present) and the static stiffness (k yy ) predicted by the perturbed Reynolds equation (lubricant inertia neglected).
In the present work, the dynamic coefficients present in Equation ( 1) are determined by correlating an applied harmonic displacement (input) with the CFD-predicted force (output).Correlation of the input and output is performed both in the time-domain through curve fitting and in the frequency domain through development of a frequency response function (FRF).
For the time-domain approach, consider the application of a known input displacement (sinusoidal) y = −A sin (ωt) to the bearing surface s 1 depicted in Figure 1.The sinusoidal displacement may alternatively be expressed as a complex exponential function as y(t) = A exp {−iωt}.If the input displacement is a harmonic sinusoidal function, the steady-state output force of this linear model will also be a sinusoidal function which can be written as a complex exponential function w y,d (t) = F exp {−i(ωt + α)}, where F is the unknown force amplitude and α is the phase shift of the harmonic force.Substituting the complex exponential forms of the displacement and force into Equation (1) leads to Equations ( 2) and (3).
Alternatively, the dynamic stiffness and damping coefficients can be determined using a frequency domain approach through the establishment of a frequency response function (FRF), H yy (−iω) as in Equation ( 4).
The dynamic stiffness and damping coefficients are then determined from the real and imaginary components of H yy as given by Equations ( 5) and (6).
In order to separate the dynamic stiffness into static stiffness and added mass components, the perturbation frequency ω is varied.Several techniques are investigated to separate the dynamic stiffness within the Results section.

Journal Bearing
The journal depicted in Figure 2 is subject to a dynamic load w(t) and the net force generated by the hydrodynamic film is f(t).If the film force f(t) is linearized through a Taylor Series expansion and expressed in terms of linearized dynamic coefficients (stiffness, damping, and added mass), the Equations of motion may be expressed in matrix form as Equation (10).
In Equation (11), M and M j represent the added mass (fluid inertia) matrix and the journal mass matrix, respectively.Unlike in an experiment, a massless journal (M j = 0) can be simulated using CFD and in conjunction with Equation ( 11) can be reduced to Equation (12).
If a sinusoidal displacement of complex form q = X(iω) exp (iωt) is applied to the journal, then the corresponding output force will also be sinusoidal, w d (t) = W(iω) exp (iωt).X(iω) and W(iω) are the complex amplitude vectors of the input displacement and the output force, respectively.Substituting these expressions into Equation ( 12), results into Equation (13).
The dynamic stiffness and damping coefficients can be evaluated at different perturbation frequencies ω i using Equations ( 15) and ( 16), respectively.
The preceding development involved the determination of the the dynamic coefficients in the fixed (X, Y) reference frame of Figure 2.However, results from the CFD model are easiest to obtain through the application of a static eccentricity along a single coordinate direction.Setting up the CFD model in this manner means that the dynamic coefficients which result will be in the (R, T) reference frame of Figure 2. A coordinate transformation is necessary to move from (O j , R, T) to (O b , X, Y).Equations ( 17) and (18) summarize the transformation of the stiffness matrix from the (X, Y) to (R, T) and from the (R, T) to (X, Y) reference frames, respectively.Identical transformations also apply to the added mass and damping matrices.
where the rotation matrix is given by The stiffness, damping, and added mass matrices will be presented in dimensionless for as given in Equations ( 20)- (22), respectively [4].
In Equations ( 20)- (22), Ω is the shaft rotational speed in radians per second, and W 0 is the static load carrying capacity of the bearing.

Hydrodynamic Fluid Film
The hydrodynamic fluid film is resolved by solving an advectiveless form ((v • ∇)v = 0) of the Navier-Stokes equations given in Equations ( 23) and ( 24) under the assumptions of a Newtonian, incompressible, isothermal fluid.
In the absence of the temporal inertia term, ∂v ∂t , the resulting set of equations are known as Stokes equations and were used to analyze the flow around spherically shaped pendulums [27] in creeping flow (Re → 0).The Stokes equations are at the basis of the genesis of the Reynolds equation, which is commonly employed in hydrodynamic lubrication.The difference resides in the fact that while the Reynolds equation is integrated across the film, and thus rendered two-dimensional, Stokes equations preserve their 3D character.
In hydrodynamic bearings resolved using CFD with ever-present temporal inertia, it is necessary to take the inertia effects into consideration through the addition of added mass terms within the mechanical model of the fluid film.The added mass elements may or may not significantly influence the dynamic stiffness depending on the bearing geometry and operating conditions.If significant, the added mass terms must be separated from the dynamic stiffness to determine the static stiffness elements.

CFD Solver
Steady and transient forms of the Navier-Stokes and Stokes equations are solved using the open-source continuum mechanics library OpenFOAM [28] version 5.0.The three dimensional form of the equations are discretized on a boundary conforming, non-overlapping computational mesh using the finite volume method.The meshing application blockMesh is used to generate the mesh using hexahedral cells."Infinitely long" models are achieved by spanning the z direction of the computational domain with a single cell in conjunction with the specification of "empty" boundary conditions (OpenFOAM terminology) on the associated axial boundaries.
The OpenFOAM application applied in the analysis is termed "pimpleDyMFoam" which solves the incompressible form of the Navier-Stokes equations (modified with removal of advective inertia) on a dynamic (moving) mesh.The solver is a pressure-based, predictor-corrector algorithm based on a combination of PISO [29] and SIMPLE [30,31] algorithms.For details regarding the specific discretization schemes and solver parameters, the interested reader is referred to [20].

Boundary Conditions
The application of the various boundary condition types to specific surfaces (patches) within the CFD models are summarized in Tables 3 and 4 for the linear slider bearing and the journal bearing, respectively.
Table 3. Linear slider bearing OpenFOAM boundary conditions.See Figure 1 for model surfaces identification names.Boundary conditions values are required for the velocity U, the pressure P, and displacement D fields directly (Dirchlet) or their gradients (Neumann).For the extremely long model, the "Front" and "Back" (not included in table) are specified as "empty".

Journal Bearing
The journal bearing CFD model is setup with a prescribed (fixed) eccentricity in the r direction in the coordinate system as introduced in Figure 2. Note that (R, T) coordinate system as introduced in Figure 2 corresponds to an (X, Y) reference system of the CFD simulation.Also, note that within the CFD model, journal rotation is not applied directly, but rather, rotational velocity is applied in the opposite direction to the bearing.Alternatively, if the journal were supplied with rotational velocity and oscillating displacement, a new boundary condition would have to be implemented within OpenFOAM.

Linear Slider Bearing
Verification of the CFD model of the linear slider bearing requires consideration of three primary factors: spatial convergence, local temporal convergence, and global temporal convergence.Spatial convergence refers to a mesh density which adequately resolves the discrete approximations of the spatial derivatives in the model (differential) equations.Local temporal convergence refers to a time step which adequately resolves the discrete approximations of the temporal derivatives in the model (differential) equations.Spatial and local temporal convergence are typical metrics associated with the accuracy of CFD results.The additional consideration of global temporal convergence is required as the dynamic CFD results will be fit to a model.The global temporal convergence is related to the number of oscillating periods of the boundary perturbation.
For all three factors, the dynamic stiffness and the damping coefficients are metrics of interest.Spatial and temporal convergence is quantified for a single perturbation frequency of 48 [rad•s −1 ].Tables 5-7, respectively summarized the impact of mesh density, time step size, and number of simulated periods on the predicted dimensionless dynamic coefficients.Dynamic coefficients predicted using both the time-domain curve-fit method and the frequency response function are presented.In general, the damping coefficients are more sensitive to the numerical parameters than the dynamic stiffness.Additionally, the FRF method is marginally more sensitive to the numerical parameters than the curve-fit method.Overall, using the mesh density of 221 × 15 with five periods of oscillation and 200 time steps per oscillation will result in numerical errors below 2% (for the curve-fit method).This setup for the mesh density and temporal parameters was chosen as it can be used to predict dynamic coefficients will less than 5% numerical error with reasonably low computational overhead (each case takes approximately 30 minutes using four threads on a Intel i7-4790 3.6 GHz CPU).In all cases, the perturbation amplitude was set to 10% of the minimum film thickness h 2 .The predicted dynamic coefficients remained unchanged reducing the oscillation amplitude from 10% to 5% of the minimum film thickness.Thus, for a perturbation amplitude of 10% of the minimum film thickness, the dynamic response of the hydrodynamic film is linear with respect to perturbation amplitude.
The nominal model setup with a mesh density of 221 × 15, five oscillation periods, and 200 time steps per oscillation, the curve-fit and FRF methods predict dynamic stiffness and damping coefficients which differ by 2.2% and 5.6%, respectively.Proceeding, the FRF results will be presented owing to the generality of the method and prevalence in experimental parameter identification techniques.

Journal Bearing
The short journal bearing model (Table 2 was ran with time step size of 1 ×10 −5 [s] with five periods of oscillation.Dynamic coefficients evaluated with a mesh density of 20,000 cells (10 × 100 × 20, "across-film" × "circumferential" × "axial") all fell within 2% of the those coefficients predicted using a mesh with 90,000 cells (15 × 200 × 30).For the long journal bearing, the mesh density was 10 × 100 × 120 and like the short bearing, the predicted dynamic coefficients fell within 2% of a finer, 15 × 100 × 120 mesh.

Linear Slider Bearing (Extremely Long)
The stiffness and damping coefficients predicted using CFD are compared against the analytic results yielded by Equations ( 25) and (26), respectively.
The dynamic stiffness of the linear slider bearing is plotted with respect to perturbation frequency in Figure 3.The perturbation frequency was varied from 4.8 to 480 [rad•s −1 ].The dynamic stiffness exhibits quadratic dependence on perturbation frequency consistent with the presence of an added mass contribution as evidenced by the form of the dynamics stiffness, K yy + ω 2 M yy .At the low frequency limit, the dynamic stiffness approaches the static value.Also plotted in Figure 3 are the CFD results with the linear slider bearing undergoing pure squeeze motion (no sliding velocity) and the analytic dynamic stiffness, i.e., the static stiffness.Under pure squeeze motion, the added mass effect is isolated and its significant contribution to the dynamic stiffness is highlighted.Moreover, it is apparent from the figure that the pure squeeze results can be utilized to separate the static stiffness and added mass coefficients from the dynamic stiffness.is the dimensionless perturbation frequency (frequency ratio) whose form is motivated by the dimensionless forms of of the static stiffness K yy and the added mass M yy given in Equations ( 7) and ( 9), respectively.
Several methods can be applied in estimating the added mass coefficients from the dynamic stiffness variation with perturbation frequency and are summarized below.
Method 1: The added mass coefficient is evaluated at each perturbation frequency from the pure squeeze CFD results.Using this technique, the added mass coefficient as well and as the static stiffness may exhibit frequency dependence.In this method, at each perturbation frequency, the dynamic stiffness predicted from the pure-squeeze CFD simulation is assumed to equal ω 2 M yy .Thus, the added mass coefficient is determined at the given excitation frequency.Then, this quantity is subtracted from the dynamic stiffness considering both squeeze and sliding, K yy + ω 2 M yy , leaving the static stiffness, K yy .
Method 2: The added mass coefficient and the static stiffness are estimated from direct curve-fitting of a quadratic function to the dynamic stiffness.Using this technique, the predicted static stiffness coefficient is K yy = 0.079 and added mass coefficient is M yy = 0.036.In this method, a function of the form p 1 + ω 2 p 2 is fit to the dynamic stiffness results using a nonlinear least squares method (SciPy [33] (v1.10) based on the algorithm of Moré [34].The curve-fit parameters p 1 and p 2 which result correspond to the static stiffness and the added mass coefficients, respectively.
Method 3: The added mass coefficient is estimated using the average difference between the CFD results and the CFD results subject to pure squeeze motion.This method is similar to Method 1, except that the pure squeeze results for the added mass coefficient are averaged over the perturbation frequency.This results in a single estimate for the added mass coefficient.Then, the added mass contribution estimated from this single added mass coefficient is subtracted from the dynamic stiffness.
Method 4: The asymptotic behavior of the frequency response function, H yy , in the high frequency limit is used to estimate the added mass coefficient [35].In this method, the added mass coefficient is estimated using the Bode magnitude plot of the H yy as depicted in Figure 4.The bode plot is consistent with the second-order linear system [36] model assumed for hydrodynamic fluid film, and the plot exhibits both low and high-frequency asymptotic behavior.The high frequency asymptote corresponding to Figure 4 is given by Equation ( 27), which is the equation for a linear line where the y-intercept is 31.18 and the slope is 37.15.
At the high frequency limit, the added mass coefficient dominates H yy as the added mass contribution is proportional to ω 2 , whereas the the damping and stiffness terms are multiplied by ω and one, respectively.The dominance of the added mass term are higher frequencies assumes that the dynamic coefficients themselves are frequency independent [16,37], or only mildly dependent on perturbation frequency.In the high frequency limit, the H yy may be approximated by Equation (28).
The added mass coefficients estimated using each of the four methods are summarized in Figure 5. Estimating the added mass coefficient using either a direct curve-fit of the dynamic stiffness (Method 2) or using a high-frequency asymptote (Method 4) provides nearly the same result.Making use of the pure squeeze CFD results (Methods 1 and 3) result in higher estimates of the added mass coefficient.As clearly indicated in Figure 5, only Method 1 retains the frequency dependence of the added mass coefficient.
If the added mass coefficient estimates in Figure 5 are used to remove the added mass components from Figure 3, the result is the static stiffness depicted in Figure 6.The negative static stiffness values predicted by Method 3 are non-physical for the linear slider bearing geometry.This indicates that the added mass coefficient is overestimated by averaging the pure squeeze CFD results over the range of perturbation frequencies considered.Figure 7 is a plot of static stiffness with respect to perturbation frequency with Method 3 removed.In this plot, Method 1 and Method 2 result in static stiffness estimates that fall within 20-30% of the analytic solution (perturbed Reynolds equation).If the added mass coefficient is considered constant and subtracted from the frequency dependent dynamic stiffness (as in Method 4), the static stiffness estimated from CFD is approximately 2× the value stemming from the perturbed Reynolds equation.It seems either Method 1 or Method 2 is most appropriate in estimating the dynamic coefficients.Methods 1 assumes that both K yy and M yy vary with perturbation frequency.Method 1, however, is computationally expensive as it requires both sliding and pure squeeze results be obtained at each perturbation frequency.Method 2 assumes that both K yy and M yy are both frequency independent and fits both parameters over the entire range of perturbation frequencies and seems to provide reasonable estimates for both static stiffness and added mass coefficients.By contrast, Methods 3 and 4 assume that M yy is constant (only), and K yy is either drastically under-or over-estimated at high perturbation frequencies.In the low frequency limit, Methods 1, 3, and 4 appear to converge to the analytic solution.
Figure 8 shows the dimensionless damping coefficient variation with perturbation frequency.The damping coefficient exhibits greater frequency dependence from 0 to 200 [rad•s −1 ] compared to the frequency range from 200 to 480 [rad•s −1 ], with the damping coefficient appearing to approach a constant asymptotic value in the higher frequency regime.The CFD-predicted damping coefficient is lower than the perturbed Reynolds equation result in the low frequency range and higher in the high frequency range with a cross-over perturbation frequency occurring somewhere between 60 and 80 [rad•s −1 ].As also evidenced in Figure 10 for the short journal bearing, the damping coefficients predicted using CFD and the perturbed Reynolds equation compare closest when the pertubation frequency coincides with the rotational frequency of the shaft.There appears to be some significance to the frequency at which the CFD-and Reynolds equation-predicted damping coefficients coincide and will be investigated in a future work.The coincidence appears to be related to the transition between sub-and super-synchronous excitation and/or the dominance of the stiffness or added mass terms in controlling the dynamic response of the hydrodynamic fluid film.Also, note that the cross-over frequency corresponding to the intersection of the low-and high-frequency asymptotes in Figure 4 8)) of linear slider bearing with respect to perturbation frequency.

Journal Bearing
Figures 9 and 10 are the dimensionless dynamic (added mass still present) and dimensionless damping coefficients, respectively, for short, L/D = 0.25, and long, L/D = 3.0, journal bearing geometries as a function of perturbation frequency ratio.The results were obtained in the absence of a cavitation model, i.e. solved subject to the full Sommerfeld condition.The perturbation frequencies were varied from 0.8× to 1.2× multiples of the rotational frequency Ω. Performing the CFD simulations near the rotational speed alleviates the need to resolve multiple time scales (high frequency perturbation) or run the simulations for long physical time with a small time step (low frequency perturbation) [38].At low perturbation frequencies, the simulation must be run for a long physical time to capture several oscillations of the perturbation boundary.However, the time step is constrained to a comparatively small value as dictated by the rotation of the shaft.For high frequency perturbations, the oscillation frequency reduces the required time step size below the value associated with resolving the rotation of the bulk fluid, i.e., the time step for the transient CFD simulation must be reduced.
In Figure 9, the cross-coupled stiffness coefficients predicted from CFD and the perturbed Reynolds equation compare well, and with the stiffness values exhibiting little variation with perturbation frequency.The negligible variation of the stiffness with perturbation frequency is consistent with the results of Ha and Yang [8].The authors [8] found little variation in both direct and cross-coupled stiffness coefficients for a load-on-pad tilting pad bearing (L/D = 0.5) subject to subsynchronous perturbation frequency ratios of the same order as the present work.
Theory based on the perturbed Reynolds equation under full Sommerfeld boundary conditions [25] predicts zero direct stiffness.The CFD predicted direct stiffness values are nonzero, but small compared to the cross-coupled stiffness values.The nonzero values likely stem from the finite nature of the harmonic perturbations applied to the journal within the CFD-based approach.Although small, 1/20 of the radial clearance of the bearing, the journal motion is finite and can only approach the limit of perturbed Reynolds equation theory in which the perturbations are truly infinitesimally small.0.80 0.85 0.90 0.95 1.00 1.05 1.10 1. 15   In Figure 10, the direct damping coefficients predicted from CFD and the perturbed Reynolds equation compare well.The invariance of the direct damping coefficients with perturbation frequency is consistent with results of Ha and Yang [8] who found the damping coefficients to be sensitive to perturbation frequency only at larger eccentricities.
The cross-coupled damping coefficients predicted from CFD are nearly zero, consistent with the perturbed Reynolds equation.The CFD-predicted direct damping coefficients compare well with the values from the perturbed Reynolds equation.The direct damping coefficients exhibit slight dependence on the perturbation frequency, but the variation is not significant in view of the numerical errors associated with the CFD results.If the dynamic stiffness values in Figure 9 are curve-fit over the perturbation frequency (Method 2 in the previous section), then Equations ( 29) and (30) for the the static stiffness and added mass coefficient matrices can be obtained.
The magnitudes of the dimensionless static stiffness essentially remain unchanged after the added mass contributions are removed.Unlike the extremely long slider bearing geometry, the role played by lubricant inertia (embodied in the added mass coefficients) in the short and long journal bearing is negligible for the perturbation frequencies and operating conditions considered.The magnitudes of the dimensionless added mass terms are misleading and are a direct result of the non-dimensionalization (Ω 2 scaling).The physical (dimensional) added mass terms are negligible in magnitude.An alternative non-dimiensionalization of the dynamic coefficients should be considered if the relative significance of the stiffness, damping, and added mass coefficients is to be assessed.
For both the short and long journal bearings considered, the CFD-predicted dynamic and static stiffness values are essentially the same, and along with the damping coefficients, compare well with coefficients predicted from the perturbed Reynolds equation.

Conclusions
This work presented a CFD-based frequency response method for the determination of the dynamic coefficients of hydrodynamic bearings.The method was derived from experimental parameter estimation techniques in which the dynamic coefficients of a bearing are determined through correlation of a known input displacement and a measured output force.The correlation takes the form of frequency response function (FRF) in which the the FRF is obtained numerically using a CFD model with a moving boundary (prescribed, harmonic journal displacement) rather than experimentally.
The CFD-based frequency response method is advantageous over the commonly applied perturbed Reynolds equation in that physics models of general applicability are more easily introduced.Examples include: fluid inertia, turbulence, void-transport based cavitation, thermal, and solid deformation.Of course, the CFD-based approach is more computationally expensive, but this concern has diminished as desktop computing power continues to increase and with more recent advent of cloud-based computing services [39].
The CFD-based frequency response method was verified by determining the dynamic coefficients of infinitely long linear slider and journal bearing geometries.The hydrodynamic fluid films and resulting bearing forces (static and dynamic) were resolved using CFD in the absence of cavitation, i.e., under full-Sommerfeld conditions.The authors are aware that cavitation is present in most bearings subject to ambient pressures at or near atmospheric pressure, and that the cavtiation is an essential, stabilizing ingredient within the bearing.However, there are circumstances in which cavitation is suppressed and the full-Sommerfeld solution has physical relevance.Bearings which operate at very high ambient pressures (pressurized system) do not experience as the sub-ambient pressures are not low enough to reach the vapor pressure of the lubricant (vaporous cavitation) or the saturation pressure of dissolved gases within the lubricant (gaseous cavitation).Additionally, if the surface dilatational viscosity of the lubricant is sufficiently large, it is possible for large subcavity pressures to persist with the cavitation zone, and the full-Sommerfeld solution is nearly reproduced [40].
The ever-presence of temporal fluid inertia within the CFD governing equations prompted the consideration of added mass coefficients within the mechanical models of the bearing hydrodynamic films.The added mass coefficients were shown to have a significant impact on the dynamic stiffness of the extremely long linear slider bearing geometry.The added mass coefficients were shown to be negligible for the short and long journal bearing geometries subject to perturbation frequencies ranging from 0.8× to 1.2× the journal rotational frequency and which correspond to approximately 1 Hz.Perturbation frequencies in this range (e.g., unbalance in a marine bearing) produce negligible added mass effects due to temporal inertia.
For the extremely long slider bearing considered, the added mass coefficient value was best-estimated either by removing the results obtained under pure-squeeze conditions (Method 1) or by curve-fitting the dynamic stiffness over a range of perturbation frequencies under the assumption of constant K and M (Method 2).Like the dynamic stiffness, the damping coefficient exhibited significant variation with perturbation frequency, a result the perturbed Reynolds equation cannot capture.
For the journal bearing geometries examined, the dynamic coefficients predicted using the CFD-based frequency method and those obtained by solving the perturbed Reynolds equation showed good coincidence.At perturbation frequencies near the rotational frequency, the dynamic coefficients were found to be nearly constant.The finding was consistent with previously published results [8].Moreover, the invariance of the dynamic coefficients with perturbation frequency appears to be a general finding owing to the disparity in the geometries and methods leading to this conclusion in the present work (plain journal bearing, numerical) and in [8] (titling pad bearing, experimental).
The results presented herein showed a discrepancy between the CFD-and Reynolds equation-predicted dynamic coefficients only for the extremely long slider bearing geometry.In a subsequent paper (Part 2), the theory developed herein (Part 1) will be applied to determine the dynamic coefficients of journal bearing geometries from CFD under more extreme operating conditions (higher rotating speeds, larger eccentricity).The CFD-predicted results will be directly compared with results obtained with the perturbed Reynolds equation with the concurrent effects of cavitation, turbulence, and advective inertia.

Figure 1 .
Figure 1.Linear slider bearing geometry and mechanical model of hydrodynamic fluid film.
2 ν ) which corresponds to a reduced sliding Reynolds number Re * u = 3.4 where Re * u = Re u h 2 b .Limiting the sliding speed of the bearing renders the flow laminar, eliminating the need, at this stage, for turbulence accounting.For a bearing with a diameter of 50.8 [mm], a sliding speed u = 1.88 corresponds to a rotational speed Ω = 74 [rad•s −1 ] or 706.8 [rpm].The rotational speed, while not characteristic in general for bearings, is in the bulk range of speeds developed within marine bearings.

Figure 2 .
Figure 2. Journal bearing (plain) geometry and mechanical model of hydrodynamic fluid film.D denotes the bearing journal diameter.Note that the coordinate system (O j , R, T) is NOT fixed at 45-deg rotation from the fixed reference frame (O b , X, Y).The origin of reference frame (O j , R, T) can vary depending on the static equilibrium position of the rotor ( , φ).

yyFigure 3 .
Figure 3. Dimensionless dynamic stiffness of linear slider bearing with respect to perturbation frequency.bω u

Figure 7 .Figure 8 .
Figure 7. Dimensionless static stiffness (Equation (7)) of linear slider bearing with respect to perturbation frequency.Method 3 removed and scaling of y-axis changed from Figure6.

Figure 9 .
Figure 9. Dimensionless dynamic stiffness coefficient of short and long journal bearing with respect to perturbation frequency ratio.The Reynolds number is fixed Re = ΩRC/ν = 323.9.RE refers to the dynamic coefficients obtained with the perturbed Reynolds equation.

Figure 10 .
Figure10.Dimensionless damping coefficient of short and long journal bearing with respect to perturbation frequency ratio.The Reynolds number is fixed Re = ΩRC/ν = 323.9.RE refers to the dynamic coefficients obtained with the perturbed Reynolds equation.

Table 4 .
Journal bearing OpenFOAM boundary conditions.Boundary conditions are required for the velocity U, the pressure P, and displacement D. "Front" and "Back" refer to the axial ends of the bearing.

Table 5 .
Dynamic coefficient variation with mesh density.Five periods and 200 time steps per period.

Table 6 .
Dynamic coefficient variation with number of time steps per period.Five periods with a 221 × 15 mesh.

Table 7 .
Dynamic coefficient variation with number of periods.Two hundred time steps per period with a 221 × 15 mesh.