The Fluid-Structure-Thermal Performance Analysis of Gas Foil Thrust Bearing by Using Computational Fluid Dynamics

: This paper presents a detailed three-dimensional (3D) thermo-elastic-hydrodynamic (TEHD) multi-physics model of the bump-type gas foil thrust bearing based on computational ﬂuid dynamics (CFD). In this model, the moving mesh technology is applied in the ﬂuid ﬂow domain, and the new boundary condition of fully developed ﬂow is applied at the inlet and outlet boundaries, which is consistent with the continuous property of ﬂuid ﬂow and has better convergence performance in CFD. The groove between pads is set as the symmetry boundary. The contact pairs with Coulomb friction and contact/separation behaviors are considered in the structure deformation and heat transfer. The simulation results indicated that the boundary pressure has a signiﬁcant inﬂuence on the foil deformation. It also revealed the heat ﬂux transfer path and temperature distribution in the gas foil thrust-bearing (GFTB) system.


Introduction
GFTB is an environmentally friendly hydrodynamic bearing that withstands axial load, which is widely used in high-speed turbomachinery due to its advantages of higher reliability, lower operating costs, and soft failure, etc. [1]. As the requirements of rotational speed increase and the operating environment became complex, the traditional calculation methods with assumptions of constant temperature cannot accurately predict the bearing performance [2].
Theoretical and experimental research considering the influence of temperature has been conducted. Scholars have established different forms of temperature field calculation models of gas foil bearings. Salehi et al. [3] simplified the energy equation by using the Couette approximation method to handle the relationship between air viscosity and temperature, but ignored the influence of the axial pressure gradient with only circumferential temperature distribution being calculated. Peng et al. [4] simplified the foil structure as a spring model to couple Reynolds equations with heat transfer equations. They considered the forced convection outside the top foil, but ignored the heat transfer to the rotor and bearing sleeve. Paulsen et al. [5] compared the differences in the results of dynamic characteristics based on isothermal and thermal models of three different types of journal foil bearing. However, the axial heat flux in the shaft was not considered, and the calculated temperature was a bit too high. Lee et al. [6] presented the heat transfer models of top foil, bump foil, cooling channel, and thrust plate in detail under the symmetric arrangement of two GFTBs and considered the effect of temperature on gas viscosity in Reynolds' equation. However, the bump foil was simplified as a lumped thermal mass model. Hu et al. [7] established a nonlinear bump stiffness model considering the friction force and obtained the performance of GFTBs by coupling the foil deformation equation with Reynolds and energy equations. Those models mentioned above calculated the flow field by solving the

Methodology
In this study, the three-dimensional TEHD computational model is established by using COMSOL software, which consists of three coupled physical fields: fluid flow, structural deformation, and heat transfer. In this section, the main problems involved in each single physical field and coupled physical fields are described, respectively.

Geometrical Model and Computational Domain of GFTB
The detailed structure of bump-type thrust bearing with six pads is shown in Figure 1, which mainly consists of a smooth surface top foil with a slight slope segment and a plane segment to provide gas lubricating conditions. The bump foil, shaped of corrugated and flat parts that are alternately connected, is located under the plane section of top foil to provide elastic support. The spacer shim is used to adjust the slope height δ h , and a backing plate is used to fix the three components above.
Lubricants 2022, 10, x FOR PEER REVIEW 3 of 20 paper is organized as follows. First, all the governing equations and principles involved in each physical field were introduced. Secondly, the model reliability was verified by comparing it with the literature results and the numerical simulations. Finally, the thermal transfer, structural deformation, and hydrodynamic performance of gas foil thrust bearings under different conditions were analyzed.

Methodology
In this study, the three-dimensional TEHD computational model is established by using COMSOL software, which consists of three coupled physical fields: fluid flow, structural deformation, and heat transfer. In this section, the main problems involved in each single physical field and coupled physical fields are described, respectively.

Geometrical Model and Computational Domain of GFTB
The detailed structure of bump-type thrust bearing with six pads is shown in Figure  1, which mainly consists of a smooth surface top foil with a slight slope segment and a plane segment to provide gas lubricating conditions. The bump foil, shaped of corrugated and flat parts that are alternately connected, is located under the plane section of top foil to provide elastic support. The spacer shim is used to adjust the slope height δh, and a backing plate is used to fix the three components above. Due to the assumption of neglecting misalignment of thrust runner, the computational domains of the entire bearing structure and fluid including six integrated pads could be simplified to one considering the rotational symmetry boundary conditions of sectoral gas film. The investigation of thermal characteristics is also conducted based on the above simplified one-pad mode, which could greatly reduce the computational workload. The final model of the physics domains is presented in Figure 2. Due to the assumption of neglecting misalignment of thrust runner, the computational domains of the entire bearing structure and fluid including six integrated pads could be simplified to one considering the rotational symmetry boundary conditions of sectoral gas film. The investigation of thermal characteristics is also conducted based on the above simplified one-pad mode, which could greatly reduce the computational workload. The final model of the physics domains is presented in Figure 2.
Lubricants 2022, 10, x FOR PEER REVIEW 3 of 20 paper is organized as follows. First, all the governing equations and principles involved in each physical field were introduced. Secondly, the model reliability was verified by comparing it with the literature results and the numerical simulations. Finally, the thermal transfer, structural deformation, and hydrodynamic performance of gas foil thrust bearings under different conditions were analyzed.

