ELSEVIER Pre ‐ Print Repository

This paper presents neural network regression models for predicting the static and dynamic reaction forces of spiral grooved gas journal bearings. The partial diﬀerential equations (PDEs) are sampled, based on a full factorial and randomly spaced parameter set. Feed-forward neural network (FNN) architectures are developed for modeling the PDEs and therefore replacing the time-consuming discrete and iterative solution procedure used to this date. A signiﬁcant speed-up factor of > 10 3 in computation time is achieved, compared to solving the PDE numerically. Furthermore, the FNN allows for multi-dimensional interpolation, which makes global system optimization easily possible. This is demonstrated by a real-case rotordynamic system optimization. By using the neural network meta-models, a complete rotordynamic system optimization time reduction of factor 300 is achieved.


Introduction
High-speed small scale turbomachinery is used in many different energy conversion systems such as domestic or commercial heat pumps [20], Organic Rankine cycles [18] or fuel cells [24]. The requirements of high rotational speeds and a high life-time expectation make gas lubricated bearings ideal for these systems. Gas bearings, such as the herringbone-groove journal bearing (HGJB), offer the advantages of oil-free operation, long lifetime, high rotational speed at relatively low frictional losses compared to other bearing types, no sealing requirements and avoidance of complex auxiliary systems. In order to develop stable compressor systems, accurate bearing models need to be available and parameter variations with thousands of computations have to be performed during an automated optimization process.
To this date, the performance of HGJB have been predicted by solving the thin-film flow equation, also called Reynoldsequation, in different ways. One of them is based on the narrow groove theory (NGT) approach, introduced by Whipple [23]. The NGT is developed by assuming an infinite number of grooves, which leads to a smooth-pressure distribution, which is governed by the NGT equation. The NGT equation intrinsically contains the bearing geometry. The NGT was further developed by Hirs [12], Malanoski and Pan [14], Vohr and Pan [21] and Vohr and Chow [22]. Fleming and Hamrock [4] used the narrow groove theory (NGT) for optimizing HGJB for maximum stability, quantified by a critical mass, which represents the mass the bearing can support before it gets unstable. A more detailed numerical analysis has been published by Bonneau and Absi [1], who solved the Reynolds-equation for HGJB based on a FEM approach. They presented solutions for different number of grooves and compared the case where the grooves are rotating and non-rotating. Iseli et al. [13] introduced the finite groove approach (FGA), which allows the computation of HGJB static and dynamic reaction forces with rotating grooves, avoiding a time-consuming transient analysis. They further presented a systematic comparison between the NGT and the FGA and limitations of the NGT were presented.
A detailed overview of the different modeling approaches is given by Gu et al. [6].
Nature of the issue. All of the above presented numerical methods involve time-intensive iterative procedures or lack the possibility of interpolating between different parameters. This is not a problem for the computation of one single bearing, but becomes problematic when optimization procedures require parameter variations in the order of thousands. The number of design variables for a rotordynamic system of small-scale turbo compressors is normally more than 10. For the design of such systems different analysis and optimization procedures are used, such as genetic algorithms or full factorial combinations and the number of computations easily exceeds millions. Therfore, one is interested in developing new methods, which are faster than already existing solvers.
Only a small number of prior work addresses this issue by using transfer functions of different types. Elrod et al. [3] introduced the step-jump response analysis, which was further developed by Miller and Green [15,16]. They analyzed the reaction motion of a self-acting cylindrical bearing by imposing an initial step-jump displacement. By using the analogy between viscoelastic elements and gas films, the bearing behavior is approximated by polynomial functions. A significant speed-up was reported, however, the usage of the derived models was limited to a simplified bearing model and to one specific geometry. Hassini and Arghir [9,10,11] used rational transfer functions for approximating the linearized bearing parameters in the frequency domain at different eccentric shaft positions. The linearized properties were computed by solving the non-linear fluid film equation. With the inverse Laplace transformation they obtained differential equations, describing the transient gas film behavior in a linear manner. They reported a speed-up in transient computation time of factor 2 [9]. Other transfer function methods, based on series and coefficient models, have been presented by Andrés and Jeung [19] for oil-bearings and by Franssen et al. [5] for aerostatic bearings. The bearing types used by Hassini and Arghir, Andrés and Jeung and Franssen et al. show moderate non-linear behavior compared to a HGJB. In order to model a HGJB, high order rational function (more than order 4) approximations are necessary. In addition, the strong non-linear behavior of HGJB in the geometrical and operational parameter space rules out the usage of simple linear, linear with interactions and linear-quadratic with interaction regression models such as normally used in a classical design of experiment analysis. A new approach is presented in this work, which is based on artificial neural network (ANN) regression models, capable of mapping several input parameters to a highly non-linear output solution space.
Goals and objectives. The goal is the derivation of artificial neural network regression models, describing the static and dynamic HGJB behavior within a given range of varying boundary conditions to decrease the computation time in comparison to existing numerical methods.
Scope of the paper. The PDEs, describing the static and dynamic behavior of HGJB, are numerically solved at full factorial and randomly spaced data points. The solution pressure fields are integrated in space in order to obtain the bearing reaction forces. Two-layer feed-forward neural network architectures are compared and selected based on a specified regression accuracy threshold and number of total neurons used. Speed-up and accuracy comparisons between the numerical analysis and the derived neural network are discussed and the possibility of multi-dimensional interpolation is presented by means of a rotordynamic case study.

