An efficient three-dimensional foil structure model for bump-type gas foil bearings considering friction

This paper presents an efficient three-dimensional (3D) structural model for bump-type gas foil bearings (GFBs) developed by considering friction. The foil structures are modeled with a 3D shell finite element model. Using the bump foil mechanical characteristics, the Guyan reduction and component mode synthesis methods are adopted to improve computational efficiency while guaranteeing accurate static responses. A contact model that includes friction and separation behaviors is presented to model the interactions of the bump foil with the top foil and bearing sleeve. The proposed structural model was validated with published analytical and experimental results. The coupled elastohydrodynamics model of GFBs was established by integration of the proposed structural model with data on hydrodynamic films, and it was validated by comparisons with existing experimental results. The performance of a bearing with an angular misalignment was studied numerically, revealing that the reaction torques of the misaligned bearing predicted by GFB models with 2D and 3D foil structure models are quite different. The 3D foil structure model should be used to study GFB misalignment.


Introduction
Gas foil bearings (GFBs) have been applied to high speed and high performance rotating machinery in many industrial fields because of their low frictional losses and oil-free operation [1,2]. Designing hydrodynamic bearings requires the analysis of important characteristics such as the hydrodynamic pressure distribution [3] and the minimum film thickness [4]. The behavior of GFBs depends on both the gas film and compliant foil structure. The foil structure usually consists of a smooth top foil that acts as a bearing surface and a support structure that provides flexibility and frictional energy dissipation. Among the various types of GFBs, the bump-type foil bearing is most widely used, and it is illustrated in Fig. 1. The bump foil structure is one of the most critical aspects of GFB designs, as it can greatly improve the bearing performance [5]. However, the complex geometry and nonlinear friction behavior lead to challenges in analyzing its mechanical characteristics.
Walowit and Anno [6] first investigated the mechanics of an individual bump without considering the friction between the top foil and the bump foil. The stiffness formula for an individual bump was obtained analytically with a planar plate-bending model. Based on this formula, Heshmat et al. [7] developed a simple elastic foundation model (SEFM) by modeling the bump foil as individual springs while neglecting the friction and the interactions of bumps. The SEFM was coupled with a hydrodynamic film to analyze the performance of a GFB. Equivalent viscous damping was introduced in the SEFM to take into account the frictional dissipation of the foil structure in the dynamic analysis [8,9]. Iordanoff [10] proposed two formulas for calculating the stiffness of the welded and free individual bumps while considering friction forces. The SEFM is widely used in GFB studies because it is easy to couple with a hydrodynamic film model [11][12][13]. By adding the bump stiffness to the stiffness matrix of the top foil finite element (FE) model, the SEFM was integrated with the top foil for the analysis of GFBs [14,15]. The direct integration of the SEFM and the top foil shell model [15] resulted in a significant overestimation of the difference between the film thicknesses at the bearing midplane and the bearing edge because the analytical formulas employed were derived for planar models. Additionally, the analytical formulas used in the SEFM usually underestimate the bump foil stiffness. As the gas film and the foil structure operate in series, an accurate determination of the structural stiffness of the bump foil is critical for the calculation of bearing stiffness coefficients.
Advanced and comprehensive foil models were developed to obtain a more accurate structural stiffness of the bump foil by considering the actual geometry, the interactions between bumps, and the frictional forces. The first analytic model that considered both the friction and interactions of adjacent bumps was developed by Ku and Heshmat [16]. It showed that the friction and interactions of bumps significantly affect the behavior of the bump foil. Le Lez et al. [17] introduced an equivalent bump foil planar mechanical model. The bump foil is modeled as a network of interacting springs, and there are only three degrees of freedom (DOF) for each bump. This model is widely used as it considers the interactions of bumps and friction with low computational cost. Following the idea of Le Lez et al. [17], Feng and Kaneko [18] developed a simple bump model by replacing each bump with two rigid links and one horizontal spring. Hryniewicz et al. [19] established an analytical model for a bump foil based on Castigliano's second theorem, and Gad and Kaneko [20] presented a similar analytical model that allows flat segments between bumps to deflect laterally. The FE method was used to model the bump foil by using planar beam elements [21] and planar solid elements [22]. Hassan and Bonello [23] introduced a modal model of the bump foil structure while considering bump interactions and foil inertia. Hu and Feng [24] presented a beam model of the bump foil by considering the rounding radius of the bump to study the effect of rounding radius on the static performance of GFBs. Recently, Arghir and Benchekroun [25] further developed the foil model introduced in [17] by modeling the closeloose contacts between the top foil and the bump foil to consider possible gaps and manufacturing errors in the foil structure. This foil model was then coupled with the gas film model to determine all operating conditions of GFBs from start-up to full speed [26]. However, all of the above models are planar models, thus they cannot reflect the three-dimensional (3D) characteristics of a bump foil.
3D FE models have been used to consider the 3D effect of bump foils. Lee et al. [27] proposed a FE model for foil structures using a shell element considering the 3D shape of the bump foil. In this model, the transverse displacements of the bump top center are assumed to be equal to those of the top foil, and friction was not included. The numerical results showed that the deflections of bumps vary along the bearing width. Lehn et al. [28] presented a shell model for the top and bump foils of thrust foil bearings. The Coulomb friction was included in the penalty contact model for the top and bump foils. It was demonstrated that curved bump strips should be used in foil thrust bearings; i.e., the 3D detailed shape of a bump foil should be considered. Most other 3D foil structure models were built using commercial FE programs to model complex contact and friction behaviors [29][30][31]. The 3D FE model can be further used to study detailed wear behaviors [32]. The drawback of these 3D FE models of foil structure is that computations are much more time consuming than those for 2D FE models.
The 3D design for the foil structure has already led to better bearing performance. For example, the third generation GFB [5], whose bump stiffness changes in more than one direction, has a better bearing load carrying capacity. Another example is that the bump foil removed from the region of minimum film thickness improves the steady behavior of a GFB rotor system with a large unbalance [33]. Furthermore, GFBs frequently have angular misalignment due to rotor bending, unbalanced torque, or installation errors. Tolerance for bearing misalignment is an important issue for rotating machinery [34]. GFBs misalignment will lead to asymmetric hydrodynamic pressure distributions and foil deflections along the bearing axial direction [35]. In this case, the model of the foil structure should have the ability to represent the 3D effect accurately.
The aim of this study is to develop an efficient 3D foil structural model for bump-type GFBs-one that considers contact and friction behaviors. The top and bump foils are both modeled using 3D shell elements. The order of the bump foil model is reduced to improve the computational efficiency. The model order reduction is based on the fact that only a few nodes of the bump foil can be loaded. Hence, the unloaded nodes can be eliminated, and only the DOFs of possibly loaded nodes remain by the Guyan reduction method [36]. The component mode synthesis method is applied to maintain the sparsity of the stiffness matrix for the reduced model. Furthermore, a general contact model that includes friction and separation behaviors is used to model the interactions of the bump foil with the top foil and bearing sleeve.
The remainder of this paper is organized as follows. Section 2 introduces the detailed modeling of the foil structure, a description of the model order reduction for the bump foil, and the modeling of the bump foil contacts. Section 3 presents the coupling between the foil structure model and the gas film model, and then the coupled elastohydrodynamics model of GFBs is established. Section 4 provides the numerical results obtained to verify the accuracy and efficiency of the proposed model. The performance of the misaligned bearing is also analyzed. Finally, Section 5 concludes the article.