Methodology
In this study, the three-dimensional TEHD computational model is established by using COMSOL software, which consists of three coupled physical fields: fluid flow, structural deformation, and heat transfer. In this section, the main problems involved in each single physical field and coupled physical fields are described, respectively.

Geometrical Model and Computational Domain of GFTB
The detailed structure of bump-type thrust bearing with six pads is shown in Figure  1, which mainly consists of a smooth surface top foil with a slight slope segment and a plane segment to provide gas lubricating conditions. The bump foil, shaped of corrugated and flat parts that are alternately connected, is located under the plane section of top foil to provide elastic support. The spacer shim is used to adjust the slope height δh, and a backing plate is used to fix the three components above. Due to the assumption of neglecting misalignment of thrust runner, the computational domains of the entire bearing structure and fluid including six integrated pads could be simplified to one considering the rotational symmetry boundary conditions of sectoral gas film. The investigation of thermal characteristics is also conducted based on the above simplified one-pad mode, which could greatly reduce the computational workload. The final model of the physics domains is presented in Figure 2.  The pressure distribution and the comparison of load capacity for the entire CFD model and symmetry CFD model are illustrated in Figure 3. It is concluded that the pressure distribution is the same for those two models and the magnitude of load capacity is also relatively similar where the errors E 0 are less than 1.6 × 10 −4 in the situations of this model. E 0 = (F 2 -F 1 )/F 1 , where F 1 (N) is the load capacity of entire model, and F 2 (N) is the load capacity of simplified symmetry model. Thus, the simplified symmetry model can be applied to subsequent analysis. The pressure distribution and the comparison of load capacity for the entire CFD model and symmetry CFD model are illustrated in Figure 3. It is concluded that the pressure distribution is the same for those two models and the magnitude of load capacity is also relatively similar where the errors E0 are less than 1.6 × 10 −4 in the situations of this model. E0 = (F2 -F1)/F1, where F1 (N) is the load capacity of entire model, and F2 (N) is the load capacity of simplified symmetry model. Thus, the simplified symmetry model can be applied to subsequent analysis.

Modeling of Fluid Flow
In the fluid component of this model, the gas film is assumed to be compressible fluid flow, the full Naiver-Stokes equations are solved by CFD software utilizing a finite element approach, and relevant fluid-governing equations are presented as follows. The continuity equation is formulated as where ρ is the density of fluid (kg/m 3 ) and u is the velocity vector (m/s). The momentum equation is formulated as

Modeling of Fluid Flow
In the fluid component of this model, the gas film is assumed to be compressible fluid flow, the full Naiver-Stokes equations are solved by CFD software utilizing a finite element approach, and relevant fluid-governing equations are presented as follows. The continuity equation is formulated as where ρ is the density of fluid (kg/m 3 ) and u is the velocity vector (m/s). The momentum equation is formulated as where τ is the viscous stress tensor (Pa) and F is the volume force vector (N/m 3 ). The energy equation of non-isothermal flow is formulated as where p is pressure (Pa), C p is the specific heat capacity at constant pressure (J/(kg·K)), T is the absolute temperature (K), q is the heat flux vector (W/m 2 ), Q contains the heat sources (W/m 3 ), and S is the strain-rate tensor.
Combining the structure feature of gas thrust bearing and the principle of dynamic pressure gas lubrication, the main gas film thickness distribution is shown in Figure 4, and its calculation equation is as follows, where h 2 is the minimum thickness of initial gas film, which locates in the plane section (m), g(θ) is the distribution of initial gas film thickness which subtracts the value of h 2 (m), β is the angle of a top foil ( • ), b is the pitch ratio, and d(r,θ) is the deformation of top foil under aerodynamic force. where is the viscous stress tensor (Pa) and F is the volume force vector (N/m 3 ). The energy equation of non-isothermal flow is formulated as where p is pressure (Pa), Cp is the specific heat capacity at constant pressure (J/(kg·K)), T is the absolute temperature (K), q is the heat flux vector (W/m 2 ), Q contains the heat sources (W/m 3 ), and S is the strain-rate tensor.
Combining the structure feature of gas thrust bearing and the principle of dynamic pressure gas lubrication, the main gas film thickness distribution is shown in Figure 4, and its calculation equation is as follows, where h2 is the minimum thickness of initial gas film, which locates in the plane section (m), g (θ) is the distribution of initial gas film thickness which subtracts the value of h2 (m), β is the angle of a top foil (°), b is the pitch ratio, and d (r,θ) is the deformation of top foil under aerodynamic force. Otherwise, modeling of the connecting groove between pads is also considered in the present work, which is shown in Figure 2b. The groove is divided into exactly the same two parts, which are connected in the flat section and the slope section of the gas film, respectively. The reason for this modeling form is that the two rotational periodic symmetry boundaries require the geometry shape and dimension to be exactly the same, to enable the computational data to be transferred completely.
Due to the fact that the magnitude of load-carrying capacity is directly affected by the minimum film thickness presented in the computational domain model, and the fact that the change in the thickness needs to deal with the geometry model again, it is difficult to parameterize the micron-scale gas film thickness, and the geometry model needs to be cleaned up and patched in the preprocessing step. Thus, the film thickness is set constant to analyze the performance of the GFTBs conveniently, and the loadcarrying capacity discussed here is also not the maximum value of foil structures. Otherwise, modeling of the connecting groove between pads is also considered in the present work, which is shown in Figure 2b. The groove is divided into exactly the same two parts, which are connected in the flat section and the slope section of the gas film, respectively. The reason for this modeling form is that the two rotational periodic symmetry boundaries require the geometry shape and dimension to be exactly the same, to enable the computational data to be transferred completely.
Due to the fact that the magnitude of load-carrying capacity is directly affected by the minimum film thickness presented in the computational domain model, and the fact that the change in the thickness needs to deal with the geometry model again, it is difficult to parameterize the micron-scale gas film thickness, and the geometry model needs to be cleaned up and patched in the preprocessing step. Thus, the film thickness is set constant to analyze the performance of the GFTBs conveniently, and the loadcarrying capacity discussed here is also not the maximum value of foil structures.