Problem definition
Parameters. A gas bearing is characterized by multiple geometrical design parameters as well as by operating conditions. The geometrical parameters of a HGJB are schematically shown in Fig. 1 and listed in Table 1  parameters in total describe a HGJB, including operational conditions, such as rotational speed and eccentricity. The Table 1 Factors of herringbone-grooved journal bearing.

Factor
Symbol Unit Groove-width ratio -Groove angle °F ilm thickness ratio -Length to diameter ratio ∕ -Land to groove ratio -Eccentricity -Compressibility number Λ groove-ridge ratio is specified as follows: the film thickness ratio: the land to groove ratio: and the static eccentricity: The land to groove ratio is kept constant through the variations with a value of 1, which corresponds to a fully grooved bearing. The number of grooves is not varied, since it is assumed to be infinite based on the Narrow Groove Theory (NGT). Each parameter is given an upper and lower limit, which are listed in Table 2. The limits given in Table 2 are Table 2 Factor ranges for regression models. Bearing reaction forces. The reaction forces of the bearing are obtained from the fluid film pressure field, enclosed by the shaft and the bushing. The pressure inside a HGJB can be described by the NGT-equation, derived by Vohr and Chow [22]. Using the same notation for the NGT coefficients as Guenat and Schiffmann [7], the NGT equation in cylindrical coordinates for an ideal gas follows: where represents the circumferential coordinate, the axial coordinate, the non-dimensional pressure, the nondimensional time, the spiral-groove angle and Λ the compressibility number (see Fig. 1 for the coordinate system convention). The compressibility number Λ is a gas bearing specific non-dimensional number and defined as follows: where is the gas dynamic viscosity, the bearing radius, the ambient pressure, ℎ 0 the nominal bearing clearance in concentric position and Ω the rotational speed. For a given bearing geometry, the static reaction forces are only influenced by the compressibility number Λ. For fixed gas and geometrical properties, an increasing compressibility number can be understood as an increase in rotational speed. Therefore, the analysis of a gas bearing over a given speed range can be obtained in a non-dimensional way by varying the compressibility number.
The whirl ratio is the ratio between the bearing excitation frequency and rotational speed Ω. The coefficients and are specified by the bearing geometry and its operational condition. The coefficient's expressions are detailed by Guenat and Schiffmann [7].
Equation (5) represents a non-linear transient partial differential equation (PDE). In order to compute the static and dynamic HGJB properties, Eq. (5) is linearized around a given eccentricity and an excitation frequency dependent small orbital motion is imposed. This is done by introducing a perturbed clearance, varying with excitation frequency : where = cos( ) and = sin( ) are the perturbed gap height variations in circumferential direction. The radial displacements Δ , Δ ≪ 1 are the small amplitudes of the excited motion around the equilibrium shaft position.The equilibrium gap height is given as follows: = 0 , + cos( ) + sin( ). (8) where , are the static eccentricities in x-and y-direction. The pressure is also expressed as a perturbed property: Inserting Eqs. (7) and (9) where the subscript 0 marks the zeroth order properties. The dynamic pressure equations are obtained by collecting all first order terms ((Δ ), (Δ )). Two equations are obtained by separating all terms, which are multiplied by Δ and all terms multiplied by Δ . The equation for can be written as follows: The dynamic pressure equation for has exactly the same form, only the subscripts are replaced by . For solving Eqs. (10) and (11) the following Dirichlet boundary conditions are imposed at the bearing edges: On the whole fluid domain the pressure field has to fulfill the following periodicity requirement: The equations are discretized by a finite difference method (FDM) and the non-linearities are treated by the itertive Newton-Raphson method. The static bearing reaction forces are computed by integrating the unperturbed pressure field, which can be written in dimensionless form as follows: The bearing stiffness and damping values are obtained in analogy to the static forces by integrating the first order pressure fields and , obtaining the complex dynamic bearing force coefficients . Based on the spectral analysis by Pan [17], are computed for different whirl ratios = ∕Ω. In this work 41 whirl ratios are computed ranging from 0.01 to 2, including = 1, which corresponds to the synchronous whirl condition where = Ω. The 8 stiffness and damping values can then be extracted from as follows: The stiffness and damping values are finally used in the equations of motion for rotordynamic stability analysis.
In order to characterize a HGJB, the static forces and and the eight stiffness and damping values , for each whirl ratio need to be known. They are normally obtained by numerically solving one non-linear PDE for the static forces and two linear PDEs for the dynamic forces for 41 whirl ratios, yielding the computation of 83 PDEs in total. In the following, a new approach is presented, which model the PDEs by ANNs. The ANNs need the numerical solutions at certain sample points in the parametere space for training, but once trained, the bearing force coefficients can be obtained directly. The main advantages of using ANNs are that the complexity of computing the HGJB properties is significantly reduced (only analytical functions need to be evaluated, instead of solving complex PDEs) and the computations can be performed fully vectorized, which results in significant speed up factors.

