The Direct-Coupling Method for Analyzing the Performance of Aerostatic Bearings Considering the Fluid–Structure Interaction Effect

: In the interest of analyzing the effect of the structural deformation root caused by gas pressure on the static features of aerostatic bearings, a ﬂuid–structure interaction (FSI) model based on oriﬁce-type aerostatic bearings is proposed that can predict the characteristics of aerostatic bearings more accurately by using the direct-coupling method (DCM). By using COMSOL Multiphysics, the governing equation matrix of the ﬁnite element model of structural deformation and gas ﬁlm pressure was solved with the integral solution method, and the oriﬁce boundary conditions were calculated with the root iteration method. At the same time, the static performance of I-shaped oriﬁce-type aerostatic bearing with various supply pressures was analyzed theoretically and tested experimentally. The results show that in comparison with the calculation results without taking account of structural deformation, the theoretical values from the model derived in this paper considering the FSI effect are closer to the experimental values. Finally, by using the orthogonal design method, FSI simulation was carried out to analyze how the key dimension factors inﬂuence the structural stiffness of the spindle, and it is concluded that the thrust bearing’s stiffness is strongly inﬂuenced by the thickness of the thrust plate.


Introduction
Due to the absence of friction, low heat generation and the good error equalization effect of the gas film, aerostatic bearings provide high rotation precision and low vibration. In ultra-precision machines and precision measuring instruments, aerostatic bearings are commonly used [1][2][3][4]. Compared with hydrostatic bearings, the biggest disadvantage of aerostatic bearings is their lower stiffness. As a part of the cutting system, the stiffness of the bearing affects the overall stiffness of the system [5]. Therefore, the most critical technical index for the design of the static gas pressure bearing is its stiffness. However, in the current design of the aerostatic bearing, there is a large difference between the theoretical design values and experimental test results. During the theoretical analysis process, the spindle structure is assumed to be rigid, and the effect of the flexible deformation of the spindle on the bearing performance is not considered. Therefore, when modeling and analyzing aerostatic bearings, the fluid-structure interaction (FSI) effect caused by structural deformation must be fully considered [6].
To date, many scholars have studied the FSI effect of aerostatic bearings. Anton van Beek et al. [7] analyzed how the performance was affected by throttled aerostatic bearings with elastically deformed surfaces by simultaneously solving the 2D Reynolds equation and 3D elastic equation of compressible fluids. Nagendra et al. [8] imported the geometric model into ANSYS Mechanical and conducted Computational Fluid Dynamics (CFD) simulations to analyze the influence of three different materials' bushings on the gas film pressure distribution. Chen et al. [9] proposed a bearing performance calculation method incorporating the FSI effect and used the gas film thickness after structural deformation to calculate the actual bearing performance. Dhande et al. [10] considered the actual bearing deformation through bidirectional fluid-structure coupling (FSI) and cavitation and concluded that cavitation would reduce the maximum bearing pressure. Based on their proposed comprehensive bidirectional FSI model, Lu et al. [11,12] conducted an investigation into the effects of structural parameters on static bearing performance and the thermal influence of gas film viscosity. Chen et al. [13] established a multiphase flow model for radial bearings that considered both heat and cavitation effects and analyzed the influences of rotational speed, eccentricity and film thickness on the elasto-hydrodynamics of journal bearings with variously shaped grooves. Yan et al. [14] established the FSIW model based on a porous aerostatic bearing and analyzed the FSI deformation law of different porous materials. Wang et al. [15] studied two journal bearings by using the free fluid-structure coupling method (FFSIM); they carried out a bifurcation analysis, studied the influence of rotational speed and rotor mass on the shaft track and generated a phase diagram and frequency response curve. Wei Wang et al. [16] analyzed the effects of gas supply pressure, Young's modulus and porous material thickness on the bearing's static performance. Gao et al. [17,18] proposed a two-wheel optimization design method for the bearing stiffness and key structure size considering FSI. Liu et al. [19] applied the CFD-FSI method to study the dynamic response of the system based on a series of dynamically unbalanced loads and various bearing materials.
The core work of analyzing aerostatic bearings' static performance is mainly to solve the specific Reynolds equation boundary conditions, thus revealing the law of the distribution of the gas film pressure and obtaining the supporting aerostatic bearing performance. However, the Reynolds equation characterized by gas lubrication is a second-order partial differential equation, which makes it difficult to obtain an exact analytical solution. Therefore, the gas lubrication problem is usually solved based on computational fluid dynamics (CFD), while the finite element method is usually used in the analysis of structural deformation. Therefore, when analyzing FSI problems, it is necessary to use different solving modules to calculate the governing equations for the fluid and solid domains and update the boundary conditions in each iteration according to different regional characteristics [20]. However, compared with the monolithic method, due to the non-conservation of energy caused by the calculation lag of the gas domain compared to the solid domain, the above method does not converge at the solid-gas interface [21]. The immersed boundary-lattice Boltzmann method is another Eulerian-Lagrangian approach used for FSI simulations, even in flexible boundary problems [22,23]. In the analysis of FSI problems, the above scholars adopted the sequential solution method to clearly analyze the gas film pressure distribution and structural deformation. However, we find that there are two obvious defects when using the sequential coupling method to solve the FSI problem for aerostatic bearings. Firstly, iteration parameters need to be constantly adjusted in the process of iteration because the calculation lag of the gas domain and solid domain tends to cause the non-convergence of iterations. Secondly, repeated modeling is required when analyzing different structural design parameters. Since the Reynolds equation represents a partial differential equation having a second order, it is also applicable to the finite element method [24]. Therefore, if the elastic deformation equation and Reynolds equation are discretized and formed into a matrix, the accuracy and efficiency of the calculation will be greatly improved.
In this paper, we propose a method for modeling FSI based on the DCM to foresee and understand the static features of aerostatic bearings. In order to verify the accuracy of the theoretical analysis, the static characteristics were tested under different gas supply pressures. In addition, we also investigated how critical structural parameters influence the static performance of thrust bearings, providing the basis for the optimal design.  Figure 1a,b show the model of the aerostatic bearing structure. The bearing adopts an I-shaped closed structure, which is composed of upper and lower thrust bearings of different sizes, radial bearings and shaft sleeves with orifice restrictors and air supply pipes. The upper, lower and journal bearings adopt three independent air supply channels. Under working conditions, the high-pressure gas in the gas supply pipeline enters the gas film area through the throttling effect of the restrictor, and after that, it returns to the air through the gas film boundary [25].  Figure 1a,b show the model of the aerostatic bearing structure. The bearing adopt an I-shaped closed structure, which is composed of upper and lower thrust bearings o different sizes, radial bearings and shaft sleeves with orifice restrictors and air supply pipes. The upper, lower and journal bearings adopt three independent air supply chan nels. Under working conditions, the high-pressure gas in the gas supply pipeline enter the gas film area through the throttling effect of the restrictor, and after that, it returns to the air through the gas film boundary [25].  Figure 1c accounts for the interaction between the spindle and gas film under high pressure external gas. F1 and F2 are the forces exerted by the gas film on two thrust plates upper and lower, respectively. F3 is the radial gas film force. As a result of the gas film force, the mandrel is mainly stretched by F3. The lower and upper thrust plates are warped by the bending moment under the action of F1 and F2 while being affected by the tension of the mandrel. Finally, the spindle in Figure 1c undergoes structural deformation as a result of the gas film force. At the same time, structural deformation leads to changes in the gap between the spindle and the bushing, which changes the actual thickness and ga film shape, thereby affecting the static performance of the bearing.