Modeling of Bearing Structure
The computational domain is shown in Figure 2; the backing plate and spacer shim material usually use structural steel, which is assumed for a fixed rigid body domain, and the foil material is a nickel-chromium alloy which has great elastic deformation properties. Thus, the main deformation components of GFTBs are the top foil and bump foil under high-speed gas pressure. The deformation of foils follows Hooke's law written as where σ is the second-order stress tensor, and ε is the second-order strain tensor, which is related to the deformation d t (r, θ) and C is the fourth-order elasticity tensor. The equilibrium equation of steady state is given by Newton's second law as follows, where f V is the body force per unit deformed volume, and ∇ x is the gradient operators taken with respect to the spatial coordinates. Meanwhile, the contact between bump foil and backing plate, as well as between the bump foil and the top foil, has a significant influence on the foil deformation too. The penalty function method is used in the multi-physics coupled contact, which is equivalent to adding a spring boundary load between the surfaces to complete the virtual contact. Furthermore, to implement loading gradually, a loading function can be defined here. In this contact model, in addition to considering the penalty type normal force, a Coulomb friction law is also added to refine the model, with a friction coefficient µ = 0.1 and a friction penalty factor f t = 0.1. The selection of the friction coefficient value is based on the values used frequently in similar GFTBs simulations or experiments analysis from other literatures, e.g., Refs. [25,26].
The contact pairs in this structure domain are shown in Figure 5. On the settings of the contact pairs, the easily deformable surface and the surface with a small bending radius are set as the target surfaces, and the surface that is difficult to deform is set as the source surface. Furthermore, the rigid body must be the source surface. The mesh size of the target surface needs to be smaller. In general, it is half the size of the source face mesh. In addition, the contact penalty factor f p could control the convergence speed, and stability, and increasing the penalty factor can speed up the solution process, but it is also easy to cause the solution to diverge. The range of the factor is set to (0.05, 0.2) in this structure model.

Modeling of Bearing Structure
The computational domain is shown in Figure 2; the backing plate and spacer shim material usually use structural steel, which is assumed for a fixed rigid body domain, and the foil material is a nickel-chromium alloy which has great elastic deformation properties. Thus, the main deformation components of GFTBs are the top foil and bump foil under high-speed gas pressure. The deformation of foils follows Hooke's law written as where σ is the second-order stress tensor, and ε is the second-order strain tensor, which is related to the deformation d r θ and C is the fourth-order elasticity tensor. The equilibrium equation of steady state is given by Newton's second law as follows, where fV is the body force per unit deformed volume, and x ∇ is the gradient operators taken with respect to the spatial coordinates. Meanwhile, the contact between bump foil and backing plate, as well as between the bump foil and the top foil, has a significant influence on the foil deformation too. The penalty function method is used in the multi-physics coupled contact, which is equivalent to adding a spring boundary load between the surfaces to complete the virtual contact. Furthermore, to implement loading gradually, a loading function can be defined here. In this contact model, in addition to considering the penalty type normal force, a Coulomb friction law is also added to refine the model, with a friction coefficient μ = 0.1 and a friction penalty factor ft = 0.1. The selection of the friction coefficient value is based on the values used frequently in similar GFTBs simulations or experiments analysis from other literatures, e.g., Refs. [25,26].
The contact pairs in this structure domain are shown in Figure 5. On the settings of the contact pairs, the easily deformable surface and the surface with a small bending radius are set as the target surfaces, and the surface that is difficult to deform is set as the source surface. Furthermore, the rigid body must be the source surface. The mesh size of the target surface needs to be smaller. In general, it is half the size of the source face mesh. In addition, the contact penalty factor fp could control the convergence speed, and stability, and increasing the penalty factor can speed up the solution process, but it is also easy to cause the solution to diverge. The range of the factor is set to (0.05, 0.2) in this structure model.