Methodology
Feed-forward neural network models. The investigation of artificial neural network (ANN) architectures is done separately for the static and dynamic response variables. The base architecture is a feed-forward neural network (FNN), which represents a non-linear function of its inputs and is composed by its neuron functions [2]. It consists of a fullyconnected input layer, a few hidden layers and a fully connected output layer. A corresponding network is schematically shown in Fig. 2, which consists of 1 input layer with 6 inputs, two hidden layers with 7 and 4 neurons and 1 output layer with 1 output. A neuron is a non-linear mapping of several inputs to one output. Each neuron of a given layer is connected to each neuron of the preceding layer. A neuron is characterized by its activation function. On layer + 1  each neuron receives inputs from all neuron outputs at layer . These inputs are weighted and summed: where is the number of neuron input values, are the weights for each input and is an additional constant term, also named "bias". The value is termed the potential of the neuron. The output of the neuron is computed by a non-linear activation function and the potential . In this work, the hyperbolic tangent sigmoid activation function is selected, which is given as follows: This activation function is used for all intermediate (hidden) layers. The final predicted value is computed by the output neuron, which has a linear activation function: where is a constant and is computed by Eq. (17).
In this work, six parameters are used as inputs and one parameter as output, i.e. a separate neural network is trained for each bearing output (static reaction forces and stiffness and damping force coefficients). In order to train the network, the weights and the bias have to be adjusted such that the output predicts the response variable. The model prediction error can be quantified by the mean-squared error (RMSE), which has to be minimized during training. The RMSE is defined as follows: where are the number of samples, and are the predicted and sample output values. The RMSE is not ideal for comparing different output parameters, since it is not a normalized property. Therefore the normalized root meansquared error NRMSE is used, which is defined as: where and are the maximum and minimum values on the whole validation data set. The nominator of Eq. (21) represents the root mean squared error, which is normalized by the sample range.
The weights and bias of the FNN are updated by the Levenberg-Marquardt optimization algorithm by means of backpropagation, which is widely described in the literature (e.g. [2]). The Levenberg-Marquardt algorithm is an iterative numerical method, where the weights and the biases of the neural network are updated by the following expression: where is a scalar, is the vector containing all weights and biases, is the network error vector (computed by the loss function) and the corresponding Jacobian. Through iteration the weights and biases are updated, leading to a reduction in network errors. Backpropagation is one way of computing the Jacobian. The inputs and outputs of the neurons are computed by one forward-pass, the loss function is evaluated and its derivative computed. This value is then backpropagated to the first layer by computing the gradient of the weighted input of each layer. The two main advantages of this process are that there are no duplicated computations and that it computes the gradient directly to the output loss function and not to additional intermediate parameters. The net variation for the response variable is shown in Fig. 3a with the corresponding accuracy of fit in Fig. 3b.     ratios are selected in order to resolve the strong non-linear variations of dynamic force coefficients (see Fig. 12). The whirl ratio could be regarded as an additional input parameter to the models, reducing the number of models to 8.
However, the training data set would be larger by factor 41 and the model complexity would increase in order to handle   whirl ratio = 1, for the compressibility number interval Λ = [1,2], is shown in Fig. 9. It can be seen, that after 50 epochs the change of RMSE is not significant and a further training yields only a small increase in accuracy. The models are validated by 55k randomly distributed data points (geometry and eccentricity) in the compressibility range of Λ = [1,40], which have not been used for training. The accuracy of the models is quantified by the NRMSE. The results are shown in Figs. 10 and 11. In Fig. 10 the maximum NRMSE of all whirl ratios in a given compressibility number range, is plotted for all dynamic bearing coefficients. It can be seen that all NRMSE are below 0.01, which has been set as the training target. In Fig. 11 the maximum NRMSE of the whole validation data set is shown at different whirl ratios. It can be seen that the models perform equally well at different whirl ratios.
A comparison of predicted and excact solution for the impedances is shown in Fig. 12. The data point was not used in the training set. It can be seen that the model shows very good agreement with the NGT computed data. By using the meta-models a significant speed-up can be achieved compared to solving the NGT numerically ( Fig. 13). A similar trend as for the meta-models for the static forces can be observed. For one single computation, the meta-model is slower than using the NGT solver, which can be attributed to the neural network initialization. However, which is shown in the following section by means of a multi-objective rotordynamic system optimization.