The Governing Equations of FSI Systems
As shown in Figure 2, the aerostatic bearing FSI model in the static state is described Vf and Vs denote the fluid region (gas film) and solid region (spindle), respectively. S0 rep resents the solid-gas boundary, Sb corresponds to the fixed displacement boundaries, and Sσ is the external load application interface. P and P0 are the gas film area pressure and ambient pressure boundary, respectively, and h is the gas film thickness. ns and nf are a pair of unit normal vectors with the same size and opposite directions on the solid-ga boundary line, which represent the unit normal vectors outside the solid and fluid bound ary elements, respectively.  Figure 1c accounts for the interaction between the spindle and gas film under highpressure external gas. F 1 and F 2 are the forces exerted by the gas film on two thrust plates, upper and lower, respectively. F 3 is the radial gas film force. As a result of the gas film force, the mandrel is mainly stretched by F 3 . The lower and upper thrust plates are warped by the bending moment under the action of F 1 and F 2 while being affected by the tension of the mandrel. Finally, the spindle in Figure 1c undergoes structural deformation as a result of the gas film force. At the same time, structural deformation leads to changes in the gap between the spindle and the bushing, which changes the actual thickness and gas film shape, thereby affecting the static performance of the bearing.

The Governing Equations of FSI Systems
As shown in Figure 2, the aerostatic bearing FSI model in the static state is described. V f and V s denote the fluid region (gas film) and solid region (spindle), respectively. S 0 represents the solid-gas boundary, S b corresponds to the fixed displacement boundaries, and S σ is the external load application interface. P and P 0 are the gas film area pressure and ambient pressure boundary, respectively, and h is the gas film thickness. n s and n f are a pair of unit normal vectors with the same size and opposite directions on the solidgas boundary line, which represent the unit normal vectors outside the solid and fluid boundary elements, respectively. To establish a static performance analysis model of the aerostatic bearing considering the influence of structural deformation, it is necessary to conduct an FSI analysis on the two processes of structural deformation and the aerostatic bearing. For the analysis of the gas film region, it is assumed that the gas is an agglutinant and incompressible gas and in a laminar flow state, while the material in the solid region is a linear elastic material. When the main shaft rotates slowly or is stationary, the motion of the gas film region satisfies the static Reynolds equation. For the analysis of solid regions based on elastic mechanics, it is necessary to consider deformation as a variable to establish the static equilibrium differential equation. After introducing the gradient operator into the FSI system, the partial differential control equation is obtained as: (1) in which σ is the solid strain tensor, b is the body force vector on the surface of the solid region, the pressure of the gas is p, and the equation constant is 3 12 , where η and h are the aerodynamic viscosity and the gas film thickness, respectively. Based on the assumption of linear elasticity, the constitutive equation can be further simplified as: where C represents a matrix of linear elastic stiffness, and ε denotes the strain tensor.
In small deformations, the relationship between displacement u and strain ε is: The boundary conditions at the solid-gas interface should meet the following requirements: (a) Kinematic boundary conditions: In the gas film region, when the gas enters, the velocity slip phenomenon is ignored; that is, as the gas enters the gas film region, its normal velocity remains continuous. At the same time, it is assumed that the Reynolds equation describes the velocity field in the gas film region. During the static analysis, the velocity at the junction of the solid and gas is set to zero. On the other hand, in the dynamic model, the gas motion must be described by the dynamic Reynolds equation. (b) Dynamic boundary conditions: The normal force is continuous at the junction of the solid and gas. At the same time, for the solid domain, the following three boundary conditions should also be met: To establish a static performance analysis model of the aerostatic bearing considering the influence of structural deformation, it is necessary to conduct an FSI analysis on the two processes of structural deformation and the aerostatic bearing. For the analysis of the gas film region, it is assumed that the gas is an agglutinant and incompressible gas and in a laminar flow state, while the material in the solid region is a linear elastic material. When the main shaft rotates slowly or is stationary, the motion of the gas film region satisfies the static Reynolds equation. For the analysis of solid regions based on elastic mechanics, it is necessary to consider deformation as a variable to establish the static equilibrium differential equation. After introducing the gradient operator into the FSI system, the partial differential control equation is obtained as: in which σ is the solid strain tensor, b is the body force vector on the surface of the solid region, the pressure of the gas is p, and the equation constant is k = h 3 12η , where η and h are the aerodynamic viscosity and the gas film thickness, respectively. Based on the assumption of linear elasticity, the constitutive equation can be further simplified as: where C represents a matrix of linear elastic stiffness, and ε denotes the strain tensor. In small deformations, the relationship between displacement u and strain ε is: The boundary conditions at the solid-gas interface should meet the following requirements: (a) Kinematic boundary conditions: In the gas film region, when the gas enters, the velocity slip phenomenon is ignored; that is, as the gas enters the gas film region, its normal velocity remains continuous. At the same time, it is assumed that the Reynolds equation describes the velocity field in the gas film region. During the static analysis, the velocity at the junction of the solid and gas is set to zero. On the other hand, in the dynamic model, the gas motion must be described by the dynamic Reynolds equation. (b) Dynamic boundary conditions: The normal force is continuous at the junction of the solid and gas. At the same time, for the solid domain, the following three boundary conditions should also be met: where u denotes the boundary condition of displacement, and p represents the boundaries of traction. The pressure at the gas boundary should meet the following conditions: p = P d in gas chamber p = P 0 on the edge of atmosphere boundary (5) In Equation (5), P d represents the boundary pressure at the exit of the orifice, and P 0 represents the boundary pressure where the gas film meets the atmosphere.
Most of the currently published references on Finite Element Method (FEM) are calculated for either displacement or pressure in a discrete form. In the Appendix A of this paper, using finite element analysis, we establish a discrete linear equation system to directly resolve the FSI issue and provide the finite element solution expression.