Modeling of the foil structure
The top foil and the bump foil are typical thinwalled structures. Therefore, the foil structures are modeled using the 3D shell elements so as to realize high computational efficiency. Without the loss of generality, the four-node mixed-interpolated shell element [37] shown in Fig. 2    The bump foil is also in contact with the bearing sleeve, and the corresponding contact forces are denoted as c,bs b ( ) f q . Additionally, the top foil is subjected to the hydrodynamic pressure p , with the corresponding generalized force vector denoted as . The equilibrium condition requires that the internal elastic forces reach a balance with the external forces, so the governing equations for the discrete foil structures are: where t K and b K denote the global stiffness matrices of the top foil and the bump foil, respectively. The contact forces, which include the normal contact force and tangential friction force, are the main sources of nonlinearity. Additionally, the friction forces significantly affect the vertical stiffness of the bump foil. Therefore, the modeling of frictional contacts is crucial and is presented in Section 2.2. The hydrodynamic pressure acting on the top foil element is given as: where a p is the ambient pressure, and n is the normal vector of the shell element.

Reduced-order model of the bump foil
Although efficient shell elements are used to model the foil structure, the FE model still has a large number of DOFs when one considers the actual 3D geometry of the bump foil. It is necessary to use the model order reduction technique to reduce computational costs and maintain small approximation errors. In fact, the bump foil is in contact with the top foil and bearing sleeve in only a few areas. This characteristic makes it possible to reduce the order and still guarantee accurate static responses with the Guyan reduction method [36]. In this study, it is assumed that the contacts of the bump foil only occur at its highest and corner nodes, as shown in Fig. 3(a). These potentially loaded nodes are treated as interface boundary nodes, and the leftover nodes are interior nodes. Based on the idea of component mode synthesis, the bump foil is divided into many substructures. These substructures are connected to adjacent substructures with shared interface boundary nodes.
Because only the boundary nodes are potentially loaded, the stiffness matrix of the substructure can be reduced to a small reduced stiffness matrix. Each substructure can be viewed as a super element. The stiffness matrix of the entire bump foil model can be obtained by assembling the stiffness matrices of substructures in a manner identical to that used with the standard FE assembling procedure. Different To maintain the sparsity of the global stiffness matrix to the greatest extent possible, the arc segment of a single bump is divided into two substructures and the flat segment is taken as a substructure, as shown in Fig. 3(a). For this substructure arrangement, the pattern of the global reduced stiffness matrix of the bump foil model is given in Fig. 3(b).
The theory of Guyan reduction is revisited briefly for clarity. The governing equation for a substructure is written as is the external force. By partitioning the DOFs of a substructure into boundary DOFs B q and interior DOFs I q , the governing equation of the substructure can be written in the following partitioned form: where the subscripts B and I refer to the boundary and interior, respectively. The second row of Eq. (3) is: In terms of coordinate transformations, the generalized coordinates vector q of a structure can be represented by boundary DOFs B q , as follows:

Modeling of frictional contacts
A general node-to-surface contact scheme is adopted here to model the contact behaviors of the bump foil. Each contact pair involves a slave node and a master surface. Taking a contact pair between the bump and top foils as an example, the highest nodes of the bump foil and elements of top foil are chosen as slave nodes and master surfaces, respectively. The schematics and nomenclature for the contact pair are illustrated in Fig. 4. The contact interaction between the bump foil and bearing sleeve can be handled in a similar way, in which the corner nodes of the bump foil are chosen as slave nodes, and the sleeve is the master surface. As the bearing sleeve is a fixed rigid plane, the procedure can be simplified and is not detailed here. All of these slave nodes are the boundary nodes of the reduced-order model of the bump foil. When a contact event is detected, the contact force The projection point P of the node Q follows the geometric conditions: where The penetration depth  of slave point Q within the master element is the key variable that indicates the contact state, and it is given by: where k is the normal contact stiffness, and e is the nonlinear exponent.
The Coulomb friction law is adopted to determine the tangential friction force which is characterized by a coefficient of friction f .  To avoid complex numerical calculations, a spring-type nonlinear force model [22] is used to simulate the Coulomb friction behaviors. The reaction force of the nonlinear spring should have a limited amplitude f n f  and oppose the displacement. The hyperbolic tangent function was adopted in this study, and the friction force is written as: where u is the tangential displacement, the shift s u is introduced to correct the sign of the friction force, and max u is the displacement tolerance in the sticking state. When sliding directions change, the shift s u is set to the current value of u . A suitable max u value should be chosen to achieve a good approximation of the Coulomb friction while avoiding numerical difficulty.

Modeling of the hydrodynamics film
The generation of pressure follows the Reynolds equation, which describes the flow of a thin film between two surfaces. For a compressible isothermal ideal fluid under steady state conditions, the Reynolds equation can be written in vector form as: x y with small angles, as shown in Fig. 5, the undeformed rigid film thickness is expressed as: where C is the radius clearance, and x  and y  are the tilt angles about the x-axis and y-axis, respectively.
The bearing is open to the environment on all sides, and hence the pressure at each edge is equal to the ambient pressure. For the single-pad gas foil bearing shown in Fig. 5, the boundary conditions for the Reynolds Eq. (11) are written as: where 0  is the leading-edge angle of the top foil. Such conditions allow for the generation of sub-ambient pressures. The full pressure condition is applied in the calculation of reaction forces and torques of the hydrodynamic film. The fluid film force, which greatly effects the behavior of a rotor system [38], is calculated as: The Reynolds Eq. (11) is discretized by the FE method. The approximate pressure field inside an element is interpolated as e p p   N p , where p N is the shape function matrix of the gas film element and e p is the element nodal pressure vector. After applying the Bubnov-Galerkin method [39,40], the weak form of the Reynolds Eq. (11) for an element e is given as: N is the gradient matrix of shape functions. The bilinear four-node elements are used to mesh the gas film. The governing equations of the gas film are obtained as a system of nonlinear algebraic equations by assembling all of the element equations.