Case study
The strength of the developed models is presented by means of a real case rotordyamic system optimization, encountered in the development phase of small scale oil-free turbomachines. In order to validate the ANN models, the  (14)). The rotordynamic stability is computed by a natural frequency analysis of the bearing-shaft system.
The pure radial motion of a shaft, modeled as a point-mass, can be described by the following system of equations: where is the mass matrix: the damping matrix: and the stiffness matrix: The following non-dimensional properties can be introduced: as the shaft displacement with as the bearing radius, as the dimensionless time, as the dimensionless mass, with the ambient pressure, ℎ 0 the bearing clearance, the dynamic fluid viscosity, Λ the compressibility number and the shaft mass. The non-dimensional damping and stiffness coefficients can be written as follows: With these definitions, Eq. (23) can be written as follows: By settinḡ =̄ Ω =̄ ̄ into Eq. (32), it follows the quadratic eigenvalue problem: which yields four eigenvalues, each coming as a complex conjugated pair of the following form: From these, the logarithmic decrement can be computed which is > 0 if the system is stable and < 0 if the system is unstable. It can be represented by the following expression: Two logarithmic decrements are computed for the two complex conjugated pairs of eigenvalues. One corresponds to a forward and one to a backward whirling shaft motion. Since the bearing impedance, and therefore the stiffness and damping matrices,̄ and̄ , are whirl ratio dependent, Γ has to be computed for each whirl ratio .
The minus signs are used since the optimization algorithm minimizes the scores and we are interested in maximum Γ , and values. The total number of individual computations per population is 120k. On an Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz processor, the computation time for one population (120k) is 147 s.