Modeling of Heat Transfer
The only heat source in this computational model comes from the shear stress of the fluid film, which is solved with the three-dimensional energy equation. After the heat flux being generated in the gas film, it flows to the rotor disk side and the thrust-bearing side through heat in conduction. On the rotor disk side, the heat flux spreads out through natural convection. On the thrust-bearing side, the heat flux from the top foil is conducted to

Modeling of Heat Transfer
The only heat source in this computational model comes from the shear stress of the fluid film, which is solved with the three-dimensional energy equation. After the heat flux being generated in the gas film, it flows to the rotor disk side and the thrust-bearing side through heat in conduction. On the rotor disk side, the heat flux spreads out through natural convection. On the thrust-bearing side, the heat flux from the top foil is conducted to the bump foil and then to the backing plate through thermal contact. Meanwhile, the heat flux is also conducted to the base plate through the air between foils, and finally it spreads out from backing plate surface to environment through natural convection. The heat flux transfer schematic is shown in Figure 6. The heat conduction and heat convection equations are as follows, where q is the heat flux by conduction (W/m 2 ), k is the thermal conductivity (W/(m·K)), ∇T is the temperature gradient (K), ρ is the density (kg/m 3 ), C p is the specific heat capacity at constant stress (J/(kg·K)), T is the absolute temperature (K), and Q contains heat sources other than viscous dissipation (W/m 3 ).
Lubricants 2022, 10, x FOR PEER REVIEW 7 of 20 the bump foil and then to the backing plate through thermal contact. Meanwhile, the heat flux is also conducted to the base plate through the air between foils, and finally it spreads out from backing plate surface to environment through natural convection. The heat flux transfer schematic is shown in Figure 6. The heat conduction and heat convection equations are as follows, where q is the heat flux by conduction (W/m 2 ), k is the thermal conductivity (W/(m·K)), T ∇ is the temperature gradient (K), ρ is the density (kg/m 3 ), Cp is the specific heat capacity at constant stress (J/(kg·K)), T is the absolute temperature (K), and Q contains heat sources other than viscous dissipation (W/m 3 ). Thermal contact is an important heat transfer route in this TEHD model, in which contact conductivity hc (W/(m 2 ·K)) is mainly influenced by the surface roughness and contact pressure coupled with structural mechanics generally, and thermal conductivity of air gap ℎ (W/(m 2 ·K)) depends on the medium properties of the gap space. The region of thermal contact pairs is the same as the region of the structure contact pairs mentioned in the last section, as shown in Figure 4. The inlet and outlet surfaces of the fluid flow domain are always considered as pressure inlet and pressure outlet, and the pressure is assumed to be constant at the same value p0 (gauge pressure). In most of the numerical simulations p0 = 0 (Pa), there are other scholars, like Gen Fu et al. [16], who assumed p0 to be slightly higher than ambient pressure. The pressure inlet and pressure outlet do not suppress reverse flows of the boundaries, which allows the simulation results of the gas film flow to obtain a more realistic flow trend.
In the present study, the fully developed flow boundary condition is applied to the inlet and outlet of channel, and this type of boundary condition has been applied to the Thermal contact is an important heat transfer route in this TEHD model, in which contact conductivity h c (W/(m 2 ·K)) is mainly influenced by the surface roughness and contact pressure coupled with structural mechanics generally, and thermal conductivity of air gap h g (W/(m 2 ·K)) depends on the medium properties of the gap space. The region of thermal contact pairs is the same as the region of the structure contact pairs mentioned in the last section, as shown in  The inlet and outlet surfaces of the fluid flow domain are always considered as pressure inlet and pressure outlet, and the pressure is assumed to be constant at the same value p 0 (gauge pressure). In most of the numerical simulations p 0 = 0 (Pa), there are other scholars, like Gen Fu et al. [16], who assumed p 0 to be slightly higher than ambient pressure. The pressure inlet and pressure outlet do not suppress reverse flows of the boundaries, which allows the simulation results of the gas film flow to obtain a more realistic flow trend.
with a porous layer, wherein the velocities and turbulent quantities were calculated separately. Batista [29] applied the fully developed water inlet boundary condition in a crossflow air-to-water fin-and-tube heat exchanger, and compared it with uniform flow velocity and temperature profiles; its numerical results coincide well with experimental measurements. Compared to traditional inlet and outlet pressure boundaries, a fully developed flow boundary could obtain more accurate numerical results, which is also consistent with the continuous property of fluid flow and has better convergence performance in CFD.