The Procedure for Solving FSI Model Based on DCM
From the above analysis, we can see that the structural deformation and gas film pressure are mutually coupled and affect each other. Due to the limited rigidity of the material of the spindle, deformation is caused by the pressure of the gas film. In the case of the spindle structural deformation, the thickness of the gas film will change at the deformed location, thereby affecting the pressure distribution in the gas film. Since structural deformation and film pressure changes occur simultaneously, both processes must be solved simultaneously. At the same time, due to the radial bearing's interaction with its structure and the thrust bearing and gas film, it can be regarded as a nonlinear problem; to solve it, the Newton-Raphson iterative algorithm was used. Figure 3 shows an overview of the iterative solution process. In this process, COMSOL Multiphysics is used to solve the governing equation of the finite element model by using the global solution method. For each iteration, the initial gas film thickness h is determined after the boundary conditions have been determined, and the influence of the gas film pressure under the gas film thickness on the structural deformation is calculated. If the maximum deformation is greater than 0.1 µm, the deformed gas film thickness is used in iterative calculations until the maximum structural deformation is less than 0.1 µm. In Equation (5), Pd represents the boundary pressure at the exit of the orifice, and P0 represents the boundary pressure where the gas film meets the atmosphere.
Most of the currently published references on Finite Element Method (FEM) are calculated for either displacement or pressure in a discrete form. In the appendix of this paper, using finite element analysis, we establish a discrete linear equation system to directly resolve the FSI issue and provide the finite element solution expression.

The Procedure for Solving FSI Model Based on DCM
From the above analysis, we can see that the structural deformation and gas film pressure are mutually coupled and affect each other. Due to the limited rigidity of the material of the spindle, deformation is caused by the pressure of the gas film. In the case of the spindle structural deformation, the thickness of the gas film will change at the deformed location, thereby affecting the pressure distribution in the gas film. Since structural deformation and film pressure changes occur simultaneously, both processes must be solved simultaneously. At the same time, due to the radial bearing's interaction with its structure and the thrust bearing and gas film, it can be regarded as a nonlinear problem; to solve it, the Newton-Raphson iterative algorithm was used. Figure 3 shows an overview of the iterative solution process. In this process, COM-SOL Multiphysics is used to solve the governing equation of the finite element model by using the global solution method. For each iteration, the initial gas film thickness h is determined after the boundary conditions have been determined, and the influence of the gas film pressure under the gas film thickness on the structural deformation is calculated. If the maximum deformation is greater than 0.1 μm, the deformed gas film thickness is used in iterative calculations until the maximum structural deformation is less than 0.1 μm.   Figure 4 shows an orifice-type aerostatic thrust bearing, in which an ε-deep pressureequalizing groove is specially set up to improve rigidity. The high-pressure gas flows into the gas film area from the orifice and finally flows into the surrounding environment from  Figure 4 shows an orifice-type aerostatic thrust bearing, in which an ε-deep pressureequalizing groove is specially set up to improve rigidity. The high-pressure gas flows into the gas film area from the orifice and finally flows into the surrounding environment from the thrust plate boundary. Using the finite element method, the Reynolds equation can be solved to calculate the orifice aerostatic thrust bearing's stiffness and the bearing capacity. From Equation (6), it is easy to see that the flow rate is correlated with the pressure drop in the orifice [26]:

The Calculation Procedure
where ∅ represents the restriction coefficient, which is usually 0.8; the area of the orifice restrictor is defined by A 0 ; the supply pressure is defined by P 0 ; the density of the air at standard atmospheric pressure is denoted by ρ a ; P a is standard atmospheric pressure; and the flow function φ is [26]: where k is the specific heat ratio of the gas, β is the throttling pressure ratio, and the critical throttle pressure ratio β k is [26]: As mentioned above, the FSI model mainly includes the direct-coupling calculation of the aerostatic bearing's boundary conditions and iterative convergence accuracy improvement. Using the direct-coupling calculation method, a FEM-based model was used to simulate the distribution of gas film pressure in the aerostatic bearings and a full-scale spindle structural FEM-based model. The structural deformation of the main shaft and the gas film pressure distribution can be calculated by using COMSOL Multiphysics. However, the establishment of the direct-coupling FSI model mainly describes the characteristics of the gas film region while ignoring the characteristics of the orifice. This is because the orifice characteristics are nonlinear. In this context, the relationship between the flow rate and pressure difference can be examined.
To ensure the convergence of the calculation process, the orifice characteristics of the specific boundary conditions must be correctly characterized. Therefore, this work performed joint calculations in MATLAB and COMSOL and introduced a programmable program to realize nonlinear iteration to characterize the flow balance at the inlet and outlet of the orifice. From Equation (6), it is easy to see that the flow rate is correlated with the pressure drop in the orifice [26]: where ∅ represents the restriction coefficient, which is usually 0.8; the area of the orifice restrictor is defined by A 0 ; the supply pressure is defined by P 0 ; the density of the air at standard atmospheric pressure is denoted by ρ a ; P a is standard atmospheric pressure; and the flow function ϕ is [26]: where k is the specific heat ratio of the gas, β is the throttling pressure ratio, and the critical throttle pressure ratio β k is [26]: As mentioned above, the FSI model mainly includes the direct-coupling calculation of the aerostatic bearing's boundary conditions and iterative convergence accuracy improvement. Using the direct-coupling calculation method, a FEM-based model was used to simulate the distribution of gas film pressure in the aerostatic bearings and a full-scale spindle structural FEM-based model. The structural deformation of the main shaft and the gas film pressure distribution can be calculated by using COMSOL Multiphysics. However, the establishment of the direct-coupling FSI model mainly describes the characteristics of the gas film region while ignoring the characteristics of the orifice. This is because the orifice characteristics are nonlinear. In this context, the relationship between the flow rate and pressure difference can be examined.
To ensure the convergence of the calculation process, the orifice characteristics of the specific boundary conditions must be correctly characterized. Therefore, this work performed joint calculations in MATLAB and COMSOL and introduced a programmable program to realize nonlinear iteration to characterize the flow balance at the inlet and outlet of the orifice.
As can be seen in Figure 5, the entire calculation process of the proposed FSI model is observed, and the calculation steps are as follows:

1.
Firstly, it is important to establish a finite element model for the spindle structure and use the Fluid module of COMSOL to perform preprocessing, such as boundary condition setting and the mesh division of the gas film. As the initial boundary condition for the gas film region, the gas pressure flowing out of the orifice into the pressure-equalizing groove is used, and the contact between the gas film and the air is used as the outer boundary, where the outer boundary condition is set to 1 bar pressure.