Coupling of the gas film and foil structure
The interactions between the gas film and foil structure are described as follows: the film pressure acts upon the top foil and causes the foil structure to deflect, and the foil deflection changes the gas film height distribution and affects the hydrodynamic pressure. Fig. 6(a) shows the schematic diagram of the numerical model of a GFB with associated nomenclature. The generalized hydrodynamic force vector that acts on a shell element of the top foil is calculated with Eq. (2). The elastic deflection that contributed to the film thickness is taken as the opposite vertical displacement of the top foil; i.e., f h v   , as shown in Fig. 6(a).
The flowchart of the solution for the elastohydrodynamic problem is illustrated in Fig. 6(b). The discrete Reynolds equation is solved by the NR method with the given foil deflections, and the foil deflection is then calculated with the updated hydrodynamic pressure. The iteration is repeated until both the hydrodynamic pressure and foil deflections converge. The eccentricities must also be updated until the hydrodynamic reaction force and the static load of the rotor achieve balance.

Validation of the structural model
Two different bump strips were modeled and simulated to verify the proposed structural model. The parameters of two bump strips are listed in Table 1. For each model, the number of elements in the axial direction is 50; the flat and arc segments of a single bump are meshed into 3 elements and 12 elements in the circumferential direction, respectively. The mesh independence had been checked by comparison with the results for double refined meshes.  Case 1 is a typical bump strip with ten bumps, and it was applied with two different distributed loads. One involved a uniform distribution with a pressure of 0.2 MPa, and the other involved an increasingdecreasing distribution with a maximum pressure of 0.4 MPa. The calculated results for the present reducedorder model were validated by comparison with the analytical model and FE analysis in Ref. [17], as shown in Fig. 7. For the present 3D model, the bump deflection was obtained from the mean deflections of all of the top nodes along the axial direction of a single bump. In addition, computational efficiencies of the reduced-order and the full-order models were compared. The number of DOFs for the orderreduced model and the full-order model are 7,905 and 38,505, respectively. The calculation with the orderreduced model required approximately one-seventh of the time needed for the full-order model (the two computational programs had both been optimized to the maximum extent possible).
In Case 2, the bump strips [22] tested under the influence of quasi-static loading and unloading were simulated to validate further the structure and friction model. In the simulations, the load was increased to the maximum value in 50 steps and then decreased to zero in 50 steps. The hysteresis force-displacement curves for bump strips with four bumps and six bumps, calculated with the present model, are plotted in Fig. 8. The results of FE models built by the commercial software ABAQUS, and the test results were used for comparison. As can be seen, the present results are in good agreement with the other two results. The frictional energy dissipation of the bump structure can be simulated correctly with the proposed model.

Validation of the coupled hydrodynamics model
The coupling of the hydrodynamics film and foil structure was verified by examining the static performance calculations, and this will be discussed in this subsection. A typical first generation GFB, for which major static performances had been tested experimentally by Ruscitto [41] in 1978, was simulated here. The parameters of the bump foil of the test bearing are listed in Table 1, Case 1. The thickness of the top foil is 0.102 mm. The bearing radius is 19.05 mm, and the radial clearance C is 31.8 μm. The ambient pressure is 5 1.0135 10 Pa  , and the gas viscosity is 5 1.95 10   Pa·s. The gas film and top foil are each divided into a mesh of 104 and 30 elements in the circumferential and axial directions, respectively. The present model did not include the stiffening factor introduced in Ref. [15] for the top foil model.
The film thickness distribution at the midplane of the bearing was calculated under a bearing load of 134.1 N at a speed of 30,000 rpm. The predicted film thickness curve is plotted in Fig. 9(a) and compared with the test result. As can be seen, the predicted result is in good agreement with test data in the region of minimum film thickness. Noteworthy sags appear on the top foil between the adjacent two bumps in this region. For the regions of large film thickness, the variations in the thickness distributions are consistent with the test result. In the beginning and end portions of the gas film, sub-ambient pressures cause the top foil to separate from the bump foil and deflect inward to the rotor, so that the film thickness decreases. The welded condition prevents the top foil from deflecting inward; hence, a wave of pressure and film thickness occur near the welded edge, as shown in Fig. 9(b) and 9(c). This phenomenon was also reported in prior work [14,33].
In addition, Fig. 10 shows the minimum film thicknesses located at the midplane and edge of the bearing under various loads at a speed of 45,000 rpm. The bearing clearance is assumed to be 21 m to achieve better predicted results, as proposed in Ref. [42]. As shown in Fig. 10, the predicted results match quite well with the test results overall, suggesting the correctness of the present model.