Mesh
Moving mesh technology is applied in this TEHD model, and the whole fluid flow domain is set as a moving mesh domain. During the calculation process, as the structure deformation is coupled, the shape of the mesh is adapted accordingly. Due to the fact that the velocity changes faster in the film thickness direction (axial direction) than that in the circumferential and radial directions, the layer number of the mesh in the direction of film thickness becomes more important than in the other two directions. To improve calculation efficiency with guaranteed accuracy, the convergence of mesh with the verify method described in Ref. [30] is illustrated in Figure 8. It is concluded that the simulation results converge when the mesh layer in the direction of gas film is increased to six and ensures the number of the mesh is sufficient in the circumferential and radial directions, which is reliable enough to applied to the following analysis. The fluid region is structured to be In the present study, the fully developed flow boundary condition is applied to the inlet and outlet of channel, and this type of boundary condition has been applied to the inlet boundary of a two-pass rectangular smooth channel in which velocity and temperature are mapped from the outlet of a periodic segment [27]. Qahtan [28] applied a fully developed flow boundary condition at the inlet channel of a wavy channel partially filled with a porous layer, wherein the velocities and turbulent quantities were calculated separately. Batista [29] applied the fully developed water inlet boundary condition in a crossflow air-to-water fin-and-tube heat exchanger, and compared it with uniform flow velocity and temperature profiles; its numerical results coincide well with experimental measurements. Compared to traditional inlet and outlet pressure boundaries, a fully developed flow boundary could obtain more accurate numerical results, which is also consistent with the continuous property of fluid flow and has better convergence performance in CFD.

Structure Deformation
For the structure domain, as presented in Figure 2a, boundary conditions are set as follows. The leading edge of the top foil and bump foil are regarded as fixed boundaries, and other edges are set as free; the left and right sides of the backing plate are set as rotational periodic symmetry boundaries, and the backing plate and spacer shim are regarded as fixed rigid bodies.

Thermal Transfer
The thermal computational domain contains the whole fluid flow and structure domain, as shown in Figure 6. The rotational periodic symmetry boundaries are the same as those two domains above. The inlet and outlet boundaries of the fluid flow domain are set as open boundaries. The upper side of the top foil is coupled with the bottom side of fluid flow, which results in a continuous temperature field. Besides the thermal contact region, the rotor disk side and all the rest of the solid surfaces are set as heat convection boundaries.

Mesh
Moving mesh technology is applied in this TEHD model, and the whole fluid flow domain is set as a moving mesh domain. During the calculation process, as the structure deformation is coupled, the shape of the mesh is adapted accordingly. Due to the fact that the velocity changes faster in the film thickness direction (axial direction) than that in the circumferential and radial directions, the layer number of the mesh in the direction of film thickness becomes more important than in the other two directions. To improve calculation efficiency with guaranteed accuracy, the convergence of mesh with the verify method described in Ref. [30] is illustrated in Figure 8. It is concluded that the simulation results converge when the mesh layer in the direction of gas film is increased to six and ensures the number of the mesh is sufficient in the circumferential and radial directions, which is reliable enough to applied to the following analysis. The fluid region is structured to be meshed through mesh control domain methods with approximately 520,000 cells, and the structure region is meshed with approximately 50,000 cells. Almost all the fluid and structural domains are gridded by structured meshes, i.e., hexahedral meshes; only the bump foil part is gridded by unstructured meshes, i.e., tetrahedral meshes. The structured mesh can conveniently achieve the boundary fitting, which is suitable for the calculation of fluid and surface stress concentrations, etc. It is generated quickly and with good quality, but it is only applicable to the regular shape of the structure. Generally speaking, the computational results of structured meshes converge more easily than unstructured meshes, and unstructured meshes have structural adaptability.
Lubricants 2022, 10, x FOR PEER REVIEW 9 of 20 meshed through mesh control domain methods with approximately 520,000 cells, and the structure region is meshed with approximately 50,000 cells. Almost all the fluid and structural domains are gridded by structured meshes, i.e., hexahedral meshes; only the bump foil part is gridded by unstructured meshes, i.e., tetrahedral meshes. The structured mesh can conveniently achieve the boundary fitting, which is suitable for the calculation of fluid and surface stress concentrations, etc. It is generated quickly and with good quality, but it is only applicable to the regular shape of the structure. Generally speaking, the computational results of structured meshes converge more easily than unstructured meshes, and unstructured meshes have structural adaptability.  Figure 9 shows the simulation results of the fluid flow trend compared with the results from Ref. [8], which all indicate clearly that most of the flow leaks out at the slope region side of the film. Figure 9b applied the normal pressure inlet and outlet boundary condition, while Figure 9c applied the fully developed pressure inlet and outlet boundary condition. The overall trend is roughly the same, but the flow becomes more complex at the inside slope region of the gas film.