2.
The FEM model adopts the global solution method to solve the pressure distribution, flow and structural deformation of the aerostatic bearing. 3.
In the COMSOL calculation in step (2), the flow rate in the gas film area is obtained, and it is introduced into MATLAB to calculate the post-orifice pressure (P ao ) of each orifice restrictor; then, the result is compared with the initially given pressureequalizing groove pressure (P op ). If the difference between them does not satisfy the convergence condition, P ao is taken as a new P op for a new round of iterations. When P ao and P op meet the convergence conditions, the iterative process ends.
As can be seen in Figure 5, the entire calculation process of the proposed FSI model is observed, and the calculation steps are as follows: 1. Firstly, it is important to establish a finite element model for the spindle structure and use the Fluid module of COMSOL to perform preprocessing, such as boundary condition setting and the mesh division of the gas film. As the initial boundary condition for the gas film region, the gas pressure flowing out of the orifice into the pressure-equalizing groove is used, and the contact between the gas film and the air is used as the outer boundary, where the outer boundary condition is set to 1 bar pressure. 2. The FEM model adopts the global solution method to solve the pressure distribution, flow and structural deformation of the aerostatic bearing. 3. In the COMSOL calculation in step (2), the flow rate in the gas film area is obtained, and it is introduced into MATLAB to calculate the post-orifice pressure (Pao) of each orifice restrictor; then, the result is compared with the initially given pressure-equalizing groove pressure (Pop). If the difference between them does not satisfy the convergence condition, Pao is taken as a new Pop for a new round of iterations. When Pao and Pop meet the convergence conditions, the iterative process ends. In order to increase the convergence speed, an effective algorithm must be developed. Currently, the root iteration method (RIM) is often used for aerostatic bearings [20]. In this paper, it is extended to the aerostatic system to improve the iterative calculation speed. Figure 6 shows the RIM with a single orifice, where the black curve indicates the flow rate through that orifice, calculated by MATLAB according to Equation (6). Equation (1) corresponds to the elastic mechanics expression for the energy expression 1 2 and its corresponding solution is related to the energy extreme. It is known that the film's mass flow rate (MFR) (q) is linearly related to the orifice outlet pressure (Pd 2 ). Therefore, the pressure P is used as the horizontal coordinate variable, and q=k +b describes the linear relationship, where k is the slope of the line. In order to increase the convergence speed, an effective algorithm must be developed. Currently, the root iteration method (RIM) is often used for aerostatic bearings [20]. In this paper, it is extended to the aerostatic system to improve the iterative calculation speed. Figure 6 shows the RIM with a single orifice, where the black curve indicates the flow rate through that orifice, calculated by MATLAB according to Equation (6). Equation (1) corresponds to the elastic mechanics expression for the energy expression 1/2f T Kf-Qf 2 , and its corresponding solution is related to the energy extreme. It is known that the film's mass flow rate (MFR) (q) is linearly related to the orifice outlet pressure (P d 2 ). Therefore, the pressure P is used as the horizontal coordinate variable, and q = kP d + b describes the linear relationship, where k is the slope of the line. From the above analysis, it is known that half of the maximum and the maximum MFR (q max ) that a single orifice can provide are used as boundary conditions. Then, COM-SOL can be used to solve the corresponding pressures Pdmax and Pdmaxhalf after the orifice. From the above analysis, it is known that half of the maximum and the maximum MFR (q max ) that a single orifice can provide are used as boundary conditions. Then, COMSOL can be used to solve the corresponding pressures P dmax and P dmaxhalf after the orifice. Finally, the expression of the gas film pressure distribution can be calculated according to the coordinates of the two points, and by combining the gas film area and the MFR formula at the orifice, the pressure and flow rate of the bearings are calculated as well.
Compared with the solution for gas lubrication with a single orifice, the outlets of multiple orifices are connected to each other through the gas film area, so the effects of deformations of different orifices are all different and affect each other. At the same time, the performance of the bearing is also affected by the change in each restrictor. During one iteration, the MFR for each orifice will be overwritten with new values, which will change the intercept, b. Therefore, RIM, as shown in Figure 7, is a more efficient algorithm. From the above analysis, it is known that half of the maximum and the maximum MFR (q max ) that a single orifice can provide are used as boundary conditions. Then, COM-SOL can be used to solve the corresponding pressures Pdmax and Pdmaxhalf after the orifice. Finally, the expression of the gas film pressure distribution can be calculated according to the coordinates of the two points, and by combining the gas film area and the MFR formula at the orifice, the pressure and flow rate of the bearings are calculated as well.
Compared with the solution for gas lubrication with a single orifice, the outlets of multiple orifices are connected to each other through the gas film area, so the effects of deformations of different orifices are all different and affect each other. At the same time, the performance of the bearing is also affected by the change in each restrictor. During one iteration, the MFR for each orifice will be overwritten with new values, which will change the intercept, b. Therefore, RIM, as shown in Figure 7, is a more efficient algorithm. Firstly, in the case of a single orifice, the maximum MFR (q max ) and half of the maximum MFR (q maxhalf ) that a single orifice can provide are used as boundary conditions. COMSOL can be used to solve the corresponding pressures Pdmax and Pdmaxhalf after the orifice. The slope k and initial intercept b0 of the Film MFR line can be calculated from the two sets of data. The gas film area and the MFR formula at the orifice are combined to calculate the pressure and flow at the operating point (i.e., q 0 and 0 ). However, because different orifices interact with each other, with multiple orifices, the q 0 and 0 calculated above do not satisfy the Film MFR equation. For this reason, in the initial calculation, we use q 0 as the MFR boundary condition and obtain the MFR equation P d0 ' , which does not satisfy the characteristics of an orifice. Therefore, while keeping the slope k constant, we perform iterative calculations based on the new Film Firstly, in the case of a single orifice, the maximum MFR (q max ) and half of the maximum MFR (q maxhal f ) that a single orifice can provide are used as boundary conditions. COMSOL can be used to solve the corresponding pressures P dmax and P dmaxhalf after the orifice. The slope k and initial intercept b 0 of the Film MFR line can be calculated from the two sets of data. The gas film area and the MFR formula at the orifice are combined to calculate the pressure and flow at the operating point (i.e., q 0 and P d0 ).
However, because different orifices interact with each other, with multiple orifices, the q 0 and P d0 calculated above do not satisfy the Film MFR equation. For this reason, in the initial calculation, we use q 0 as the MFR boundary condition and obtain the MFR equation P d0 , which does not satisfy the characteristics of an orifice. Therefore, while keeping the slope k constant, we perform iterative calculations based on the new Film MFR equation. Each pair of q i and P di+1 will get a new intercept b i+1 during the iteration. When b i+1 is almost unchanged, the iteration ends.

Simulation Results
According to the spindle structure in Figure 1, its FEM model can be established and meshed. Figure 8a-c show the FEM models with three different mesh densities. By calculating the static performance of the upper thrust bearing with different mesh densities, the calculation accuracy of the algorithm in this paper can be verified.
After running the calculation program, the simulation results considering FSI were obtained, as presented in Figure 9. During the simulation analysis, the gas supply pressure was set to 5.5 Bar, thrust bearings were designed with a gas film thickness of 23 µm on both sides, the designed gas film thickness of the journal bearing was 15 µm, the orifices measured 200 mm in diameter and the rotating shaft weighed 150 kg, and nitride steel was selected as its material.
When bi+1 is almost unchanged, the iteration ends.