Performance of the misaligned bearing
The characteristics of a GFB experiencing angular misalignment are described in this subsection. The GFB parameters in Section 4.2 are used here. In this study, the rotation center of the shaft is chosen as the equilibrium position of the shaft under a given load and with perfect alignment. With asymmetrical distributions of hydrodynamic pressure, the reaction torques of the bearing are calculated as follows: Two coupled GFB models with 2D and 3D foil structures were compared in the analysis of shaft misalignment. The 2D foil structure is modeled using the planar Timoshenko beam elements [43], with the assumption that variations in deflection in the axial    has the largest value at the same tilt angle. As expected, the angular stiffness (slope of the torque-angle curve) of the two models is quite different. The angular stiffness calculated by the GFB model with the 2D foil structure model is greater. This is because the gas film and the foil structure work in series, and the axial rigidity of the 2D foil structure is assumed to be infinite whereas that of the 3D foil structure is finite. In terms of the pressure distribution, the peak pressure of the GFB model with the 2D foil structure is larger at the same tilt angle, as shown in Fig. 12. As the shaft tilts about the y-axis, the minimum film thickness decreases, and the peak pressure increases and moves toward the edge of the bearing. This leads to more deformations in the 3D foil structure model near the edge of the bearing. However, the 2D foil structure model cannot reflect the variations in pressure in the axial direction because the load of the gas film acting on the 2D foil structure is the mean pressure across the bearing width. Thus, the minimum film thickness of the 2D foil structure model is smaller than that of the 3D   Moreover, the signs of x M calculated by the GFB models with 2D and 3D foil structures are opposite for the shaft tilting about the y-axis, as shown in Fig. 11(b). This is caused by the two different sub-ambient pressure distributions. When omitting the sub-ambient pressure in the calculation of torques, the torque x M calculated by the GFB model with the 2D foil structure becomes negative. As depicted in Fig. 12, the sub-ambient pressure distribution of the GFB model with a 3D foil structure model is almost symmetric about midplane of the bearing ( whereas for the GFB model with the 2D foil structure model, it is not. The pressure distributions are related to the film thickness distributions. The undeformed rigid heights of the two GFB models are the same, and both are antisymmetric about the midplane. For the 2D structure model, the deflections in the axial direction are constant; the film thickness is antisymmetric about the midplane, as shown in Fig. 13(a), resulting in an asymmetric pressure distribution. For the 3D structure model, the top foil locally separates from the bump foil in areas where sub-ambient pressure develops. As shown in Fig. 13(b), the film thicknesses in the regions 0-90 and 270-360 are relatively symmetrical about the midplane, resulting in a fairly symmetrical sub-ambient pressure distribution.

Conclusions
The 3D design of the compliant foil structure of GFBs has been employed recently to obtain better bearing performance. This paper presents an efficient 3D structural model for bump-type GFBs.
The foil structures are modeled as a 3D FE shell model. Using the bump foil mechanical characteristic, the Guyan reduction and component mode synthesis methods are adopted to improve the computational efficiency while guaranteeing accurate static responses. A general contact model that includes friction and separation behaviors is presented to model the interactions of the bump foil with the top foil and bearing sleeve. The proposed reduced-order model of the bump structure was verified by analytical and experimental results. In terms of efficiency, calculations with the reduced-order model are seven times faster than those with the full-order model. The coupled elastohydrodynamics model of GFBs was established by the integration of the proposed structural model with the hydrodynamic film. The coupled GFB model was validated by comparisons with existing experimental results. The performance of a bearing exhibiting angular misalignment was studied numerically. The numerical results show that there are clear differences between the reaction torques obtained with 2D and 3D foil structure models. The axial rigidity of the foil structure affects the angular stiffness of the bearing. In addition, the top foil locally separates from the bump foil along the axial direction in the sub-ambient pressure region, resulting in a fairly symmetrical sub-ambient pressure. Therefore, the 3D foil structure model should be used for the analysis of bearing misalignment.
Other types of support structures for GFBs, e.g., compression springs [44], also exhibit the characteristic of fewer loaded nodes than unloaded nodes. Hence, the proposed method of model order reduction can be also applied to consider the 3D effects of the foil structure with relatively low computational cost. Furthermore, although this work is focused on static analysis, the proposed foil model can be further extended to dynamic simulations of GFBs performed while taking fixed-interface normal modes into account.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/. Yongpeng GU. He received his B.S. degree in 2012 and his Ph.D. degree in 2017, both from Tsinghua University, China. In 2017, he joined the School of Aerospace Engineering in Tsinghua University as a postdoctoral researcher. His current research interests include gas foil bearings, dynamics of rotor-bearing systems, and multibody dynamics.