Validation of CFD Model
The fluid structure coupled simulation results in the CFD model are compared with the results calculated by the numerical simulation method mentioned in Ref. [31], using the GFTBs geometry and operating parameters in Table 1. As Figure 10 shows, there is a positive linear trend in the relationship between bearing load-carrying capacity and rotational speed in both numerical results and CFD results. Meanwhile, the CFD results are higher than the results of the numerical simulation, which might be caused by the simplifications of the Reynolds equation and structure deformation in the numerical method. Therefore, the CFD model built in this study is reliable for subsequent analysis of the performances of GFTBs.  Figure 9 shows the simulation results of the fluid flow trend compared with the results from Ref. [8], which all indicate clearly that most of the flow leaks out at the slope region side of the film. Figure 9b applied the normal pressure inlet and outlet boundary condition, while Figure 9c applied the fully developed pressure inlet and outlet boundary condition. The overall trend is roughly the same, but the flow becomes more complex at the inside slope region of the gas film.   The fluid structure coupled simulation results in the CFD model are compared with the results calculated by the numerical simulation method mentioned in Ref. [31], using the GFTBs geometry and operating parameters in Table 1. As Figure 10 shows, there is a positive linear trend in the relationship between bearing load-carrying capacity and rotational speed in both numerical results and CFD results. Meanwhile, the CFD results are higher than the results of the numerical simulation, which might be caused by the simplifications of the Reynolds equation and structure deformation in the numerical method. Therefore, the CFD model built in this study is reliable for subsequent analysis of the performances of GFTBs.   Table 1 lists the bearing structure and operating parameters used in the present simulation.   Table 1 lists the bearing structure and operating parameters used in the present simulation.  Figure 11 illustrates the pressure distribution of gas film obtained by the fluid structure coupled simulation, with a rotational speed of 30,000 rpm for different average pressures p 0 (Pa) at inlet and outlet boundaries which are all set as fully developed flow. Figure 11a shows the pressure distribution at p 0 = 6000 Pa. Along the circumferential direction, the pressure increases rapidly from the groove symmetry boundaries at inlet, and reaches the maximum value when the gas flow reaches the junction of slope section and flat section. Then, it decreases quickly to the same pressure value with the inlet when it reaches groove symmetry boundary at outlet. In the radial direction, the high-pressure area is concentrated near the connecting line of the slope section and the flat section and close to the outer radius region. The pressure in the area closest to the inlet and outlet boundaries is close to the atmospheric pressure, but the boundary pressure near the highest pressure at the outer radius is slightly higher than atmospheric pressure due to the fully developed flow. concentrated near the connecting line of the slope section and the flat section and close to the outer radius region. The pressure in the area closest to the inlet and outlet boundaries is close to the atmospheric pressure, but the boundary pressure near the highest pressure at the outer radius is slightly higher than atmospheric pressure due to the fully developed flow.

Fluid-Structure Coupled Simulation
Due to the pressure distribution discipline with different inlet and outlet average boundary pressure are the same, as described in the previous paragraph. To show the difference in their numerical values more clearly, Figure 11b presents the pressure along the circumferential sectional line at the pressure highest point of the gas film where r = 17 mm, approximately, with different inlet and outlet average boundary pressures. It can be seen that the pressure increases almost uniformly in the form of an arithmetic progression on this line as the boundary pressure increases.  Figure 12 shows the deformations of top foil with different average boundary pressures at inlet and outlet listed in Figure 11b. The maximum top foil deformation occurs at the outer radius edge of the slope section in most instances, whereas the boundary pressure is high enough, as shown in Figure 12b-d, which is the consequence of no bump foil support under the slope section. However, in the case of insufficient boundary pressure, as shown in Figure 12a, the maximum top foil deformation appears at the outer radius edge near the trailing edge, which is the location that often first appears as an abrasion. Hence, the initial average boundary pressure has a significant effect on the deformation of the top foil. Due to the pressure distribution discipline with different inlet and outlet average boundary pressure are the same, as described in the previous paragraph. To show the difference in their numerical values more clearly, Figure 11b presents the pressure along the circumferential sectional line at the pressure highest point of the gas film where r = 17 mm, approximately, with different inlet and outlet average boundary pressures. It can be seen that the pressure increases almost uniformly in the form of an arithmetic progression on this line as the boundary pressure increases. Figure 12 shows the deformations of top foil with different average boundary pressures at inlet and outlet listed in Figure 11b. The maximum top foil deformation occurs at the outer radius edge of the slope section in most instances, whereas the boundary pressure is high enough, as shown in Figure 12b-d, which is the consequence of no bump foil support under the slope section. However, in the case of insufficient boundary pressure, as shown in Figure 12a, the maximum top foil deformation appears at the outer radius edge near the trailing edge, which is the location that often first appears as an abrasion. Hence, the initial average boundary pressure has a significant effect on the deformation of the top foil.
Under the action of the pressure transmitted by top foil, the deformation of the bump foil as a support structure was also investigated, as shown in Figure 13. When the inlet and outlet average boundary pressure is relatively small, the free end of the bump foil will be upturned, and the deformation is more significant than the front bumps, which might be caused by the peak pressure of gas film acting on the front bumps (combining the pressure distribution in Figure 11). The leverage of the force causes the rear bump to upturn, and the pressure of the trailing edge is not large enough to push it down, as Figure 13a shows. When the inlet and outlet average boundary pressure increases, the upturn displacement decreases gradually, and the first bump at the region of maximum pressure is pushed downwards more evidently; meanwhile, the following two bumps (2 nd and 3 rd ) will be deformed upward under the leverage of the force. The last two bumps (4 th and 5 th ) are pushed downwards by the top foil under the large boundary pressure. Under the action of the pressure transmitted by top foil, the deformation of the bump foil as a support structure was also investigated, as shown in Figure 13. When the inlet and outlet average boundary pressure is relatively small, the free end of the bump foil will be upturned, and the deformation is more significant than the front bumps, which might be caused by the peak pressure of gas film acting on the front bumps (combining the pressure distribution in Figure 11). The leverage of the force causes the rear bump to upturn, and the pressure of the trailing edge is not large enough to push it down, as Figure  13a shows. When the inlet and outlet average boundary pressure increases, the upturn displacement decreases gradually, and the first bump at the region of maximum pressure is pushed downwards more evidently; meanwhile, the following two bumps (2 nd and 3 rd ) will be deformed upward under the leverage of the force. The last two bumps (4 th and 5 th ) are pushed downwards by the top foil under the large boundary pressure. The effect of rotational speed on the pressure of GFTB is also investigated, as in Figure 14. The pressure distributions along the circumferential sectional line is highest pressure point where r = 17 mm. The results indicates that the pressure ri increases as the rotational speed increases and the positions of the peak pressure not change with the rotational speed. The effect of rotational speed on the pressure of GFTB is also investigated, as shown in Figure 14. The pressure distributions along the circumferential sectional line is at the highest pressure point where r = 17 mm. The results indicates that the pressure rise rate increases as the rotational speed increases and the positions of the peak pressures does not change with the rotational speed. The effect of rotational speed on the pressure of GFTB is also investigated, as shown in Figure 14. The pressure distributions along the circumferential sectional line is at the highest pressure point where r = 17 mm. The results indicates that the pressure rise rate increases as the rotational speed increases and the positions of the peak pressures does not change with the rotational speed.