Simulation Results
According to the spindle structure in Figure 1, its FEM model can be established and meshed. Figure 8a-c show the FEM models with three different mesh densities. By calculating the static performance of the upper thrust bearing with different mesh densities, the calculation accuracy of the algorithm in this paper can be verified. After running the calculation program, the simulation results considering FSI were obtained, as presented in Figure 9. During the simulation analysis, the gas supply pressure was set to 5.5 Bar, thrust bearings were designed with a gas film thickness of 23 μm on both sides, the designed gas film thickness of the journal bearing was 15 μm, the orifices measured 200 mm in diameter and the rotating shaft weighed 150 kg, and nitride steel was selected as its material. After running the calculation program, the simulation results considering FSI were obtained, as presented in Figure 9. During the simulation analysis, the gas supply pressure was set to 5.5 Bar, thrust bearings were designed with a gas film thickness of 23 μm on both sides, the designed gas film thickness of the journal bearing was 15 μm, the orifices measured 200 mm in diameter and the rotating shaft weighed 150 kg, and nitride steel was selected as its material.  Figure 9a illustrates the thrust plate's deformation after accounting for FSI effects. Under the action of the gas film force, the thrust plate bends and deforms. The deformation range of the upper thrust plate is 0.1-3.3 µm, and the deformation range of the lower thrust plate is 0.1-2.4 µm. It is at the point where the outer diameter of the upper thrust plate meets the air volume that the thrust plate experiences its greatest deformation, and the deformation amount is 3.3 µm. The minimum deformation occurs at the contact of the mandrel of the thrust plate, which is mainly due to the elongation of the shaft under the action of the gas film force, and the deformation is 0.1 µm. As shown in Figure 9b, there is a variation in the thrust bearing's gas film thickness under FSI conditions. Clearly, the change trend of the gas film thickness is the same as that of the thrust plate. At the same time, in a closed thrust bearing, the gas film clearance is the total clearance of the upper and lower thrust bearings; at the contact point between the outer diameter of the lower thrust plate and the air environment, the gas film thickness will deform the most, and the deformation is 4.8 µm (deformation of the upper thrust plate is 2.7 µm, and bottom thrust deformation is 2.1 µm).
Based on FSI, Figure 9c shows the thrust bearing's pressure distribution. It can be proven that the pressure distribution is the same as that of the traditional orifice-type aerostatic thrust bearing, and the gas film is filled with high-pressure gas from the orifice and finally decays to atmospheric pressure at the gas film boundary. Simultaneously, to analyze the FSI effect on the thrust bearing's performance, we calculated the pressure distribution in the thrust bearing without considering FSI, and Figure 9d illustrates the results. Clearly, the overall pressure of the gas film is greater without considering FSI. This is because after considering the effect of FSI, the gas film gap behind the orifice increases since the thrust plate is bent, and it deforms under the influence of the gas film force, so at the same gas supply pressure, the pressure behind the orifice decreases, leading to a decrease in the overall pressure of the gas film.
The calculation results in Figure 9 demonstrate that the structural deformation of the thrust plate is quite small, so the nonlinear factor can be ignored when modeling the solid phase, and the deformation is assumed to be linear. At the same time, in the case of small deformation, gas wetting still satisfies the assumption of the Reynolds equation.

Experimental Validation
In order to test the accuracy of the theoretical analysis, we measured the stiffness of the thrust bearing. Under different gas supply pressures, the upper thrust plate deforms, and the thrust bearing flow rate varies. As shown in Figure 10a, it is an orifice-type aerostatic bearing in the experiment. The material of the spindle is nitride steel, the thickness of the upper thrust plate is 60 mm, and the outer diameter is 400 mm. The thickness of the lower thrust plate is 45 mm, and the outer diameter is 350 mm. The outer diameter of the shaft is 200 mm. Gas supply pressure settings are 5, 5.5, 6 and 6.5 bar. In this experiment, two probes were used to measure the stiffness and deformation of the thrust bearing through an inductance micrometer (TESA-TT80 inductance sensor, measuring range ±1 mm, resolution 10 nm). The thrust bearing flow was measured with an SMC-PF2M721-02-D flow meter. Figure 10b is the device used to measure stiffness and deformation. The upper thrust plate's center and edge displacements were measured by using probe A and probe B, respectively. The thrust plate's deflection was determined by comparing the reading values of probe A and probe B.