Results.
The Pareto curve of the rotordynamic system optimization at generation 10 is shown in Fig. 14 and suggests a clear trade-off between stability and load capacity. The corresponding design variables along the Pareto front are represented in Fig. 15, which suggest that increased load capacity is achieved by reducing the effect of the grooves (reduced groove width and depth) and by increasing the ∕ ratio. Figure 14: The Pareto curve for the optimized rotordynamic system in terms of minimum stability (at = 0 and 0.5) and maximum load capacity (for = 0.5) at generation 10. Gray region corresponds to unstable region.
A minimum logarithmic decrement of < 0 (gray region in Fig. 14) represents systems of unstable behavior and are not interesting in practice. The remaining geometries are chosen based on a trade-off between stability safety margin and achievable load capacity at eccentricity = 0.5.
Model performance and accuracy. The meta-model approach provides significant computation time improvements, compared to the direct finite difference method (FDM). For 10 generation a total number of 1.32 million single computations have to be performed. By using the meta-models, the computation of 10 generations, including the computation of the rotordynamic stability, takes a total time of 27 min. The same amount of computations, using the FDM with 4 parallel processes, would take (one single computation takes 1.5 s) 5.7 days, which corresponds to a factor ≈ 300, compared to the meta-model based optimization.
In order to assess the accuracy of the ANNs compared to the classical FDM approach in predicting the rotordynam-  Table 3. The relative error between FDM and meta-model approach for the minimum stability is 1.1 % and the error for predicting the static radial bearing reaction force is <0.1 %. (a) Groove-width ratio . (c) Groove depth ratio .

Summary and Conclusion
Data-driven meta-modeling of the static and dynamic analysis of HGJBs are presented in view of speeding up bearing computation time and enhancing bearing analysis possibilities.
Feed-forward neural networks (FNN) are used, offering many degrees of freedom and the possibility of fitting highly non-linear data sets. The models are tested and derived for predicting the static and dynamic HGJB behavior in different compressibility number intervals and eccentric shaft positions.
It can be concluded, that with two layer neural network architectures, a meta-model for the static reaction forces and can be trained with a NRMSE <7 × 10 −4 . For the dynamic force coefficients and a NRMSE <1 × 10 −2 is achieved.
The ANNs offer significant computational speed-up, in particular towards a high number of computations. At 10 4 computations the speed-up factor levels-off at > 10 5 for the static force models and > 10 3 for the dynamic force models. This means a significant reduction in bearing computation time, allowing for large (in the range of millions) parameter variations in minutes.
The ANNs are trained in parameter intervals, which gives the possibilities for multi-dimensional interpolation, avoiding discrete point limitations. The models can predict the bearing reaction forces at arbitrary eccentric positions, which makes them ideally suited for being used in a rotordynamic context. The subdivision into compressibility number intervals leads to a reduction of parameter space size, which makes it possible to use simple and fast FNN architectures.
Furthermore, it allows for easily expanding the models towards larger compressibility numbers.
The ANNs can be smoothly integrated in a rotordynamic context, allowing for time-efficient, optimizations. A relative difference of <1.1 % has been observed between the ANN based model and the classical FDM based approach for the target optimization parameter values. The meta-models lead to a speed-up factor for rotordynamic computations of ≈ 300, which reduces the computation time from several days to minutes.