Fluid-Structure-Thermal Coupled Simulation
The thermal investigation is conducted to analyze the effect of increased boundary pressure and rotational speed in this section. Under the rotational speed of 30,000 rpm,

Fluid-Structure-Thermal Coupled Simulation
The thermal investigation is conducted to analyze the effect of increased boundary pressure and rotational speed in this section. Under the rotational speed of 30,000 rpm, Figure 15 shows the pressure along the circumferential sectional line at highest pressure point where r = 17 mm, with different boundary pressures. Compared with that result in Figure 11b, it has a higher overall pressure curve when considering the thermal influences.
Lubricants 2022, 10, x FOR PEER REVIEW 14 of 20 Figure 15 shows the pressure along the circumferential sectional line at highest pressure point where r = 17 mm, with different boundary pressures. Compared with that result in Figure 11b, it has a higher overall pressure curve when considering the thermal influences. Correspondingly, the load-carrying capacity of GFTBs would also change with the different boundary pressures and rotational speeds, as illustrated in Figure 16. It can be concluded that the load-carrying capacity tends to increase linearly with the growth of rotational speed in the isothermal model; meanwhile, the increase is more pronounced at larger boundary pressure p0. When considering the influence of temperature, the loadcarrying capacity of GFTBs can be improved to some extent compared with the isothermal model, and as the rotational speed increased, the improvement of the load capacity also rose. It is an approximately ten percent increase where rotational speed is 30,000 (rpm), which may be caused by the increase in temperature. Figure 17 shows the maximum temperature influenced by boundary pressures. The change of boundary pressure generally causes a minor impact on the maximum temperature of GFTB components (maximum about 6%). The temperature of bump foil and backing plate rise a few degrees when the boundary pressure increases. On the other hand, the temperature of gas film and top foil decreased by a few degrees, which might be caused by the fully developed boundary condition that accelerates the heat exchange of boundaries.
In comparison, as shown in Figure 18, the rotational speed significantly influences the maximum temperature of GFTB due to the heat generated within the gas film by viscous shear effect, which is directly related to the gradients of fluid velocity in gas film thickness direction. The temperatures of gas film and top foil increase rapidly with the increase of rotational speed, and the temperatures of bump foil and backing plate also increase at a relatively lower rate. This is because when the heat flow transfers from gas film to the top foil, bump foil, and backing plate, on one hand the thermal contact resistant reduces the heat transfer efficiency and on the other part of the heat is lost through natural convection from structure surfaces.  Correspondingly, the load-carrying capacity of GFTBs would also change with the different boundary pressures and rotational speeds, as illustrated in Figure 16. It can be concluded that the load-carrying capacity tends to increase linearly with the growth of rotational speed in the isothermal model; meanwhile, the increase is more pronounced at larger boundary pressure p 0 . When considering the influence of temperature, the loadcarrying capacity of GFTBs can be improved to some extent compared with the isothermal model, and as the rotational speed increased, the improvement of the load capacity also rose. It is an approximately ten percent increase where rotational speed is 30,000 (rpm), which may be caused by the increase in temperature. Figure 17 shows the maximum temperature influenced by boundary pressures. The change of boundary pressure generally causes a minor impact on the maximum temperature of GFTB components (maximum about 6%). The temperature of bump foil and backing plate rise a few degrees when the boundary pressure increases. On the other hand, the temperature of gas film and top foil decreased by a few degrees, which might be caused by the fully developed boundary condition that accelerates the heat exchange of boundaries.       In comparison, as shown in Figure 18, the rotational speed significantly influences the maximum temperature of GFTB due to the heat generated within the gas film by viscous shear effect, which is directly related to the gradients of fluid velocity in gas film thickness direction. The temperatures of gas film and top foil increase rapidly with the increase of rotational speed, and the temperatures of bump foil and backing plate also increase at a relatively lower rate. This is because when the heat flow transfers from gas film to the top foil, bump foil, and backing plate, on one hand the thermal contact resistant reduces the heat transfer efficiency and on the other part of the heat is lost through natural convection from structure surfaces.   Figure 19. The left side is under the rotational speed of 6000 rpm, and the right side is under the rotational speed of 30,000 rpm, and both have the inlet and outlet boundary pressure of p0 = 6000 Pa. Correspondingly, Figure 20 presents the  Figure 19. The left side is under the rotational speed of 6000 rpm, and the right side is under the rotational speed of 30,000 rpm, and both have the inlet and outlet boundary pressure of p 0 = 6000 Pa. Correspondingly, Figure 20 presents the characteristic temperature distribution of GFTB on the same condition. It can be seen from Figure 20a,b that the temperature is distributed evenly in the thin gas film and accumulated at the outer radial trailing edge, which is because the rotor disk is simplified as a natural convection. Moreover, the lowest temperature region appears at the groove of the inner and outer radial boundary where the ambient air enters, which means the inner and outer radial boundaries of the groove are the actual gas inlet boundaries, and the rest of pressure boundaries are the actual gas outlet boundaries. This is due to the fact that the backflow is not suppressed at the pressure inlet and outlet boundaries.
The temperature distributions of top foil at different rotational speeds are similar. Figure 20c,d show that the highest temperature region appears at the outer radial trailing edge of top foil, which is the same as the results in Ref. [11]. Considering the fluid flow diagram in Figure 9c, most of the flow leaks out at the slope region of the gas film, which also takes away some heat leading to a low-temperature area at the slope section of top foil.
The temperature rise of bump foil and backing plate is mainly due to thermal contacts (i.e., bump foil/top foil and bump foil/backing plate), which contain air gap heat conduction and a diffuse radiation of gray bodies between parallel planes associated with surface roughness, contact pressure, air gap properties, etc. It can be clearly seen that the heat transfer through the contact pairs in Figure 20e-h, combines the pressure distribution and foil deformation in Figure 19, and the heat transfer efficiency of the contact pairs are greatly affected by the contact pressure. There is also a clear line of a temperature drop on the top foil, at the contact pair coincidence area where the heat transfer efficiency is most pronounced especially in Figure 20d. lated at the outer radial trailing edge, which is because the rotor disk is simplified as a natural convection. Moreover, the lowest temperature region appears at the groove of the inner and outer radial boundary where the ambient air enters, which means the inner and outer radial boundaries of the groove are the actual gas inlet boundaries, and the rest of pressure boundaries are the actual gas outlet boundaries. This is due to the fact that the backflow is not suppressed at the pressure inlet and outlet boundaries. The temperature distributions of top foil at different rotational speeds are similar. Figure 20c,d show that the highest temperature region appears at the outer radial trailing edge of top foil, which is the same as the results in Ref. [11]. Considering the fluid flow diagram in Figure 9c, most of the flow leaks out at the slope region of the gas film, which also takes away some heat leading to a low-temperature area at the slope section of top foil.
The temperature rise of bump foil and backing plate is mainly due to thermal contacts (i.e., bump foil/top foil and bump foil/backing plate), which contain air gap heat conduction and a diffuse radiation of gray bodies between parallel planes associated with surface roughness, contact pressure, air gap properties, etc. It can be clearly seen that the heat transfer through the contact pairs in Figure 20e-h, combines the pressure distribution and foil deformation in Figure 19, and the heat transfer efficiency of the contact pairs are greatly affected by the contact pressure. There is also a clear line of a temperature drop on the top foil, at the contact pair coincidence area where the heat transfer efficiency is most pronounced especially in Figure 20d. The simulated top foil temperature is compared with other similar research results as shown in Figure 21. It can be seen that the temperature distribution trends are similar, and the temperature peak both appears at the trailing edge near out radius, but the overall temperature distributions in reference [11] is relatively higher than the results in this work, which is mainly caused by the facts that the rotor speed in reference is much higher, and the size of the bearing is also relatively larger. Thus, the TEHD model built in present work is reasonable for the analysis.
The simulated top foil temperature is compared with other similar research results as shown in Figure 21. It can be seen that the temperature distribution trends are similar, and the temperature peak both appears at the trailing edge near out radius, but the overall temperature distributions in reference [11] is relatively higher than the results in this work, which is mainly caused by the facts that the rotor speed in reference is much higher, and the size of the bearing is also relatively larger. Thus, the TEHD model built in present work is reasonable for the analysis.  Top foil temperature comparison with Ref. [11]. (a) Top foil temperature increase for a rotor speed of 120 krpm in Ref. [11]. (b) Top foil temperature distribution in present model.

Conclusions
A detailed three-dimensional thermo-elastic-hydrodynamic (TEHD) coupling model has been developed to analyze the fluid-structure-thermal performance of bump-type GFTBs in the current study. The fully developed flow boundary conditions were applied Figure 21. Top foil temperature comparison with Ref. [11]. (a) Top foil temperature increase for a rotor speed of 120 krpm in Ref. [11]. (b) Top foil temperature distribution in present model.