PEER REVIEW 11 of 18
stiffness dropped significantly after considering FSI and only grew from 1293 N/μm to 2073 N/μm as the pressure increased, and the higher the supply pressure, the greater the difference. Clearly, the thickness of the gas film is greatly affected by the structural deformation caused by the gas supply pressure. With an increase in gas supply pressure, the deviation between the actual gas film thickness of the bearing and the ideal value designed by the optimal stiffness criterion increases, and the theoretical stiffness of the thrust bearing decreases correspondingly without considering the FSI effect.   Figure 11a shows the static stiffness of the thrust bearing under different gas supply pressures. It can be seen that, without considering the FSI effect, in the case of an increase in gas supply pressure from 5 Bar to 6.5 Bar, the theoretical stiffness of the thrust bearing improves from 1607 N/µm to 2738 N/µm. According to the results, the thrust bearing's stiffness dropped significantly after considering FSI and only grew from 1293 N/µm to 2073 N/µm as the pressure increased, and the higher the supply pressure, the greater the difference. Clearly, the thickness of the gas film is greatly affected by the structural deformation caused by the gas supply pressure. With an increase in gas supply pressure, the deviation between the actual gas film thickness of the bearing and the ideal value designed by the optimal stiffness criterion increases, and the theoretical stiffness of the thrust bearing decreases correspondingly without considering the FSI effect. (a) (b) Figure 11. Impacts of gas supply pressure: (a) stiffness and (b) flow rate. Figure 11b shows the flow of the thrust bearing under different gas supply pressures. The results show that after considering the influence of FSI, the flow rate of the bearing increases significantly, and the deviation increases with the increase in the gas supply pressure. In this case, the thickening of the thrust bearing's gas film is caused by the bending of the thrust plate under the gas film force, thereby increasing the flow rate of the thrust bearing. At the same time, Figure 11 also shows changes in the stiffness and flow of thrust bearings with different grid densities. It can be seen that as the grid density increases, the calculated bearing stiffness and flow will change slightly, but the effect on the overall results is not big. Therefore, the grid density has little influence on the algorithm, which means that a coarser grid can be selected to reduce the calculation time during theoretical analysis.
The comparison of the experimental stiffness of the thrust bearing with the theoretical stiffness considering FSI is given in Figure 12a. It can be seen that the error between the theoretical stiffness and the experimental value is less than 7%, and the calculation accuracy is greatly improved compared with the calculation results without considering FSI. At the same time, with an increase in gas supply pressure, the changing trend of the theoretical stiffness is similar to the experimental value. As shown in Figure 12b, the theoretical analysis and experimental results reveal that the upper thrust plate experienced its maximum deformation. It can be seen that the theoretical calculation error is only 4.3%  The results show that after considering the influence of FSI, the flow rate of the bearing increases significantly, and the deviation increases with the increase in the gas supply pressure. In this case, the thickening of the thrust bearing's gas film is caused by the bending of the thrust plate under the gas film force, thereby increasing the flow rate of the thrust bearing. At the same time, Figure 11 also shows changes in the stiffness and flow of thrust bearings with different grid densities. It can be seen that as the grid density increases, the calculated bearing stiffness and flow will change slightly, but the effect on the overall results is not big. Therefore, the grid density has little influence on the algorithm, which means that a coarser grid can be selected to reduce the calculation time during theoretical analysis.
The comparison of the experimental stiffness of the thrust bearing with the theoretical stiffness considering FSI is given in Figure 12a. It can be seen that the error between the theoretical stiffness and the experimental value is less than 7%, and the calculation accuracy is greatly improved compared with the calculation results without considering FSI. At the same time, with an increase in gas supply pressure, the changing trend of the theoretical stiffness is similar to the experimental value. As shown in Figure 12b, the theoretical analysis and experimental results reveal that the upper thrust plate experienced its maximum deformation. It can be seen that the theoretical calculation error is only 4.3% compared with the experimental value. Figure 12c illustrates the theoretical and experimental flow rates of thrust bearings under varying supply pressures. By comparing the results in Figure 11b with the simulations without the FSI effect, it is evident that the theoretical model derived in this paper can reduce the flow calculation error from 35.3% to 10.7%.
The above results prove that the FSI model established in this paper can more accurately analyze the aerostatic bearing's static performance. On the other hand, it also explains why there is a large error between the calculation results of the traditional theoretical model and the experimental values. However, although the calculation error is greatly reduced, in actual engineering, it is still affected by non-ideal factors, such as the surface roughness of the gas film and processing and assembly errors. compared with the experimental value. Figure 12c illustrates the theoretical and experimental flow rates of thrust bearings under varying supply pressures. By comparing the results in Figure 11b with the simulations without the FSI effect, it is evident that the theoretical model derived in this paper can reduce the flow calculation error from 35.3% to 10.7%. The above results prove that the FSI model established in this paper can more accurately analyze the aerostatic bearing's static performance. On the other hand, it also explains why there is a large error between the calculation results of the traditional theoretical model and the experimental values. However, although the calculation error is greatly reduced, in actual engineering, it is still affected by non-ideal factors, such as the surface roughness of the gas film and processing and assembly errors.

Structural Size Optimization Design
Since thrust bearings have a higher load capacity and stiffness index than journal bearings, the non-negligible deformation of the thrust plate is a leading factor in the FSI problem. According to the aforementioned FSI simulation results and experimental verification, it can be seen that the thrust plate is affected by FSI and undergoes bending deformation, which leads to the gas film thickness difference between the design value and the actual value, thereby changing the performance of the thrust bearing. Therefore, structural deformation is coupled with the flow field in the gas film.
The main factors affecting the deformation of the spindle structure are the structural stiffness and gas film forces, which are responsible for determining the thickness of gas films. Among them, the aerostatic bearing is designed based on the gas film force and cannot be changed arbitrarily. Hence, as a means of ensuring the accuracy of the theoretically designed optimal gas film thickness in engineering design, the structural stiffness of the spindle must be increased. Figure 13 shows the critical dimensional factors affecting the spindle stiffness of a typical I-shaped structure. The thrust plate's thickness H, the thrust plate's outer diameter D1 and the outer diameter D2 of the mandrel directly determine the overall stiffness of the shaft. FSI simulation was carried out by using the method for orthogonal design, in which three factors correspond to four levels, as shown in Table  1, with the purpose of analyzing the effects of the above parameters on the static performance of the aerostatic bearing. The L16 (43) orthogonal table was chosen to carry out the FSI simulation, and the thrust stiffness K (N/μm) of the aerostatic bearing was used as the target. Orthogonal design parameters and their FSI simulation results are listed in Table  2. The gas supply pressure was set to 5.5 Bar.

Structural Size Optimization Design
Since thrust bearings have a higher load capacity and stiffness index than journal bearings, the non-negligible deformation of the thrust plate is a leading factor in the FSI problem. According to the aforementioned FSI simulation results and experimental verification, it can be seen that the thrust plate is affected by FSI and undergoes bending deformation, which leads to the gas film thickness difference between the design value and the actual value, thereby changing the performance of the thrust bearing. Therefore, structural deformation is coupled with the flow field in the gas film.
The main factors affecting the deformation of the spindle structure are the structural stiffness and gas film forces, which are responsible for determining the thickness of gas films. Among them, the aerostatic bearing is designed based on the gas film force and cannot be changed arbitrarily. Hence, as a means of ensuring the accuracy of the theoretically designed optimal gas film thickness in engineering design, the structural stiffness of the spindle must be increased. Figure 13 shows the critical dimensional factors affecting the spindle stiffness of a typical I-shaped structure. The thrust plate's thickness H, the thrust plate's outer diameter D 1 and the outer diameter D 2 of the mandrel directly determine the overall stiffness of the shaft. FSI simulation was carried out by using the method for orthogonal design, in which three factors correspond to four levels, as shown in Table 1, with the purpose of analyzing the effects of the above parameters on the static performance of the aerostatic bearing. The L16 (43) orthogonal table was chosen to carry out the FSI simulation, and the thrust stiffness K (N/µm) of the aerostatic bearing was used as the target. Orthogonal design parameters and their FSI simulation results are listed in Table 2. The gas supply pressure was set to 5.5 Bar.
The average values of the simulation results in Table 2 were calculated according to the method for same-level range analysis (R), and finally, among the above three structural factors, the performance of the thrust bearing is most affected by the outer diameter of the thrust plate D 1 , whereas the thrust plate's thickness H has the greatest impact on the performance of the thrust bearing. The thrust bearing's performance is minimally affected. In Figure 14, the relationship between various influencing factors and stiffness is given. Figure 14a shows that the influence of the thrust plate's thickness H on stiffness increases linearly when it is less than 40 mm, and the influence gradually decreases when it is greater than 40 mm. As the thrust plate's thickness increases, its structural rigidity increases, and the deformation of the thrust plate caused by the gas film pressure decreases. Figure 14b shows that the stiffness of the thrust bearing increases linearly as the outer diameter D 1 of the thrust plate increases, because the areas with high pressure of the bearing increase as the outer diameter of the thrust bearing increases. At the same time, as shown in Figure 14c, as the outer diameter D 2 of the shaft increases, the areas with high pressure on the thrust bearing decrease, resulting in a decrease in the bearing stiffness.    be comprehensively designed according to conditions such as performance and control accuracy requirements.
given. Figure 14a shows that the influence of the thrust plate's thickness H on stiffness increases linearly when it is less than 40mm, and the influence gradually decreases when it is greater than 40 mm. As the thrust plate's thickness increases, its structural rigidity increases, and the deformation of the thrust plate caused by the gas film pressure decreases. Figure 14b shows that the stiffness of the thrust bearing increases linearly as the outer diameter D1 of the thrust plate increases, because the areas with high pressure of the bearing increase as the outer diameter of the thrust bearing increases. At the same time, as shown in Figure 14c, as the outer diameter D2 of the shaft increases, the areas with high pressure on the thrust bearing decrease, resulting in a decrease in the bearing stiffness.

Conclusions
In order to analyze and reveal the influence of fluid-solid coupling on the static gas pressure spindle as a result of the gas film force and its influence on the static performance of the aerostatic bearing, a fluid-solid coupling calculation model was constructed in this study. For the purpose of improving the accuracy of the calculations and the convergence speed, DCM is introduced to solve the gas film flow field and structure field equations, and the numerical solution quickly converges. According to the proposed FSI-based model, the influence of structure size and the effect of supply pressure on the aerostatic bearing's performance were investigated.
The main results achieved are drawn in the following manner: 1.
The FSI model established in this work based on the DCM method can greatly improve the calculation accuracy of the aerostatic bearing's performance. The model uses the finite element method to simultaneously calculate the solid field and the fluid field. Compared with the traditional theoretical model, the calculation error is greatly reduced, and at the same time, the method is less affected by the grid division, which greatly reduces the calculation speed. Therefore, compared with the separation approximation method, greater levels of accuracy and efficiency are achieved in the calculation.

2.
To establish the FSI model of the I-shaped gas static pressure bearing, the DCM method is proposed in this paper, and how the FSI affects the performance of the thrust bearing under different supply pressures is analyzed based on COMSOL, including deformations in the thrust plate and the stiffness and flow rate of the thrust bearing. Experimental evidence confirms the theoretical analysis, and the results show that the FSI model calculates values that are closer to experimental measurements. Compared with simulation results without taking FSI into account, the variation trend of the stiffness of the thrust bearing in the simulation analysis considering the FSI effect is closer to the experimental results with supply pressure increases.

3.
Based on the FSI model, this paper further analyzes the influence of the critical structural dimensions of the I-shaped spindle in terms of the thrust bearing's static performance. According to the results, the thrust plate's thickness has the greatest influence on the thrust stiffness, and the aerostatic bearing design can be optimized on the basis of these results.
In this paper, the FSI model mainly focuses on using the I-shaped orifice-type aerostatic bearing, and in order to optimize the FEM model details, it is necessary to consider the spindle's sensitive positions. At the same time, this method can be extended to other throttling forms of aerostatic bearings, such as porous aerostatic bearings with more obvious effects of fluid-solid coupling problems.
Author Contributions: Conceptualization, Y.W., W.C. and B.W.; methodology, W.C. and Z.Q.; software, Y.W. and W.C.; validation, Y.W.; formal analysis, W.C.; investigation, Y.W. and W.C.; resources, Y.W. and Z.Q.; data curation, Y.W. and W.C.; writing-original draft preparation, W.C.; writingreview and editing, Q.Z.; visualization, Q.Z. and W.C.; supervision, Z.Q. and B.W.; project administration, Z.Q. and B.W.; funding acquisition, Z.Q. and B.W. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The data supporting the results reported by the authors can be requested by e-mail.

Conflicts of Interest:
The authors declare no conflict of interest.