# **Gas Bearings Modelling, Design and Applications**

Edited by Terenziano Raparelli, Federico Colombo, Andrea Trivella and Luigi Lentini

Printed Edition of the Special Issue Published in *Applied Sciences*

www.mdpi.com/journal/applsci

## **Gas Bearings: Modelling, Design and Applications**

## **Gas Bearings: Modelling, Design and Applications**

Editors

**Terenziano Raparelli Federico Colombo Andrea Trivella Luigi Lentini**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Terenziano Raparelli Politecnico di Torino Italy Luigi Lentini Politecnico di Torino Italy

Federico Colombo Politecnico di Torino Italy

Andrea Trivella Politecnico di Torino Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/Gas Bearings).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5607-9 (Hbk) ISBN 978-3-0365-5608-6 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## *Editorial* **Special Issue "Gas Bearings: Modelling, Design and Applications"**

**Federico Colombo \*, Luigi Lentini, Terenziano Raparelli and Andrea Trivella**

Politecnico di Torino, Department of Mechanical and Aerospace Engineering, Corso Duca degli Abruzzi 24, 10129 Turin, Italy

**\*** Correspondence: federico.colombo@polito.it

#### **1. Introduction**

Gas bearings are widely employed in high-precision devices and in high-speed applications, such as in micro turbomachinery and micro machining tools. In metrology and high-tech industry, the absence of contact between the surfaces that are in relative motion is appreciated due to the possibility of obtaining high-precision motion. Recent trends in small-scaling turbomachinery demand increasingly higher rotational speeds; the advantage is the ability to operate with higher power to weight ratios. Additionally, in this case, gas bearings represent the only solution due to their simplicity and long operating life, with no need for maintenance.

This issue collects research and review papers covering the modelling, design, and application of gas bearings. Each published paper distinguishes itself from the others in the sense that it deals with the problem from an original point of view. Some papers analyze several prototypes to improve their performance, while others concentrate more on the mathematical modelling of bearings, developing analytical or numerical formulations with lumped or distributed parameters. These different approaches show how complex the analysis of gas bearings can be, but they also show how fascinating the process of finding simple solutions to complex problems is.

#### **2. Outline of the Issue**

The papers published in the Special Issue are here summarized in the order of their publication dates. Most of them investigate gas bearings using theoretical and experimental approaches. Numerical models, after they have been experimentally validated, are important tools for the design and optimization of gas bearings. Current models are able to consider multi-physical problems such as fluid–structure interaction, thermal effects, and porous materials. In some cases, tailor-made software is written to solve the equations, and in other cases, commercial software is employed. The experimental results confirm the validity and accuracy of multi-physical numerical models.

In high-precision devices, such as fly cutting machine tools or rotary tables, the deformation of solid parts is often not negligible, as it influences the performance of the aerostatic bearings. In [1], the fluid–structure interaction effect is considered to investigate the static and dynamic performance of an aerostatic spindle. In particular, the thickness of the thrust plate is optimized to reach a compromise between the reduction of the spindle mass and the increase in the thrust bearing stiffness.

Foil gas bearings are particularly useful in turbomachinery, as they are able to sustain high-speed blowers, compressors, and turbochargers without the requirement of an external air supply. Paper [2] describes the thermo-hydrodynamic model of a bump foil thrust gas bearing. The cooling models of the foil structure and thrust plate as well as the influence of temperature on the thermal expansion of the bearing structure are considered. The effects of the bearing speed, thrust load, and external cooling gas on the bearing temperature field are calculated and analyzed.

**Citation:** Colombo, F.; Lentini, L.; Raparelli, T.; Trivella, A. Special Issue "Gas Bearings: Modelling, Design and Applications". *Appl. Sci.* **2022**, *12*, 9048. https://doi.org/10.3390/ app12189048

Received: 8 September 2022 Accepted: 8 September 2022 Published: 8 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Paper [3] reviews the state-of-the-art of gas foil bearings in China, many of which have been developed for use in high-speed turbo machinery. Different types of foil bearings are described (bump, multi-leaf, wire, viscoelastic, spring and protuberant), as are their different applications (e.g., fuel cell air compressors, blowers, turbo-expanders, and oil-free turbochargers). The addressed challenges concern thermal management, rotor-dynamic stability, and wear-resistant coatings.

High-tech industry products such as solar panels, chips, and displays are composed of very thin layers of brittle materials that need to be handled with care. An original solution to handle the thin substrates without mechanical contact is proposed in [4]. The system works by floating the substrate on a thin film of air and by creating a viscous traction force on it. The actuator is improved to reach a bandwidth of 300 Hz, while the position control loop of the substrate reaches a bandwidth of 10 Hz with a position error of less than 13 μm. It is a good example of the so-called tribotronics discipline.

It was shown in the literature that the radial error motion of an aerostatic journal bearing is mainly influenced by the shape errors of the shaft. The influence of the shaft roundness and the cylindricity errors on the aerostatic spindle's rotation accuracy is evaluated in [5]. The research proposes an index that is more effective in predicting the spindle's rotation accuracy than the roundness and cylindricity.

Paper [6] presents an elegant analytical model of a porous gas journal bearing after a literature survey on the topic. The analytical model is able to predict the flow and the dynamic force coefficients. An FE model is compared to the analytical one, and a perfect match is found. A validation of the static results with experimental data from the literature is also carried out. The paper represents a useful design guide for the quick selection of physical parameters, both for static performance and for stability analysis.

A partial arc annular-thrust porous journal bearing that can carry both radial and axial loads is investigated in [7]. The analysis is carried out using commercial multi-physical software in order to evaluate the pressure distribution for low eccentricity as well as for low tilting angles and tilting speeds. The dynamic coefficients of tilting motion are then determined using the finite difference technique.

Finally, the modal parameters of an aerostatic spindle are identified experimentally based on the impulse response in paper [8]. Different methods are proposed to identify the damped natural frequencies and the damping ratios of the conical and cylindrical modes of the spindle. The same modal parameters are also evaluated with a non-linear numerical model, and a comparison with experimental results is provided to validate the model.

#### **3. Future of Gas Bearings**

The different gas bearing solutions require accurate mathematical models and cannot be easily standardized due to the complexity and variety of the problems. Despite gas lubrication being a well-established discipline since the 1970s, the ever-rising demands of modern technology offer new challenges to obtain higher precision and higher rotational speeds, stimulating researchers to broaden the frontiers of gas lubrication.

**Author Contributions:** Writing—original draft preparation, F.C.; Writing—review and editing F.C., L.L., T.R. and A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The editors would like to acknowledge all of the authors for their high-quality contributions to this Special Issue. The MDPI management and staff are also acknowledged for their continuous support.

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

#### **References**


## *Article* **A Two-Round Optimization Design Method for Aerostatic Spindles Considering the Fluid–Structure Interaction Effect**

**Qiang Gao 1,2, Siyu Gao 1,\*, Lihua Lu 1,\*, Min Zhu <sup>3</sup> and Feihu Zhang <sup>1</sup>**


**Abstract:** The fluid–structure interaction (FSI) effect has a significant impact on the static and dynamic performance of aerostatic spindles, which should be fully considered when developing a new product. To enhance the overall performance of aerostatic spindles, a two-round optimization design method for aerostatic spindles considering the FSI effect is proposed in this article. An aerostatic spindle is optimized to elaborate the design procedure of the proposed method. In the first-round design, the geometrical parameters of the aerostatic bearing were optimized to improve its stiffness. Then, the key structural dimension of the aerostatic spindle is optimized in the second-round design to improve the natural frequency of the spindle. Finally, optimal design parameters are acquired and experimentally verified. This research guides the optimal design of aerostatic spindles considering the FSI effect.

**Keywords:** aerostatic spindle; fluid–structure interaction; restrictor geometrical parameters; structural dimension; optimal design

#### **1. Introduction**

Aerostatic spindles are widely employed in various kinds of precision equipment for their distinct advantages of ultra-high precision, negligible friction, and wide rotating speed range [1]. However, the drawbacks of low capacity, low stiffness, and poor high-speed stability restrict their possible application significantly. Many factors have a great impact on the overall performances of aerostatic spindles, such as the geometrical parameters of the aerostatic bearing, the supply pressure, the rotating speed, and the restrictor type, etc. Hence, it is difficult to optimize the design parameters of aerostatic spindles comprehensively.

Numerous research studies have been conducted to facilitate the design of aerostatic bearings by investigating the impact of the air film geometrical dimensions on their static performance [2,3]. Li et al. [4] numerically studied the performances of the orifice-type restricted aerostatic thrust bearings with varying geometrical parameters based on computation fluid dynamics, and the effect tendency of geometrical parameters to the bearing performances was acquired. Belforte [5] and Du et al. [6] investigated the stiffness, load-carrying capacity (LCC) of aerostatic thrust bearings and journal bearings with pressure-equalizing groove (PEG), and found that the PEG can enhance the static performance of aerostatic bearing dramatically.

Many studies also have been carried out on the optimal design approach of air bearings [7–9]. Wang et al. [10] proposed an optimal design approach based on a genetic algorithm, and a hybrid selection scheme was adopted when selecting the mating groups, which is more efficient. Federico et al. [11,12] optimized the orifice diameter, orifice position, and orifice number of inherently compensated rectangular aerostatic pad bearings based

**Citation:** Gao, Q.; Gao, S.; Lu, L.; Zhu, M.; Zhang, F. A Two-Round Optimization Design Method for Aerostatic Spindles Considering the Fluid–Structure Interaction Effect. *Appl. Sci.* **2021**, *11*, 3017. https:// doi.org/10.3390/app11073017

Academic Editor: Terenziano Raparelli

Received: 12 March 2021 Accepted: 25 March 2021 Published: 28 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on a multi-objective genetic algorithms optimization approach. To optimize the design parameters and supply pressure of a rectangular aerostatic thrust bearing, Nikhil et al. [13] proposed a Pareto optimization method to avoid the variation of the optimal results with different initial populations; the trial points are therefore distributed in the design space uniformly.

However, the structural deformation caused by pressurized air film was neglected in the aforementioned researches. Since the performance of the aerostatic bearings is significantly impacted by the air film thickness, this phenomenon may induce great error when calculating the performance of the aerostatic, and the initial design requirements may not be satisfied by the actual performance of the aerostatic bearing. This phenomenon is especially dramatic in the design of aerostatic bearings with high stiffness and high LCC, such as the spindle of fly cutting machine tools, rotary tables, etc., and has drawn increasing attention in recent years [14,15]. Lu et al. [16,17] numerically and experimentally investigated the performance of an air spindle considering the fluid–structure interaction (FSI) effect. They found that the LCC curve declines slightly due to the FSI effect, and the axial stiffness declines with the decrease of the thrust thickness plate. In a study by Gao et al. [18], the performance of an aerostatic spindle was investigated by a two-way FSI simulation model. It showed that the air film thickness changes about 7.65 μm due to the structure deformation and the stiffness of the spindle declines about 34%.

According to the above literature review, it can be seen that the design of aerostatic bearing includes the design of the air bearing parameters and the design of the crucial dimensions of the solid parts. Even though many optimal design approaches of air bearings have been proposed, the FSI effect is seldom considered, and the crucial structural dimensions are usually decided empirically. Therefore, the aerostatic spindle cannot be optimized comprehensively from its LCC, volume flow rate (VFR), stiffness, natural frequency, etc.

To overcome this issue, a two-round optimization design method for the aerostatic spindle with consideration of the FSI effect is proposed in the current study. The design parameters for the aerostatic bearing and the solid parts can be optimized by the proposed approach. As an application case, an aerostatic spindle is optimized. The impact of the crucial parameters of air film and the structural dimensions on the performances of the aerostatic spindle are investigated based on the finite element method (FEM). The performance of the aerostatic spindle is enhanced from the aspects of LCC, VFR, stiffness, natural frequency, etc. Furthermore, experiments were implemented to validate the reliability of the calculation results.

#### **2. A Two-Round Optimization Design Method of Aerostatic Spindles Considering the FSI Effect**

According to the existing design philosophy of aerostatic bearings [19–21], the design of air bearing is usually first carried out based on its performance requirements and design restrictions, and then the crucial solid parts are designed empirically based on the dimensions of air film. The elastic deformation of crucial parts caused by high air pressure and its impact on the overall performance of the aerostatic spindle are ignored. However, when the FSI effect is considered, the design of the aerostatic spindle consists of the design of the geometrical parameters of air bearing and the design of key structural dimensions. Both aspects impact the overall performance of the spindle system significantly. To facilitate the optimal design of aerostatic spindles considering FSI, a two-round optimization design method is proposed in this research, as shown in Figure 1. Its basic idea is to optimize the geometrical parameters of aerostatic bearing and the key structural dimensions of solid parts sequentially so that the static and dynamic performance of the spindle system can be optimized comprehensively.

In the first-round design stage, the static performances of the aerostatic bearing are improved by optimizing its geometrical parameters. The FEM model of aerostatic bearing is built in this stage. Then, the influence of the crucial design parameters (such as the orifice diameter, the air film thickness, and the orifice number) on the static performance of the aerostatic bearing can be acquired. Finally, the optimized design parameters for air bearing

design can be acquired by selecting a good combination of LCC, stiffness, and volume flow rate (VFR).

**Figure 1.** A two-round optimization design method for aerostatic spindles considering the fluid–structure interaction (FSI).

In the second-round design stage, the static and dynamic performances of the aerostatic spindle are further investigated with consideration of the FSI effect. A two-way FSI simulation model is built to calculate the structural deformation of the spindle and its impact on the static performance of the aerostatic bearing. The crucial structural dimensions that greatly influence the structural rigidity and weight of solid components are analyzed. Furthermore, to clarify the effect of structure dimension on the natural frequency of the spindle, modal analysis with varying structure dimensions is carried out, and finally, the optimal structural dimensions are acquired.

In the following sections, the optimal design of an aerostatic spindle is implemented as a case study to elaborate the detailed procedures of the proposed two-round optimization design method.

#### **3. The Configuration of an Aerostatic Spindle**

Figure 2a presents the configuration of an aerostatic spindle that is employed in an ultra-precision fly cutting machine tool. Its rotating parts consist of a shaft and two thrust plates. Two air thrust bearings and a radial bearing are adopted to separate the rotating parts from the shaft sleeve and resist the external load. In ultra-precision fly cutting machining, the axial direction of the spindle is the error-sensitive direction. Therefore, the axial vibration and tilt vibration of the rotating parts has a significant impact on the machined surface quality [22]. It indicates that enhancing its axial stiffness and angular stiffness is beneficial to reduce its vibration amplitude and the machined surface quality.

Figure 2b demonstrates the geometrical parameters of the aerostatic thrust bearing. The restrictor of aerostatic thrust bearings is a pocketed orifice-type restrictor. The inner diameter of the thrust bearing is *D*1, and the outside diameter of the thrust bearing is *D*3. There are *n* orifices distributes equally along the circumferential direction of the thrust bearing, and the diameter is *D*2. *d* and *h* are the diameter of the orifice restrictor and the thickness of air film, respectively. *d*<sup>1</sup> and *h*<sup>1</sup> are the diameter and the depth of the recess, respectively. Table 1 shows the initial values of these parameters.

**Figure 2.** The configuration of an aerostatic spindle and the design parameters of its thrust bearing. (**a**) The configuration of the aerostatic spindle and (**b**) the design parameters of the thrust bearing.

**Table 1.** The design parameters of the aerostatic thrust bearing.


#### **4. The First-Round Optimal Design of Aerostatic Thrust Bearings**

The static performances of an aerostatic bearing are greatly impacted by the geometrical parameters of the orifice-type restrictor, such as the orifice diameter, the air film thickness, and the orifice number. Therefore, its performances can be improved by optimizing these parameters. In this section, a numerical model of the aerostatic thrust bearing is built first based on the FEM. Additionally, then, the impact of the above three parameters on the static performances of aerostatic thrust bearing is calculated based on the FEM model. Finally, the static performances of the aerostatic thrust bearing are improved by optimizing these three design parameters.

#### *4.1. FEM Modeling of Aerostatic Bearings*

The flow status of an orifice restrictor is simplified as the compressible flow through an ideal nozzle, and the orifice is simplified as a point. When the air flows through the orifice, the air pressure drops from *Ps* to *Pd* immediately. Considering the difference between the actual mass flow rate of the orifice restrictor and the ideal theoretical mass flow rate, a constant called discharge coefficient *Cd* = 0.8 is introduced to correct the mass flow rate. The mathematical models for calculating the theoretical mass flow rate of the orifice restrictor are Equations (1) and (2) [20].

$$m = \mathbb{C}\_d A \, p\_s \sqrt{\frac{2\rho\_a}{p\_a}} \psi\_{s\prime} \tag{1}$$

$$\psi\_{s} = \left\{ \begin{array}{c} \left[ \frac{k}{2} \left( \frac{2}{k+1} \right)^{(k+1)/(k-1)} \right]^{1/2}; \left( \frac{p\_{d}}{p\_{s}} \le \beta\_{k} \right) \\\left\{ \begin{array}{c} \frac{k}{k-1} \left[ \left( \frac{p\_{d}}{p\_{s}} \right)^{2/k} - \left( \frac{p\_{d}}{p\_{s}} \right)^{(k+1)/k} \right] \\\end{array} \right\}; \left( \frac{p\_{d}}{p\_{s}} > \beta\_{k} \right) \end{array} \right\}; \left( \beta\_{k} = \left( 2/(k+1) \right)^{(k+1)/k} \right) \tag{2}$$

where *A* = *πd*2*/4* is the cross-section area of orifice, *Pa* is the atmosphere pressure, *ρ<sup>a</sup>* is the density of air under standard state, *ψ<sup>s</sup>* is flow rate function, and *k* is the specific heat ratio of air.

The fluid domain of the aerostatic thrust bearing is described using the governing equation in the cylindrical coordinate system [20] below:

$$
\frac{\partial}{\partial \overline{\xi}} \left( \overline{h^3} \frac{\partial \overline{p^2}}{\partial \overline{\xi}} \right) + \frac{\partial}{\partial \theta} \left( \overline{h^3} \frac{\partial \overline{p^2}}{\partial \theta} \right) + \overline{r^2} \overline{\mathbb{Q}} \delta\_i = \Lambda \overline{r} \frac{\partial}{\partial \theta} \left( \overline{l} \overline{p} \right), \tag{3}
$$

where *h* is the dimensionless thickness of air film, *p* is the dimensionless pressure, *ξ* is the dimensionless coordinate in the radial direction, *r* is the coordinate in the radial direction, and *θ* is the coordinate in the circumferential direction. *Q* is the mass flow rate factor of the orifice restrictor. *δ<sup>i</sup>* is a constant coefficient with a value of 1 or 0 (in the computational domain where there is an orifice, it equals 1; otherwise, it equals 0.) Λ is the dimensionless bearing number.

$$
\overline{h} = \frac{h}{h\_{\text{ff}}}, \overline{\tau} = \frac{r}{r\_0}, \overline{\xi} = \ln \overline{\tau}, \overline{p} = \frac{p}{p\_s}, \overline{\mathbb{Q}} = \frac{24 \eta r\_0^2 p\_a}{h\_m^3 p\_s^2 p\_a} \rho \overline{v}, \Lambda = \frac{12 \eta v \rho\_{\overline{\theta}2} r\_0}{h\_m^2 p\_s}, \tag{4}
$$

where *h* is the thickness of air film. *r0* is the characteristic length, *hm* is the reference air film thickness, *p* is the gas pressure distribution function, *ρ* is the air density, *η* is air dynamic viscosity, <sup>∼</sup> *v* is the velocity of orifice gas flow, and *vθ*<sup>2</sup> is the velocity components in the circumferential direction. In this paper, the performance of aerostatic thrust bearing is analyzed under motionless conditions, and therefore, Λ = 0.

The pressure square function is defined as follows:

$$
\overline{f} = \overline{p^2}(\overline{\xi}, \theta) = \overline{f}(\overline{\xi}, \theta). \tag{5}
$$

Taking the pressure square function *f* as the variable, the functional Φ can be constructed based on Equation (3) as follows:

$$\Phi\left(\overline{f}\left(\overline{\xi},\theta\right)\right) = \int\_{\Omega} \left\{ \frac{\overline{h^3}}{2} \left[ \left( \frac{\overline{\partial f}\left(\overline{\xi},\theta\right)}{\partial\overline{\xi}^\dagger} \right)^2 + \left( \frac{\overline{\partial f}\left(\overline{\xi},\theta\right)}{\partial\theta} \right)^2 \right] - \overline{Qf}\left(\overline{\xi},\theta\right)\delta\_i \right\} d\overline{\xi}d\theta \tag{6}$$

where Ω denotes the computational domain of the air film. It has been mathematically proved that, under certain boundary conditions, the extremum of Equation (6) can only be acquired by a certain pressure square function *f* that is the solution of Equation (3).

The triangular finite element is employed to solve Equation (6) in this paper. *i*, *j*, *m* are the three nodes of triangular finite elements. The interpolation function of triangular finite element is defined as

$$f = \mathbf{N}^{\mathbf{c}^T} \mathbf{f}^{\mathbf{c}} \tag{7}$$

$$\mathbf{f}^{\epsilon} = \begin{bmatrix} f\_i & f\_j & f\_{\mathbf{w}\_i} \end{bmatrix}^T \tag{8}$$

$$\mathbf{N}^{\mathbf{c}} = \begin{bmatrix} N\_i & N\_j & N\_m \end{bmatrix}^T \tag{9}$$

where **f**<sup>e</sup> represents the pressure square of the nodes. **N**<sup>e</sup> is the shape function of the triangular element, which can be calculated by Equation (10).

$$\mathbf{N}\_{l} = \frac{1}{2\Delta e} (a\_{l} + b\_{l}\theta + c\_{l}\overline{\xi})\_{\prime} (i = i, j, m) \tag{10}$$

$$\begin{cases} a\_i = \theta\_j \overline{\xi\_{\mathbb{S}^m}} - \theta\_m \theta\_j \\\ b\_i = \overline{\xi\_j} - \overline{\xi\_m} \\\ c\_i = \theta\_m - \theta\_j \end{cases} \tag{11}$$

where Δ*e* is the area of the triangular element. By changing the subscripts of Equation (11) in the order of *i*, *j*, *m* sequentially, and *Nj* and *Nm* can be obtained.

By substituting Equation (7) into Equation (6) and applying it to the finite elements of the entire computational domain, Equation (12) can be derived as follows:

$$\begin{split} \Phi(f) &= \sum\_{\epsilon=1}^{m} \frac{1}{2} \int\_{\Lambda \mathcal{C}} \overline{h^3} \left[ \left[ \frac{\partial}{\partial \overline{\mathbf{x}}} \left( \mathbf{N}^{\epsilon T} \mathbf{f}^{\epsilon} \right) \right]^2 + \left[ \frac{\partial}{\partial \overline{\mathbf{z}}} \left( \mathbf{N}^{\epsilon T} \mathbf{f}^{\epsilon} \right) \right]^2 \right] d\_5^{\overline{\mathbf{x}}} d\theta \\ &\quad - \sum\_{\epsilon=1}^{m} \int\_{\Lambda \mathcal{C}} \overline{\mathbf{Q}} f d\_5^{\overline{\mathbf{x}}} d\theta \delta\_i \end{split} \tag{12}$$

The partial derivative with respect to *f* should be zero to obtain the solution of Equation (6); thus, we can obtain Equation (13). The solution of Equation (13) is the pressure square function *f* , which makes the functional Φ reach the extremum.

$$\frac{\partial \Phi}{\partial f\_i} = \sum\_{\mathbf{c} \in \Lambda\_i} \left( c\_i \mathbf{c}^{\epsilon T} + b\_i \mathbf{b}^{\epsilon T} \right) \mathbf{f}^{\epsilon} \int\_{\Lambda \epsilon} \overline{h^3} d\overline{\xi} d\theta / \left( 2\Delta \epsilon \right)^2 - k\_1 \mu\_r \dot{m}\_r \delta\_i = 0 \tag{13}$$
 
$$(i = 1, 2, \dots, n)$$

Additionally, Equation (13) can be rewritten in the matrix form as follows:

$$\mathbf{K}\mathbf{F} = \mathbf{T},\tag{14}$$

where **F** = [*f1 f2 fn*] <sup>T</sup> is a column vectors with *<sup>n</sup>*×1 dimensions. **<sup>T</sup>** is composed of three types of elements: 0, *K*<sup>1</sup> . *mrμ<sup>r</sup>* and − *c* (*k*) *<sup>i</sup> c* (*k*) *<sup>j</sup>* + *b* (*k*) *<sup>i</sup> b* (*k*) *<sup>j</sup>* <sup>×</sup> <sup>Δ</sup>*<sup>k</sup> h*3*dξdθ*/(2Δ*k*) <sup>2</sup> . . *mr* is the real mass flow rate of the orifice restrictor. **K** is a square matrix with *n*×*n* dimensions, and its elements can be acquired by

$$K\_{ij} = \sum\_{c \in \Delta\_i \land \Delta\_j} (c\_i c\_j + b\_i b\_j) \int \overline{h^3} d\overline{\xi} d\theta (2\Delta \varepsilon)^2. \tag{15}$$

Then, the LCC, the stiffness, and the VFR of the aerostatic thrust bearing can be calculated by the equations below. The LCC of the aerostatic thrust bearing is the sum of the integration of the pressure on each element. Then, the stiffness of the aerostatic thrust bearing can be calculated by Equation (17). The VFR can be acquired by Equation (18).

$$\mathcal{W}\_{\rm LCC} = p\_s r\_0^2 \sum\_{c=1}^{m} \int\_{\Delta \mathbf{z}} \overline{p} d\_s^{\overline{\mathbf{z}}} d\theta \tag{16}$$

$$K\_{\rm stiffges} = \frac{W\_{\rm LCC}(h + \Delta h) - W\_{\rm LCC}(h)}{\Delta h} \tag{17}$$

$$Q\_{\rm VFR} = m/\rho.\tag{18}$$

Figure 3 presents the FEM model of the aerostatic thrust bearing. Since the thrust air film is symmetric in the circumferential direction, a basic sector is built and the periodic boundary conditions are employed to save calculation time. The computational domain is equally divided into *n<sup>θ</sup>* parts along the circumference direction, while it is divided into *n<sup>ξ</sup>* parts along the radial direction.

#### *4.2. Optimal Design of the Geometrical Parameters of the Restrictor*

To optimize the orifice diameter, the air film thickness, and the orifice number of the thrust bearing, their impact on the performance of aerostatic thrust bearing is investigated in this section. Except for these three variables, the values of other geometrical parameters of the air bearing are the same as shown in Table 1.

Figure 4 presents the performance curves of aerostatic thrust bearings with varying orifice diameter *d* and air film thickness *h*. It can be seen that with the increase of the air film thickness, the stiffness curve ascents first and then descends, which has a peak. However, the LCC curves decline, and the VFR curve increases monotonically. Moreover, with the increase of the orifice diameter, the LCC and the VFR of the aerostatic thrust bearing increase dramatically. However, for the stiffness curves, the maximum stiffness declines greatly with the increase of orifice diameter. Moreover, the air film thickness, which is corresponding to the peak position of stiffness curves, declines gradually along with the decline of orifice diameter. It indicates that a higher stiffness can be acquired by the combination of smaller orifice diameter and a smaller air film thickness.

**Figure 3.** The finite element method (FEM) model of the air thrust bearing.

**Figure 4.** The performance curves of aerostatic thrust bearing with varying *d* and *h.*

Figure 5 presents the performance curves of the aerostatic thrust bearing with different orifice numbers *n*. It shows that the LCC and stiffness of aerostatic thrust bearing can be improved by increasing the orifice number. However, it can also cause a dramatic increase in VFR. In addition, the LCC and stiffness only rise slightly when the orifice number increases to 12, while the VFR increases greatly. This means that the performance of aerostatic thrust bearing cannot be enhanced continuously with increasing orifice number.

**Figure 5.** The performance curves of aerostatic thrust bearing with varying orifice number *n.*

According to the above investigation, the combination of smaller orifice diameter and a smaller air film thickness helps improve the stiffness, and increasing the orifice number

can also improve the stiffness. To improve the performance of aerostatic thrust bearing, the orifice diameter of 0.15 mm, the air film thickness of 13 μm, and the orifice number of 8 are selected in this section. Compared with the initial design, the LCC is improved by 28.3%, from 781.0 N to 1001.8 N, the stiffness is improved by 44.2%, from 72.7 N/μm to 104.8 N/μm, and the VFR is reduced by 27.2%, from 9.2 L/min to 6.7 L/min.

#### **5. The Second-Round Optimal Design of Crucial Structural Dimensions**

The geometrical parameters of air film determine the design of structural dimensions. After the design parameters of the aerostatic thrust bearing are optimized, the structural dimensions thrust plates are designed empirically based on the dimensions of thrust air film. For example, the outer diameter of air thrust film *D*<sup>3</sup> should be less than the diameter of thrust plates. However, the thrust plate thickness has a dramatic impact on its bending rigidity and it is generally designed empirically. When the FSI effect is considered, the bending deformation of the thrust plate caused by the air film force could increase the air film thickness and change the stiffness of the aerostatic bearing. Generally, improving the thickness of the thrust plate is beneficial to strengthen its bending rigidity. However, it will also improve the weight of rotating parts, which decreases the natural frequency of the spindle. Therefore, the thrust plate thickness that greatly influences the structural rigidity and weight of rotating parts should be further optimized from the perspective of the natural frequency of the spindle.

#### *5.1. FSI Modeling of Aerostatic Thrust Bearing*

To calculate the performance of aerostatic thrust bearing with consideration of elastic deformation of the thrust plate, a two-way FSI modeling method [23] of aerostatic thrust bearing is adopted, as illustrated in Figure 6. It consists of the FEM model of the thrust plate and the FEM model of the aerostatic thrust bearing. The elastic deformation of the thrust plate can be calculated by the FEM model of the thrust bearing. The pressure distribution of the air film can be acquired by the FEM model of the aerostatic thrust film. Since the modeling method of the aerostatic thrust bearing has already been detailed in the previous section, it is not detailed again in this section, and only the modeling method of the thrust plate is explained below.

**Figure 6.** The FSI model of the aerostatic thrust bearing. (**a**) The FEM model of aerostatic thrust bearing; (**b**) the pressure distribution of air film; (**c**) the pressure load applied on thrust plate; (**d**) the FEM model of thrust plate; and (**e**) the structure deformation of the thrust plate.

Considering the structure feature of the thrust plate, a two-dimensional FEM model is built in this section to save calculation time, and the shell element is employed to simulate the bending deformation of the thrust plate. It is usually to simulate the thin-walled to moderately thick-walled structures, which can acquire excellent calculation accuracy and high computational efficiency. As shown in Figure 6d, the inner edge of the thrust plate is constrained by fixed support, and the symmetrical boundary condition is adopted to further save calculation time. The air film pressure is applied on the bearing surface of the thrust plate.

The governing equation for the FEM model of the thrust plate can be written as below.

$$\left[ \left[ K \right] \left\{ \mu \right\} = \left\{ F \right\} \right] \tag{19}$$

where {*u*} and [*K*] represent the displacement vector and stiffness matrix, respectively. {*F*} donates the pressure load applied on the bearing surface of the thrust plate. According to the Mindin–Reissner shell theory, the stiffness matrix of the shell element can be acquired by Equation (20) [24].

$$\mathbf{K}^{\varepsilon} = \frac{h\_{\varepsilon}^{\gamma}}{12} \int \int \mathbf{B}\_{b}^{T} \mathbf{D}\_{b} \mathbf{B}\_{b} dx dy + h\_{\varepsilon} \int \int \mathbf{B}\_{s}^{T} \mathbf{D}\_{s} \mathbf{B}\_{s} dx dy \tag{20}$$

where *he* is the thickness of the shell element. **D** is the elastic matrix as follows:

$$\mathbf{D}\_{b} = \frac{E}{1-\mu^{2}} \begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0 & 0 & (1-\mu)/2 \end{bmatrix}, \mathbf{D}\_{s} = \begin{bmatrix} \kappa G & 0\\ 0 & \kappa G \end{bmatrix} \tag{21}$$

where *E* denotes the elasticity modulus. *μ* is the Poisson ratio. *κ* denotes the shear correction factor. *G* is the shear modulus. **B** can be calculated by Equations (22) and (23).

$$\mathbf{B}\_{b} = \begin{bmatrix} \mathbf{B}\_{b1} & \mathbf{B}\_{b2} & \dots & \mathbf{B}\_{n} \end{bmatrix}, \mathbf{B}\_{s} = \begin{bmatrix} \mathbf{B}\_{s1} & \mathbf{B}\_{s2} & \dots & \mathbf{B}\_{sn} \end{bmatrix} \tag{22}$$

$$\mathbf{B}\_{bi} = \begin{bmatrix} 0 & 0 & \frac{\partial N\_i}{\partial x} \\ 0 & -\frac{\partial N\_i}{\partial y} & 0 \\ 0 & -\frac{\partial N\_i}{\partial x} & \frac{\partial N\_i}{\partial y} \end{bmatrix}, \mathbf{B}\_{si} = \begin{bmatrix} \frac{\partial N\_i}{\partial y} & -N\_i & 0 \\ \frac{\partial N\_i}{\partial x} & 0 & -N\_i \end{bmatrix} \tag{23}$$

where *N* is the shape function of rectangular elements.

$$N\_i = (1 + \zeta\_i \zeta\_i)(1 + \eta\_i \eta)/4(i = 1, 2, 3, 4) \tag{24}$$

where *ξ<sup>i</sup>* and *η<sup>i</sup>* are the coordinates of nodes in the local coordinate system of the shell element.

A crucial part of the two-way FSI modeling is to build the bidirectional data exchange interface between the FEM model of the aerostatic thrust bearing and the FEM model of the thrust plate. To realize it, the data exchange procedure is designed and it can be decomposed into the following steps:

(1) Based on the FEM model of the aerostatic thrust bearing, as shown in Figure 6a, conducting a steady-state analysis of air film and the pressure distribution of air film can be acquired, as shown in Figure 6b;

(2) Exporting the air film pressure data to the two-dimensional FEM model of the thrust plate, as shown in Figure 6d, and applying it as pressure load, as shown in Figure 6c;

(3) Conducting steady-state structure analysis based on the two-dimensional FEM model of the thrust plate, in which the deformation of the thrust plate can be acquired, as shown in Figure 6e;

(4) Transferring the deformation data of the thrust plate to the FEM model of the aerostatic thrust bearing and updating the air film thickness of each element;

The bidirectional data exchange can be accomplished by the above procedures. By implementing the above process repeatedly, the static performance of aerostatic thrust bearing can be acquired.

#### *5.2. Optimal Design of the Structural Dimensions of the Aerostatic Spindle*

Figure 7 presents the axial stiffness and the first-order natural frequency of the aerostatic spindle with varying thrust plate thickness and air film thickness considering the FSI effect. It can be derived from Figure 7a that with the increase of thrust plate thickness, the maximum stiffness increases sharply in the range of 5–15 mm and then increases slowly in the range of 15–30 mm. This phenomenon can be explained as follows. The structural rigidity of the thrust plate is directly related to its thickness. By increasing its thickness, its structural rigidity will be improved, and the deformation of the thrust plate will be smaller. When the thickness of the thrust plate is larger than 15 mm, the stiffness trends to steady gradually with the increase of the thrust plate. It is because that the structural rigidity of the thrust plate is high enough to resist the air film force. The stiffness is improved slightly by further increasing its thickness. Therefore, the range of 15–30 mm is recommended for the thrust plate thickness for higher stiffness.

**Figure 7.** The performance of aerostatic spindle with varying thrust plate thickness *L* and air film thickness *h*. (**a**) The axial stiffness of the spindle with different thrust plate thickness *L* and air film thickness *h* and (**b**) the first-order natural frequency of the aerostatic spindle with varying thrust plate thickness *L* and air film thickness *h*.

However, the mass of the thrust plate is also directly related to its thickness. To increase the natural frequency, the designer generally hopes to reduce the mass of the rotating parts in the design stage. It means that the thickness of the thrust plate cannot be increased blindly. As shown in Figure 7b, with the increase of the thrust plate thickness, the first-order natural frequency first increases sharply and then declines quickly. A maximum region for the first-order natural frequency can be observed around the thrust plate thickness of 10–15 mm. This phenomenon can be explained as follows. According to the calculation result shown in Figure 7a, the stiffness first increases sharply and then trends to steady, while the mass increases linearly with the increase of the thrust plate thickness. Generally, the natural frequency is proportional to the stiffness and is inversely proportional to the mass. In the region where the stiffness increases sharply, the natural frequency increases with the increase of the thrust plate thickness. In the region where the stiffness increases slowly, the natural frequency decreases with the increase of the thrust plate thickness.

Considering the stiffness and natural frequency comprehensively, a thrust plate thickness of 20 mm is selected in this paper. Comparing with the initial design, the first-order natural frequency of the aerostatic spindle increases by 45.3%, from 167.4 Hz to 243.3 Hz.

#### **6. Experimental Validation of the Calculation Result**

To validate the above calculation result, the performance of the aerostatic spindle with optimized design parameters is measured in this section. Since the machining error and assembly error is unavoidable, the actual geometry size of the aerostatic bearings may deviate from the design value. The performance of aerostatic bearing is susceptible to orifice diameter and air film thickness. To minimize the errors introduced by the manufacturing and assembly process on the measurement results, the orifice diameter and air film thickness are first measured in this section.

#### *6.1. Measurement of Orifice Diameter*

An optical image measuring instrument (Hexagon Optiv Classic 321GL tp) is employed in this research for orifice diameter measuring, as shown in Figure 8a. It has a CCD (Charge Coupled Device) image system to capture the image of the tiny object to be measured, which is capable to measure the contour and dimensions of various complex parts accurately and efficiently. Figure 8b presents a tested image of the orifice. As can be seen from the figure, the round outline of the orifice is clear, and no obvious burr can be observed. These orifices were acquired by micro-drilling machining. The image of the orifice indicates that excellent shape accuracy can be acquired by micro-drilling machining.

**Figure 8.** Experimental setup for orifice diameter measuring. (**a**) The Experimental apparatus for orifice diameter test and (**b**) the tested image of an orifice.

Table 2 lists the tested data of eight orifices. It shows that the actual diameter of the orifice diameter has a good consistency. The nominal size of the orifice diameter is 0.15 mm, and all the measured orifice diameter is around 0.15 mm with a dimension deviation less than ±5 μm, which indicates that the orifices have high dimensional accuracy.

**Table 2.** Experimental result of the orifice diameter.


#### *6.2. Measurement of Air Film Thickness*

Generally, to ensure the static performance of the aerostatic bearing and the motion accuracy of a spindle, a serial of stringent machining accuracy requirements (including flatness, parallelism, and perpendicularity, etc.) are proposed for each part of the spindle system (the shaft, sleeve, and thrust plates). Nevertheless, the real air film thickness may deviate from the design value due to unavoidable machining errors and assembly errors.

Figure 9a presents the experimental setup for air film thickness measuring. By measuring the displacement of the thrust plate under the condition of the air supply valve closing and opening, the actual air film thickness can be acquired. A capacitive displacement sensor (Micro-Epsilon capaNCDT6500) with a resolution of 1 nm and a measurement range of 50 μm is employed in this research to record the displacement of the thrust plate. The experimental setup is mounted on a vibration isolation worktable to avoid the impact of external vibration.

Figure 9b shows the displacement curve of the thrust plate before and after the opening of the air supply valve. It can be seen that the thrust plate moves upward due to the support of pressurized air film. The above measuring process is implemented at eight points. These points are equally spaced along the circumferential direction of the

thrust plate, and the tested data is detailed in Table 3. It can be inferred that the measured average air film thickness is about 12.84 μm, which is slightly less than the design value of 13 μm. It may result from the machining error of the shaft, sleeve, and thrust plate.

**Figure 9.** Experimental setup for air film thickness measuring and test result. (**a**) The Experimental setup for air film thickness measuring and (**b**) a test result of the air film thickness test.

**Table 3.** Measuring result of the air film thickness.


#### *6.3. The Axial Stiffness Test of the Aerostatic Spindle*

To verify the calculation result, the axial stiffness of the spindle with different thrust plate thicknesses is tested in this section. The stiffness test is conducted by applying different loads along the axial direction of the spindle and measuring the axial displacement of the spindle by a displacement sensor. The experimental apparatus mainly consists of a computer, a displacement sensor, and an aerostatic spindle, which is shown in Figure 10. The displacement of the thrust plate can be acquired by the displacement sensor with varying axial loading.

**Figure 10.** Experimental apparatus for the stiffness test.

The experimental result is compared with the simulation result, as shown in Figure 11. It can be seen that the experimental curve and the simulation curve show the same trend. According to the test date, the axial stiffness of the spindle with the thrust plate thickness

of 10 mm is 176.2 N/μm, which is much lower than the ideal stiffness of 209.6 N/μm. It proves that the stiffness will be declined due to structural deformation. By increasing the thrust plate thickness to 20 mm, the stiffness increases to 201.8 N/μm. It indicates that the structural rigidity of the thrust plate can be improved by increasing its thickness. With further increase of the thrust plate thickness to 30 mm, the axial stiffness of the spindle increases slightly to 204.2 N/μm. It means that the stiffness is close to the ideal stiffness with the further increase of the thrust plate thickness. The experimental result validates the accuracy of the FSI model.

**Figure 11.** Comparison between the simulation result and experimental result.

#### **7. Conclusions**

This paper proposes a two-round optimization design method for aerostatic spindles considering the FSI effect. The optimal design of an aerostatic spindle is implemented as a case study to detail the design process of the proposed approach. A set of optimal design parameters are acquired, and experiments are implemented to verify the simulation results. The main conclusions are drawn as follows:

1. A two-round optimization design method is proposed to facilitate the design of the aerostatic spindle with consideration of the FSI effect, and an aerostatic spindle is optimized as a case study that validates the effectiveness of the proposed design method;

2. In the first-round design stage of the two-round optimization design method, a set of optimized design parameters of the aerostatic thrust bearing is acquired by analyzing the performance of the aerostatic thrust bearing with varying geometrical parameters of the restrictor. Compared with the initial design, the LCC is improved by 28.3%, the stiffness is improved by 44.2%, and the VFR is reduced by 27.2%;

3. In the second-round design stage of the two-round optimization design method, a recommended range for the design of thrust plate thickness is acquired by calculating the stiffness and natural frequency of the spindle with different thrust plate thicknesses. Compared with the initial design, the first-order natural frequency of the aerostatic spindle increases by 45.3%;

4. Experiments are conducted, and the experimental results agree well with the simulation results, which verifies the accuracy of simulation results.

**Author Contributions:** Conceptualization, Q.G.; methodology, L.L.; software, Q.G. and S.G.; validation, M.Z. and F.Z.; formal analysis, Q.G. and S.G.; investigation, Q.G. and S.G.; data curation, M.Z.; writing—original draft preparation, Q.G.; writing—review and editing, L.L.; visualization, Q.G.; funding acquisition, Q.G. and S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Heilongjiang Postdoctoral Foundation, grant number LBH-Z19149, the Leading Wild Goose Program Scientific Research Project of Harbin Institute of Technology, grant number XNAUEA5640202220-07, the National Natural Science Foundation of China, grant number 51705501, and the State Key Laboratory of Mechanical Transmissions in Chongqing University, grant number SKLMT-KFKT-201801.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on reasonable request from the corresponding author.

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

#### **References**


## *Article* **Thermal Characteristics Study of the Bump Foil Thrust Gas Bearing**

**Xiaomin Liu 1,2, Changlin Li 1,2,\*, Jianjun Du 1,2,\* and Guodong Nan <sup>1</sup>**


**Abstract:** In this paper, a thermo-hydrodynamic model of the bump foil thrust gas bearing is developed, which solves the coupled gas film three-dimensional energy equation, non-isothermal Reynolds equation, and the foil deformation equation. The effects of bearing speed, thrust load, and external cooling gas on the bearing temperature field are calculated and analyzed. The test rig of foil thrust gas bearing was built to measure the bearing temperature under different working conditions. Both simulation and experiment results show that there exist temperature gradients on the top foil both in the circumferential and radial directions. The simulation results also shows that the top foil side of the gas film has the highest temperature value in the entire lubrication field, and the position of highest temperature moves radially inward on the thrust plate side as the rotor speed increases. The gas film temperature increases with the increasing rotor speed and bearing static load, and rotor speed has greater effects on the temperature variation. Cooling air flow passing through the bump foil is also considered in the simulations, and the cooling efficiency decreases as the mass of gas flow increases.

**Keywords:** bump foil gas bearing; thermo-hydrodynamic model; temperature field; thermal characteristics analysis

#### **1. Introduction**

Dynamic bump foil gas bearings have been developed since 1960s, and numerous scholars have conducted the researches of related lubrication theories and experiments [1,2]. Nowadays, requirements of larger power density performance are put forward for the turbomachinery, such as blowers, air compressors, and turbochargers, which indicates that gas dynamic thrust bearings need to work under the conditions of higher speeds, larger axial loads, and harsher environments [3]. Compared with traditional gas bearings, bump foil type gas bearings have many significant advantages, such as long operating life, high reliability, large carrying capacity, high speed and high temperature resistance, and good impact resistance. As the rotational speed gradually increases, the elastic foil structure of the bearing deforms, thereby forming the corresponding air film thickness, which has strong adaptability. Therefore, more and more researchers pay attention to the dynamic bump foil gas bearing, which are widely used in a variety of high-speed rotating machinery [4]. Under operating conditions, the ultra-high speed of the rotating thrust plate drives the gas in the clearance to move fast. Due to the viscous friction and gas compression effect between the gas films, the temperature of the bearing will increase sharply. It is insufficient to maintain an ideal temperature inside the bearing only by relying on the bearing self-cooling effect. With the accumulation of high temperature, the air film gap will decrease obviously due to the effect of thermal expansion of the structure, and the wear-resistant coating of the top foil will be damaged, which will eventually cause disastrous damages to the whole rotor system.

In recent years, many impressed studies on the thermal characteristics of gas bearings have been published. Salehi et al. applied the Couette approximation method to solve the

**Citation:** Liu, X.; Li, C.; Du, J.; Nan, G. Thermal Characteristics Study of the Bump Foil Thrust Gas Bearing. *Appl. Sci.* **2021**, *11*, 4311. https:// doi.org/10.3390/app11094311

Academic Editor: Terenziano Raparelli

Received: 8 April 2021 Accepted: 30 April 2021 Published: 10 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

temperature field distribution of journal foil gas bearing by coupling the simplified energy equation and the Reynolds equation. However, the calculation model did not consider the heat transfer between the foil structure and rotor [5]. Lee et al. coupled solved the nonisothermal Reynolds equation, gas film thickness equation that considers foil deformation and thermal expansion, and the gas film energy transfer equation. The simulation results of temperature field distribution under the symmetric arrangement of two thrust bearing are obtained, and the influences of rotor speed, bearing load, and cooling air pressure on the bearing temperature are explored [6]. Dickman conducted the thermal-influenced loading experiments based on three corrugated foil thrust bearings with the same size. It was pointed out that the bearing load capacity increases almost linearly with the increase in rotating speed without thermal failure and decreases with the increasing rotor speed when the thermal runaway occurred [7]. GAD et al. calculated the temperature field of thrust bearing by using Couette approximation method, and pointed out that the gas leakage of bearing clearance mostly occurs in the inclined area and leads to the lower temperature. In addition, the study pointed out that more than 70% of the heat generated by gas film can be removed by controlling the cooling flow rate. Meanwhile, the simulation results were compared with Dickman's experimental results, the bearing loading curves in comparison are more consistent when the rotate speed is lower, and there are certain errors when rotor speed is higher. The author speculates that bearing thermal runaway occurred at higher speeds, resulting in a lower load capacity with speed increased [8]. Andreas Lehn et al. present the formulas for calculating the thermal resistances between the bump foil and top foil, and between the bump foil and bearing housing. The authors also considered the thermal expansion of the thrust disc, and the temperature field distribution was solved by multi-field coupling methods. It was pointed out that the bearing load capacity decreased with the increase in rotating speed when the speed increased to a certain value. This trend was due to the bending deformation of thrust disc caused by the high gas film temperature, which leads to the decrease in gas film thickness and the final thermal runaway [9]. Li Changlin established conduction and convection heat transfer models for the bearing sleeve and rotor of the journal bump foil gas bearing, which coupled the non-isothermal Reynolds equations, gas film energy transfer equations, and foil deformation equations. The simulation results compared the THD model with the isothermal model, analyzed the static characteristics and dynamic characteristics of the gas foil bearing [10]. Luo Yixin applied the commercial software to simulate the temperature field of thrust bearing, and studied the influences of different boundary conditions, rotor speeds, inlet pressure, and temperature on bearing thermal characteristics [11]. Peng and Khonsari simplified the foil structure as a line spring model, considered the forced convection cooling outside the flat foil, calculated the multi-dimensional distribution of gas film temperature, but ignored the heat transfer from temperature to rotor and sleeve [12]. Aksoy's research results show that the maximum temperature of the film increases with the grow of load and rotational speed, in which the influence of the rotational speed is greater. At the same time, it is found that the temperature gradient along the axial direction of the rotor can easily cause structural and thermal instability [13]. Andres and Kim found that the viscosity of the gas film increased with the temperature rises, which enhanced the effect of dynamic pressure. The influence of the temperature of the gas film on the structure stiffness of the foil is small. The closer to the shaft end, the greater the gas film temperature gradient, and the axial temperature gradient grew with the increase in the rotational speed [14]. Talmage and Carpino found that the axial temperature gradient will cause the flat foil warping in the end region of the bearing, resulting in uneven axial distribution of the film gap, which makes the flat foil easy to grind with the rotor [15]. Kim and Ki studied the complete heat transfer in a 50KW turbine. It was found that the dynamic stiffness and damping coefficient of the bearing would decrease due to the decrease in the elastic modulus of the foil structure at high temperature [16]. Keun and Andres et al. heated the rotor without cooling gas flow, and the bearing failure occurred when the rotor temperature reached 251 ◦C [17]. Therefore, the heat control and temperature control strategies, such as injecting

cooling gas flow in high temperature environment are very necessary. In this study, a detailed 3D thermal model of bump foil thrust gas bearing is considered. The temperature field distribution of the thrust bearing is solved by coupling the three-dimensional energy equation, non-isothermal Reynolds equation and gas film thickness equation. The cooling models of foil structure and thrust plate, as well as the influence of temperature on the thermal expansion of bearing structure are considered. The thermal characteristics of the bearing under different working conditions are analyzed. Finally, the test rig of the foil thrust gas bearing is built to measure the bearing temperature distributions of bearing under different working conditions.

#### **2. Numerical Simulation Research**

The aerodynamic bump foil thrust bearing mainly consists of three parts: the top foil is used to provide lubricating surface; the bump foil plays the role of elastic support; the bearing housing is applied to fix the top foil and the bump foil. The structure of bump foil thrust gas bearing is shown in Figure 1.

As shown in Figure 2, the entire structural system is divided into six areas along the gas film thickness direction to describe the heat flow. From area A to area G, heat is generated inside the gas film (area C), and one part of the heat is transferred to the ambient or cooling air (area A) through the thrust disc (area B). Another part is transferred to the top foil (area D), and then to the bearing housing (area F) by heat conduction or heat convection (area E) and finally transferred to the open air (area G). When there is cooling air passing through the bump foil channel, only a small part of the heat is transferred to the surrounding air through the bump foil and bearing sleeve.

**Figure 2.** Heat flow transfer model.

#### *2.1. Calculation Process and Parameters*

In this paper, the energy equation of the gas film in cylindrical coordinate system is established as follow:

$$\begin{split} \rho(T)c\_p(T) \left( v\_r \frac{\partial T}{\partial r} + \frac{v\_\theta}{r} \frac{\partial T}{\partial \theta} + v\_z \frac{\partial T}{\partial z} \right) &= k\_d(T) \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial T}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 T}{\partial \theta^2} + \frac{\partial^2 T}{\partial z^2} \right) + \\ \left( v\_r \frac{\partial p}{\partial r} + \frac{v\_\theta}{r} \frac{\partial p}{\partial \theta} \right) + \mu(T) \left( \left( \frac{\partial v\_r}{\partial z} \right)^2 + \left( \frac{\partial v\_\theta}{\partial z} \right)^2 \right) \end{split} \tag{1}$$

where *cp* is the specific heat capacity of gas at constant pressure (J/(kg·K)); *ka* is the thermal conductivity of the gas (W/(m·K)); *μ* is the gas viscosity which changes with temperature, according to the reference [5], the relationship between *μ* and temperature is expressed as:

$$
\mu = a \left( T - T\_{ref} \right) \tag{2}
$$

where the slope *<sup>a</sup>* = 4 × <sup>10</sup><sup>−</sup>8, when the unit of temperature *<sup>T</sup>* is degrees Celsius, *Tref* <sup>=</sup> −458.75. The energy equation is nondimensionalized as follow:

$$\begin{split} \frac{1}{\rho}\overline{\rho}(T)\left(\overline{v}\_{r}\frac{\partial T}{\partial \overline{\rho}} + \frac{\overline{v}\_{\theta}}{\overline{\rho}}\frac{\partial T}{\partial \overline{\theta}} + \frac{1}{\overline{\rho}}\overline{v}\_{z}\frac{\partial T}{\partial \overline{z}}\right) &= k\_{1}\left(\frac{1}{\overline{\rho}}\frac{\partial T}{\partial \overline{\theta}} + \frac{\partial^{2}T}{\partial \overline{\theta}^{2}} + \frac{1}{\overline{\rho}^{2}}\frac{\partial^{2}T}{\partial \theta^{2}}\right) \\ + k\_{2}\frac{1}{h^{2}}\frac{\partial^{2}T}{\partial \overline{z}^{2}} + k\_{3}\left(\overline{v}\_{r}\frac{\partial \overline{\rho}}{\partial \overline{\theta}} + \frac{\overline{v}\_{\theta}}{\overline{\rho}}\frac{\partial \overline{\rho}}{\partial \theta}\right) &+ k\_{4}\frac{1}{h^{2}}\left(\left(\frac{\partial \overline{v}\_{r}}{\partial \overline{z}}\right)^{2} + \left(\frac{\partial \overline{v}\_{\theta}}{\partial \overline{z}}\right)^{2}\right) \end{split} \tag{3}$$

where *k*1, *k*2, *k*3, *k*<sup>4</sup> are expressed as follows:

$$\begin{array}{ll}k\_1 = \frac{k\_{air}}{\rho\_0 c\_p \omega r r\_2^2} & k\_2 = \frac{k\_{air}}{\rho\_0 c\_p \omega \delta\_h^2} \\ k\_3 = \frac{p\_a}{\rho\_0 c\_p \left(T\_0 - 273.15 - T\_{ref}\right)} & \\ k\_4 = \frac{\mu \omega r\_2^2}{\rho\_0 c\_p \left(T\_0 - 273.15 - T\_{ref}\right) \delta\_h^2} \end{array} \tag{4}$$

and the dimensionless coefficients in Equation (3) are defined as follows:

$$\begin{array}{l} \overline{\rho} = \frac{\rho}{\rho\_0}, \overline{\tau} = \frac{r}{\overline{\kappa}\_2}, \overline{\tau} = \frac{z}{h} = \frac{z}{h\delta\_h}, T = \frac{T - T\_0}{\overline{\tau\_0} - 2\overline{\tau}\overline{\lambda}.15 - T\_{ref}}\\ \overline{p} = \frac{p}{p\_a}, \overline{\upsilon\_r} = \frac{v\_r}{\omega r\_2}, \overline{\upsilon\_\theta} = \frac{v\_\theta}{\omega r\_2}, \overline{\upsilon\_z} = \frac{v\_z}{\delta\_h \omega} \end{array} \tag{5}$$

The velocity components *vr*, *vθ*, *vz* could be solved by the continuity equation of hydrodynamics and the boundary conditions assumptions in Reynold's equation.

#### *2.2. Non-Isothermal Reynolds Equation of Gas Film*

When the rotor is running at a higher speed, the heat generated by gas viscous dissipation and compression will increase the gas film temperature. Meanwhile, gas properties such as density and viscosity will change, which can have important impacts on the bearing performance. The non-isothermal gas film Reynolds equation applied in this paper is formulated as:

$$\frac{1}{r}\frac{\partial}{\partial r}\left(rh^3p\frac{\partial p}{\partial r}\right) + \frac{1}{r^2}\frac{\partial}{\partial \theta}\left(h^3p\frac{\partial p}{\partial \theta}\right) = 6\mu(T)\omega\frac{\partial(ph)}{\partial \theta} \tag{6}$$

where *μ*(*T*) represents gas viscosity. In the calculation process, the average temperature of the gas film is used to update its viscosity.

#### *2.3. Gas Film Thickness Equation*

According to the structure form of the thrust bearing, considering the influence of temperature on the thermal expansion of the bearing structure, the gas film thickness distribution is shown in Figure 3, and its equation is as follows:

$$h = h\mathfrak{z} + \mathfrak{z}(\theta) + w(r, \theta) - \Delta\_{\text{thermal}}\tag{7}$$

$$\mathcal{S}(\theta) = \begin{cases} & \delta\_h[1 - \theta/(b\beta)] & 0 \le \theta < b\beta \\ & 0 & b\beta \le \theta \le \beta \end{cases} \tag{8}$$

$$w(r,\theta) = \frac{F(r,\theta)}{K\_{\mathsf{b}} \times \Delta l}, \qquad K\_{\mathsf{b}} = \frac{Et\_{\mathsf{b}}^3}{2(1-v^2)l^3} \tag{9}$$

where *h2* is the gas film thickness in the plane area (m); *g*(*θ*) is the initial gas film thickness distribution which minuses the value of *h2* (m); *β* is the angle of top foil (◦); *b* is the pitch ratio; *w*(*r*, *θ*) is the deformation of top foil under aerodynamic force, and *K*<sup>b</sup> is the stiffness of the bump foil per unit of transverse length (N/m2) [18]; Δthermal is the total thermal expansion of all bearing structures (m).

**Figure 3.** Gas film thickness distribution diagram.

#### *2.4. Bearing Load Capacity and Friction Loss Power*

Load capacity (*F*) and friction loss power (*W*) are the main static characteristics to measure bearing performance, their calculation formulas are as follows:

$$\mathcal{F} = \frac{F}{p\_a r\_2^{\prime 2}} = \int\_0^{\S} \int\_{r\_1/r\_2}^1 (\mathfrak{P} - 1) \mathfrak{F} d\mathfrak{F} d\theta \tag{10}$$

$$\overline{T\_r} = \frac{T\_r}{p\_a h\_2 r\_2} = \int\_0^\beta \int\_{r\_1/r\_2}^1 \frac{\overline{r}\overline{h}}{2} \frac{\partial \overline{p}}{\partial \theta} + \frac{\Lambda}{6} \frac{\overline{r}^3}{\overline{h}} d\overline{r} d\theta \tag{11}$$

$$
\overline{W} = \frac{2\pi n}{60} \times \overline{T\_r} \tag{12}
$$

where *n* is the rotational speed of the thrust plate, *Tr* is friction torque.

#### *2.5. Heat Transfer Model of Bearing Components*

2.5.1. Heat Transfer Model of Top Foil Side

In this paper, the equivalent thermal resistance is applied to calculate the overall heat transfer processes when cooling air is introduced or not.

Figure 4 is the heat transfer model on the top foil side. The figure on the left shows that when cold air is not passed into the foil, the heat generated in the gas film passes through the top foil through conduction. Part of the heat is conducted to the base plate through the bump foil, and the other part is conducted to the base plate through the air in the gap between foils. As the base plate is in close contact with the bearing housing, the heat is transferred to the bearing housing through heat conduction, and finally spread out through natural convection. The figure on the right shows that when the cooling air flow is introduced into the foil, the heat generated in the gas film transferred through a different way. The heat passes through the top foil through conduction and part of it is transferred to the bump foil through heat conduction, and another part is taken away by the cold air directly. The heat transferred to the bump foil is conducted to the base plate partly, and the rest of it is taken away by the cold air which passed through the upper and lower surfaces of the bump foil. Then, part of the heat transferred to the base plate is taken away by the cooling air flow on the upper surface of the base plate through convection, the other part is transferred to the bearing housing through heat conduction, and finally spread out through natural convection. Correspondingly, the equivalent thermal resistance models of these two situations are established, respectively.

**Figure 4.** Heat transfer paths near the top foil.

Figure 5 is the schematic diagram of thermal resistance model established based on the ways of heat transfer in the foil structures, where *T*<sup>1</sup> is the air film temperature on the top foil side; *T*<sup>0</sup> is the ambient temperature. The remaining parameters are the corresponding equivalent thermal resistance.

When there is no cooling air flow injected, the equivalent total thermal resistance is calculated as follows:

$$R\_a = R\_-t + \frac{1}{\frac{1}{R\_-fp} + \frac{1}{R\_-fb}} + R\_-b + R\_-w + R\_-e \tag{13}$$

when the cooling air flow is introduced into the bump foil, the equivalent total thermal resistance is calculated as follows:

$$R\_a = R\_-t + \frac{1}{\frac{1}{R\_-bp + \frac{1}{\frac{1}{R\_-bhc} + \frac{1}{R\_-t}}} + \frac{1}{R\_-c}}\tag{14}$$

where *R\_i* is calculated as follows:

$$R\_-i = R\_-b + \frac{1}{\frac{1}{R\_-bb} + \frac{1}{R\_-w + R\_-\epsilon}}\tag{15}$$

the calculation formula for each thermal resistance is as follows:

$$\begin{aligned} R\_{\underline{\boldsymbol{t}},\boldsymbol{t}}R\_{\underline{\boldsymbol{h}}\boldsymbol{p}}\boldsymbol{R}\_{\underline{\boldsymbol{t}}}\boldsymbol{t}b, R\_{\underline{\boldsymbol{h}}\boldsymbol{h}}\boldsymbol{R}\_{\underline{\boldsymbol{h}}}\boldsymbol{w} &= \boldsymbol{t}\_{\boldsymbol{i}}/k\_{\boldsymbol{i}}\boldsymbol{A}\_{\boldsymbol{i}}\\ R\_{\underline{\boldsymbol{c}}}\boldsymbol{c},\boldsymbol{R}\_{\underline{\boldsymbol{h}}}\boldsymbol{h}\boldsymbol{c},\boldsymbol{R}\_{\underline{\boldsymbol{h}}}\boldsymbol{h}b, R\_{\underline{\boldsymbol{c}}}\boldsymbol{c} &= 1/A\_{\boldsymbol{j}}h\_{\underline{\boldsymbol{j}}} \end{aligned} \tag{16}$$

where *i,j* are the serial number of the thermal resistance; *t* is the thickness on heat transmission direction (m); *k* is the thermal conductivity of the material (W/K·m); *h* is the coefficient of heat convection(W/m2·K); *<sup>A</sup>* is the heat transmission area (m2). Under normal circumstances, the parameters *t* and *A* are determined by the geometric features of the components, exceptionally, the *Ai* of *R\_bp* is contact area between the bump foil and base plate, and the *Aj* of the *R\_hbc* is approximated obtained by straight inclined segments of the bump foil.

**Figure 5.** Equivalent thermal resistance model diagram.

When the gas film temperature reaches equilibrium, the heat transferred into the top foil from the gas film will be in dynamic balance with the heat transferred from the top foil side to the outside. Therefore, the heat balance equation between *T*<sup>1</sup> and *T*<sup>0</sup> considering the total thermal resistance of the foil structures *Ra* will form the:

$$-k\_d A \frac{\partial T\_1}{\partial\_Z} = \frac{T\_1 - T\_0}{R\_d} \tag{17}$$

#### 2.5.2. Heat Transfer Model of the Thrust Plate

Figure 6 is the assembly diagram of a single side thrust bearing. Considering that the thrust plate has high rotating speed, it is assumed that the surface temperature of the thrust plate is consistent along the circumferential direction in the calculation process.

The governing equation of the one-dimensional thrust plate heat conduction model in the radial direction is formulated as:

$$\frac{1}{r}\frac{\partial}{\partial r}\left(k\_{\rm disc}r\frac{\partial T\_{\rm disc}}{\partial r}\right) + \frac{h\_{\rm conv\\_disc}}{t\_{\rm disc}}(T\_{\rm disc} - T\_0) + \overset{\bullet}{q}\_{film\\_disc} + \overset{\bullet}{q}\_{disc\\_sur} = 0\tag{18}$$

where *kdisc* is thermal conductivity of thrust plate (W/K·m); *Tdisc* is the thrust plate surface temperature (K); *hconv\_disc* is the convection coefficient of air and thrust plate (W/K·m); *tdisc* is the thrust plate thickness (m).

**Figure 6.** Single side gas thrust bearing configuration.

The last two items represent the heat flux density from the gas film to the thrust plate and from the thrust plate to the ambient environment through the inner and outer radial, respectively. The heat transferred from the air film to the thrust plate is expressed as:

$$
\dot{q}\_{film\\_disc} = \frac{k\_{air}}{t\_{disc}} \frac{dT}{dz} \tag{19}
$$

In this paper, when considering the heat transfer condition on the other side of the thrust plate, the heat dissipation model of the backside of the rotor thrust plate is regarded as the convection heat dissipation of the rotating plate in static air. When calculating the heat dissipation coefficient of convection between the back of the thrust plate and the environment, the calculation formula in reference [19] is used in this article:

$$Nu = 0.015 \text{Re}\_{\omega}{}^{0.8} \tag{20}$$

$$h\_{conv\\_disc} = N\_u \times k / L \tag{21}$$

where *Nu* represents the Nusselt number of the fluid, which is used to describe the convective heat transfer intensity, and Reω represents the Reynolds number, and *k* is the thermal conductivity of thrust plate (W/K·m), and *L* is the characteristic length (m).

In order the calculate the convective heat dissipation from the thrust plate outer surface to the surrounding air, the convection coefficient on the outer surface of a rotating cylinder is used in this paper [20]:

$$h\_{conv\\_out} = 0.133 \text{Re}\_2^{2/3} \text{Pr}^{1/3}(k\_{air}/D\_R) \tag{22}$$

where Re2 is the Reynolds number with *DR* being the characteristic length, and *DR* is equal to 2*r*2; *r2* is the outer radius of the thrust bearing (m); Pr is the Prandtl number of air; *kair* is the thermal conductivity of gas (W/K·m).

. *q*disc\_sur is the heat flux density transferred from the thrust plate to the environment in the radial direction, which is calculated as follows:

$$
\dot{q}\_{\text{disc\\_surr}} = \frac{h\_{\text{conv\\_out}}}{\Delta r} (T\_{\text{disc}} - T\_0) \tag{2.3}
$$

2.5.3. Temperature Boundary Condition in the Gas Film Inlet Area

For bump foil thrust gas bearings, the pressure in the gas film area is higher than the ambient atmospheric pressure, thus leading to the leakage of air from the radial boundaries. Meanwhile, the same amount of air Q*suc* is brought into the gas film at the top foil leading edge, which will mix with recirculation air Q*rec* to achieve the conservation of the total amount of gas Q*in*.

The temperature of the inlet mixed air is calculated as follows based on the energy balance condition of the controlled gas volume:

$$T\_{in} = \frac{T\_{\rm rec}Q\_{\rm rec} + T\_{\rm succ}Q\_{\rm succ}}{Q\_{\rm rec} + Q\_{\rm succ}} \tag{24}$$

where *Tin* is the mixed gas temperature at the leading edge of the top foil (K); *Trec* is the circulating air temperature at the rear edge of the top foil (K); *Tsuc* is the temperature of the intake cooling air (K).

#### **3. Thermal Analysis**

*3.1. Calculation Process and Parameters*

Figure 7 shows the calculation algorithm of the air-thermoelastic coupling simulation, and Table 1 lists the bearing parameters used in the numerical and experimental studies.

**Figure 7.** Algorithm of the air-thermoelastic coupling simulation.


**Table 1.** This is a table. Tables should be placed in the main text near to the first time they are cited.

#### *3.2. Analysis of Influence Factors of Bearing Temperature*

Figures 8 and 9 present the top foil deformations and the dimensionless gas film pressure distributions when the rotor speed is 50 krpm and the initial minimum gas film thickness is 7 μm. Figure 8a,b show the calculation results of the isothermal model and THD model, respectively. It can be seen that the basic variation trends of the top foil deformations calculated by different models are the same, but the top foil deflection calculated by the THD model is larger. This is because the gas film pressure is larger when the temperature effect is considered. The top foil deformation in the inclined section is larger than that of the plane section due to the absence of the corrugated bump foil underneath. The larger foil stiffness in the plane section helps increase the bearing load capacity.

**Figure 8.** Deformation distribution of top foil under different thermal model. (**a**) Isothermal model; (**b**) THD model.

It can be seen from Figure 9 that the pressure distribution under the two calculation models is basically the same. Along the circumferential direction, the gas film pressure increases sharply from the inlet area, and reaches the maximum value when the gas film transits from the inclined area to the horizontal area. Then it slowly decreases, and finally the air pressure drops to atmospheric pressure in the gas outlet area. In the radial direction, both the inner and outer outlet area are atmospheric pressure, so the pressure presents a trend of increasing first and then decreasing from the inner radius to the outer radius boundaries. There is continuous leakage of gas from the radial boundaries when the bearing is under working condition because of the pressure distribution characteristics. Compared with the isothermal model, the gas film pressure calculated by the THD model is larger than that by the isothermal model. This is due to the fact that the gas viscosity

continues to increase as the temperature rises, and the gas film thickness is reduced to a certain extent due to the thermal expansions of the thrust plate and the foil structures.

**Figure 9.** Dimensionless gas film pressure distribution under different thermal model. (**a**) Isothermal model; (**b**) THD model.

Figure 10 shows the change of bearing load capacity and friction loss power with the initial minimum gas film thickness calculated by the isothermal and THD models. It can be seen from the Figure 10a that bearing load capacity decreases with the increase in the initial minimum gas film thickness. When the initial minimum gas film thickness is increased to 30 μm, the gap between the thrust plate and the top foil is too large. Thus, the lubricating gas film is difficult to form, and the bearing almost has no load capacity. It can be seen from Figure 10b that the frictional loss power decreases as the initial minimum gas film thickness increases, which is because the gas film pressure is very small when the initial minimum gas film thickness is large. However, there is no obvious difference between the values of friction loss power calculated by the isothermal model and the THD model. When the initial minimum gas film thickness is the same, the bearing load capacity and friction loss power both increased when the rotor speed goes up, and the difference between the values of the bearing load capacity calculated by the THD model and the isothermal model increased significantly when the speed reach to higher value. This is due to the rapid rise of the gas film temperature caused by the increased rotational speed, and the gas viscosity increases with the temperature rising up.

**Figure 10.** Static characteristics under different thermal model. (**a**) Load capacity; (**b**) Friction loss power.

Figure 11 is the temperature distribution diagram of top foil under the rotor speed of 50 krpm and the static load of 100 N. Similar with the distribution of gas film pressure, the surface temperature of the top foil experiences, firstly, increasing along the rotor rotating direction and then maintaining at a certain value. In the entrance area, the temperature of the ambient gas that is brought into the gas film is relatively low, and it will rise sharply after mixing with the recirculation hot air. The top foil temperature in the plane section is almost the same, which is obviously higher than that in the inclined region. In the radial direction, the temperature increases first and then decreases slightly as the radius increases. This is because the larger radius indicates the larger linear velocity of the gas film and the larger amount of heat generated due to viscous friction. The decreasing trend near the edge can be explained as the heat is dissipated by the forced convection on the thrust plate outer radius surface.

**Figure 11.** Top foil temperature distribution.

Figure 12 shows the relationship between bearing temperature and the static load at the rotor speed of 50 krpm. The temperatures of top foil and thrust plate are the average values in the simulations. It can be seen from the figure that as the bearing load increases, the gas film temperature increases almost linearly. Additionally, when the bearing load increases, the differences between the maximum air film temperature, the top foil temperature and the thrust plate temperature are almost kept constant. The increase in bearing static load means the gas film is compressed, and the heat generated by the viscous dissipation and the compression in the gas film is larger, resulting in the rise of overall temperature.

**Figure 12.** Relationship between bearing temperature and load.

Figure 13 shows the variation of bearing temperature influenced by rotor speed under the bearing static load of 100 N. As the bearing speed increases, the temperatures of gas film, top foil, and thrust plate all increase exponentially. Therefore, rotor speed has a greater influence on the bearing temperature compared with the influence of static load. At low rotor speeds, the difference between the top foil and the thrust plate average temperatures is relatively small, and the maximum temperature of gas film is also relatively small. As the rotational speed increases, the differences between the temperatures of the top foil, thrust plate, and the maximum gas film temperature also increase.

Figure 14 shows the temperature distributions at various radial positions of the thrust plate at different rotor speeds under the same static load of 100 N. It shows that the temperature of the thrust plate first increases and then decreases in the radial direction, and the highest temperature position is closer to the outer radius boundary. With the increase in rotor speed, the distance from this position to the outer radius boundary tends to become larger.

**Figure 14.** Relationship between thrust plate temperature and speed.

This is because part of the heat is dissipated through forced convection from the cylindrical surface at the thrust plate outer radius, so the gas film temperature will decrease in this area to some extent. As the rotor speed increases, the convective heat exchange effect between the outer edge of the thrust plate and the environment becomes much stronger.

Figure 15 shows the relationship between bearing temperature and the flow rate of cooling gas under the rotor speed of 100 krpm and static load of 100 N. The results show that when the flow rate of cooling gas gradually changes from 0 to 0.2 m3/min, the bearing temperature drops rapidly. If the gas flow rate is further increased, the cooling effect will be weakened gradually. Therefore, in engineering applications, the appropriate flow of cooling gas passed into the foil structure should be designed to decrease the bearing temperature with the maximum efficiency.

**Figure 15.** Relationship between bearing temperature and cooling gas flow.

#### **4. Experimental Study**

In this paper, the test specimen of the bump foil thrust gas bearing is manufactured and assembled, as shown in Figure 16. The K-type thermocouple is applied to measuring the bearing temperature. A number of thermocouples are stick to the top foil on the back side with copper foil tape. Figure 17 shows the positions of these thermocouples for measuring the local temperature. The thermocouple wire is buried in the bump channel.

The same number and position of thermocouples are pasted under the two centrosymmetric top foils. The temperature measuring points 1–4 are arranged at the radius of 30 mm, and points 5–6 are arranged at the radius of 25 mm. The temperature measuring point 1 is arranged under the slope section, and the remaining temperature measurement points are in the loading region. The assembly of the bearing is completed by a spot-welding machine.

The test bench in this paper is shown in Figure 18. The actual working process of the bearing can be realized by actively controlling the speed of the motorized spindle and by loading the thrustmeter on the other side of the thrust bearing. The displacement between the rotating plate and the bearing is measured by the eddy current displacement sensor.

**Figure 16.** Physical image of experimental bearing.

**Figure 17.** Temperature measuring point bonding.

**Figure 18.** Thrust bearing temperature test bench.

In this experiment, the axial load, which is provided by the SN-500 pointer thruster produced by SUDOO company, is applied on the side of the thrust bearing to simulate the actual axial force of the system. The structure of pointer thruster is shown in Figure 19, where the limit load is 500 N and the resolution is 2 N, the pointer thruster is fixed to the antivibration platform through the supports, and the axial force to the bearing is applied by rotating the button.

**Figure 19.** Picture of pointer thruster.

Figure 20 shows the relationship between the tested top foil temperature and the bearing static load at the rotor speed of 15 krpm. During the loading process, the bearing temperature increases continuously. It can be seen from the figure that the measured temperatures on points 2, 3, and 4 are significantly higher than the temperatures on points 5 and 6, and the temperature on point 1 in the slope section is lower compared with the temperatures measured on the points of the same radius. These distribution characteristics of the foil thrust bearing temperature are consistent with the simulation results. Compared with the simulation results, the measured top foil temperatures are relatively lower, which is mainly due to the unclear boundary conditions in the experiments which should be clarified in the simulations.

**Figure 20.** Relationship between bearing top foil temperature and load.

Figure 21 shows the relationship between the measured top foil temperature and rotor speed under the static load of 80 N. It shows that the temperature distribution regularities of the measuring points are similar with the regularities demonstrated in Figure 20, and the overall trend of the measured temperature is the same with the simulation results, which is the temperature will rise up almost linearly with the rotational speed going up. However, the experimental data are relatively higher than the simulation data, which is probably caused by the unclear boundary conditions in the experiments which should be clarified in the simulations or due to the measurement error during the experiment.

**Figure 21.** Relationship between bearing top foil temperature and rotational speed.

#### **5. Conclusions**

In order to study the thermal characteristics of the aerodynamic bump foil thrust bearing, the air-thermoelastic coupling model is developed in this paper with the heat transfer models of the bearing and thrust plate, and the thermal expansion effect of the bearing structure being considered. The temperature field of the gas thrust bearing is solved through the multi-field coupling algorithm. Finally, the experiments of bearing temperature measurement are carried out under different working conditions. Based on the simulation results, the following conclusions are obtained:


The current simulation ignored some factors, which makes the simulation results are not exactly the same with the experiment results. The further research will consider more accurate structure model and thermal model, including the surface roughness of the foils and the friction between them, and more complex heat transfer models, etc.

**Author Contributions:** Conceptualization, J.D. and X.L.; methodology, X.L.; software, X.L.; validation, X.L., C.L., J.D. and G.N.; formal analysis, J.D.; investigation, G.N.; resources, J.D.; data curation, J.D.; writ-ing—original draft preparation, X.L.; writing—review and editing, X.L.; visualization, C.L.; supervision, C.L.; project administration, J.D.; funding acquisition, C.L. and J.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work described in this paper was supported by National Key R&D Program of China (Grant No. 2018YFB2000100) and the National Natural Science Foundation of China (Grant No. 52005126).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on reasonable request from the corresponding author.

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

#### **Nomenclature**


#### **References**


## *Review* **Application of Gas Foil Bearings in China**

**Yu Hou \*, Qi Zhao, Yu Guo, Xionghao Ren, Tianwei Lai and Shuangtao Chen**

School of Energy and Power Engineering, Xi'an Jiaotong University, Xi'an 710049, China; wheredream2012@stu.xjtu.edu.cn (Q.Z.); yu\_guo@stu.xjtu.edu.cn (Y.G.); rxh0206@stu.xjtu.edu.cn (X.R.); laitianwei@mail.xjtu.edu.cn (T.L.); stchen.xjtu@mail.xjtu.edu.cn (S.C.) **\*** Correspondence: yuhou@mail.xjtu.edu.cn

**Abstract:** Gas foil bearing has been widely used in high-speed turbo machinery due to its oil-free, wide temperature range, low cost, high adaptability, high stability and environmental friendliness. In this paper, state-of-the-art investigations of gas foil bearings are reviewed, mainly on the development of the high-speed turbo machinery in China. After decades of development, progress has been achieved in the field of gas foil bearing in China. Small-scale applications of gas foil bearing have been realized in a variety of high-speed turbo machinery. The prospects and markets of high-speed turbo machinery are very broad. Various high-speed turbomachines with gas foil bearings have been developed. Due to the different application occasions, higher reliability requirements are imposed on the foil bearing technology. Therefore, its design principle, theory, and manufacturing technology should be adaptive to new application occasions before mass production. Thus, there are still a number of inherent challenges that must be addressed, for example, thermal management, rotor-dynamic stability and wear-resistant coatings.

**Keywords:** foil bearing; turbo machinery; gas bearing; industrial application

#### **1. Introduction**

Gas foil bearing has the distinct advantages of cleanliness, low friction, low temperature rise, and cost effectiveness [1]. These advantages make it applicable in most advanced turbo machinery systems, for instance, aerospace, energy and power engineering, etc. [2–4].

Compared with sliding bearing and ball bearing, gas foil bearings have a high running speed and a wider operating range, from tens of thousands rpm (round per minute) to hundreds of thousands rpm. Its power ranges from several kilowatts to more than 100 kilowatts, as shown in Figure 1. It is especially suitable for micro–small high-power density turbo machinery, such as fuel cell air compressors, turbo-expanders, oil-free turbochargers, etc. In comparison, sliding bearings are favored to be used in low-speed applications due to their oil lubricating properties, which are an important support component in large capacity power machinery [5]. As for magnetic bearings, their speed is between that of foil bearings and sliding bearings. Most are used in turbine generator, aeroengine, high-speed CNC (computer numerical control) machine tools, wind power equipment and other cases.

For gas foil bearings, their stability and reliability depend largely on elastic support. Under the elastic top foil, an elastic supporting structure provides additional damping, which conforms to the variable load and speed [6]. Through the damping effect, vibration energy can be absorbed due to Coulomb friction dissipation between the foils and bearing housing. Even if there is impact force or unstable whirl motion, the rotor maintains its stability [7].

In order to improve the bearing performance, different types of gas foil bearings have been proposed. According to the underlying elastic structure, the most-used gas foil bearings are bump foil type, multi-leaf type, protuberant type, and so on [8]. The lubrication mechanisms of a gas foil journal and thrust bearings are basically the same. Herein, a gas foil journal bearing is chosen as an illustration, shown in Figure 2. The basic

**Citation:** Hou, Y.; Zhao, Q.; Guo, Y.; Ren, X.; Lai, T.; Chen, S. Application of Gas Foil Bearings in China. *Appl. Sci.* **2021**, *11*, 6210. https://doi.org/ 10.3390/app11136210

Academic Editors:

Terenziano Raparelli, Andrea Trivella, Luigi Lentini and Federico Colombo

Received: 30 May 2021 Accepted: 30 June 2021 Published: 5 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

structure of a bump foil bearing is shown in Figure 2a [9], which consists of a flexible flat foil and an elastic bump foil. The bumps of the supporting foil are evenly distributed along the circumferential direction. The top foil is stacked over corrugated foil, and both foils are fixed at one end and are free at the other end. As for multi-leaf foil bearing, it is composed of a bearing housing and several cantilever foils [10], as shown in Figure 2b. The leading edge of the cantilever foil is fixed in the bearing, while the trailing edge is freely overlapped on the next cantilever foil. The extension direction of the cantilever foil is consistent with the rotation direction of the rotor. A wire can also be used for elastic support. In the wire supporting foil bearing, a certain number of copper wires are fixed on the back of the foil, as shown in Figure 2c [11]. The rigidity of the foil element can be adjusted by changing the number of copper wires and thicknesses of the copper wires. A viscoelastic supporting gas foil bearing is another kind of foil bearing, as shown in Figure 2d. In 1999, the Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University proposed a viscoelastic supporting gas foil bearing [12]. The foil element is a combination of a metal foil and a piece of viscoelastic supporting sheet. The spring supported foil bearing is a new type of foil bearing, as shown in Figure 2e. Hunan University presented a theoretical and experimental study on a novel nested compression spring [13,14]. The stiffness and damping of the support can be adjusted by varying the number and size of springs. In the protuberant foil structure, as shown in Figure 2f, there are multiple layers of protuberant foil between the top foil and the bearing housing [15]. The supporting foil element is composed of a top foil and two protuberant foils. In the assembly process, the bearing performance and stability of the bearing can be improved and the wear on the bearing surface is reduced due to the effective elastic pre-deformation of the protuberant foils [16,17].

**Figure 1.** Application range of different kinds bearings.

**Figure 2.** Typical gas foil journal bearings. (**a**) bump foil bearing; (**b**) multi-leaf foil bearing; (**c**) wire supporting foil bearing; (**d**) viscoelastic supporting gas foil bearing; (**e**) spring supported foil bearing; (**f**) protuberant foil bearing.

Before large scale industrial application of the gas foil bearing can be achieved, there are still some practical problems that should be considered. For example, gas foil bearings need to operate stably between takeoff speeds and landing speeds. During starting and stopping processes, the foil will wear poorly in the long run due to the contacting dry friction and the large friction torque. On the other hand, larger power density turbomachinery might lead to a larger axial force, which imposes harsh working requirements on the thrust bearing [18]. In addition, due to the interaction of the foil support structure and the nonlinearity, challenges still exist in accurate theoretical analyses and numerical predictions on foil bearing performance.

The research on gas foil bearings has been carried out relatively early abroad, and mainly used in micro gas turbines [19,20], turbo-expanders [21,22], turbochargers [23–25], fuel cell air compressors machines [26,27], air cycle machines [28], and other options [29]. In recent years, in order to catch up with advanced technology, gas foil bearing technology has been developing rapidly in China. In the 1970s, domestic scholars conducted research on double-row orifice hydrostatic gas journal bearings. Subsequently, a large number of theoretical and experimental studies were carried out in other forms of radial and thrust gas bearings. A series of in-depth studies were conducted in Xi'an Jiaotong University and other institutions, such as, Harbin Institute of Technology, and Beijing Precision Engineering Institute Aircraft industry, etc. [30] For example, a low-temperature turboexpander using rubber-stabilized hydrostatic gas bearings was successfully developed by Xi'an Jiaotong University in 1981. Currently, a series of theoretical and experimental research experiments on gas foil bearings has been conducted at Xi'an Jiaotong University [31], Hunan University [32,33], the Harbin Institute of Technology [34], Nanjing University of Aeronautics and Astronautics [35–37], the GREE company, and other institutions [38]. Although some progress has been achieved in foil bearings and their applications, a certain number of breakthroughs are still anticipated before large-scale application. In this paper, the progress of domestic gas foil bearings and their applications in high-speed turbo machinery are reviewed, which could provide a reference and guidance for further developments and applications of gas foil bearings.

#### **2. Applications of Gas Foil Bearings**

#### *2.1. Application in Low Temperature High-Speed Turbo-Expander*

The working temperature of a cryogenic expander is always lower than 120 Kelvin. In order to enlarge the mass flow rate and maintain a relatively high adiabatic efficiency, its rotor needs to run at a high speed under low temperatures. By using the gas foil bearing, the unstable whirl motion of the rotor bearing system can be suppressed and effectively attenuated. Its stability and reliability can be improved notably. In the 1990s, theoretical and experimental research on air foil bearings for turbo-expanders was initiated by the Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University. The rated air volume of the turbo-expander is 600 Nm3/h@1.15 MPa, and it can work between 0–110 krpm (overspeed to 150 krpm) according to the air supply pressure [39]. The research work on the wire support, viscoelastic support, bump foil support, protuberant support, and multi-leaf foil bearings has been carried out successively, and some progress has been made.

In China, foil bearings were first used in cryogenic turbo-expanders, and then gradually extended to other fields. In 1998, a wire support elastic structure foil bearing was first proposed [40]. In the experiment, a stable speed of 120 krpm or more was reached, and the maximum amplitude was less than 12 μm. The experimental results showed that the turbo-expander ran smoothly at an overspeed of 120% at 118 krpm. In 1999, a turbo-expander supported by viscoelastic foil bearings achieved an excellent performance of 40% over-speeding at 148 krpm [21]. The test used MoS2 powder and two methods of plating Cr and Tic to enhance the hardness of the foil, and over 50 start–stop performance tests were carried out. The test results showed that the surface-treated elastic foil thrust bearing had better start–stop performance and stability. The Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University also conducted a theoretical study on viscoelastic foil bearings [41] and compared the experimental results with wire supported foil bearings. It was found that the stiffness distribution of the viscoelastic bearings was more uniform and had a better stability during ultra-high-speed operations [42]. In 2006, considering the influence of friction, a simplified formula of structural stiffness was obtained, and the structural stiffness of bump foil bearing was analyzed [43]. By comparing two types of foil journal bearings with foil thicknesses of 0.05 mm and 0.07 mm, the vibration characteristics and stability of a turbo-expander were studied. The maximum speed of the turbo-expander using 0.05-mm bump foil bearings was 93,336 rpm. For the turbo-expander using a 0.07-mm bump foil bearing, rotor whirl motion appeared when its maximum speed reached 93,161 rpm. The results showed that when the turbo-expander reached 93 krpm, the maximum amplitude of the rotor was less than 20 μm [44]. In 2012, in order to analyze the feasibility of the application of protuberant foil journal bearings in a turbo-expander, another experimental study was carried out. The results showed that the maximum speed of the turbo-expander reached 99,044 rpm, and the low-frequency whirl motion was small during the entire speed-up and speed-down process [45–47]. In the total protuberant foil bearing test (both journal and thrust bearings are protuberant type), the gas foil bearing ran smoothly and had high repeatability during the operation of the turbo-expander [15]. Subsequently, the multi-leaf foil journal bearing and thrust bearing were simultaneously applied to this high-speed turbo-expander. Transient acceleration and high-speed operation tests of the rotor bearing system were carried out, and the transient characteristics and the stability of the rotor bearing system were analyzed. The experimental results showed that the starting friction torque and running resistance torque of the multi-leaf foil journal

bearings were small, and the turbo-expander could reach 93,900 rpm. The rotor bearing system had the advantages of a smaller main frequency amplitude (<6 μm) and a low frequency whirl amplitude (<0.5 μm). The rotor locus was clear and had good repeatability [48,49]. For a cryogenic working fluid turboexpander, as shown in the Figure 3, a helium turbo-expander (500 W) supported by foil journal bearings was designed with a rated speed of 220 krpm. The speed could exceed the design speed by 12% in the test with a normal air temperature. In addition, a series of hydrogen turbo-expanders supported by gas foil journal bearings was also developed, with power ranging from more than 20 to 40 kW. The turbo-expander (39.7 kW) with a rated speed of 74.5 krpm was tested in an air environment. The speed could reach 60 krpm, and the amplitude was less than 0.023 μm. Due to their large thrust force, hydrostatic bearings are currently used as thrust bearings in hydrogen and helium turbo-expanders. Meanwhile, a carbon dioxide turbo expander with a rated speed of 125,000 rpm has been designed, as shown in Figure 4. Its rated speed could reach 80,722 rpm in the experiment. The inlet and outlet pressures of the turbo expander were 10 MPa and 7.3 MPa, respectively.

**Figure 3.** Helium turbo-expander and hydrogen turbo-expander by Xi'an Jiaotong University.

**Figure 4.** CO2 turbo-expander by Xi'an Jiaotong University.

#### *2.2. Application in Air Cycle Machine (ACM)*

An air cycle machine is always used to adjust the environment in an aircraft cabin, which is the core component of the air conditioning system. It is also widely used in armored vehicles abroad for their high reliability, long service life, and oil-free characteristics. Rotor systems with gas foil bearings have become popular. At present, the main domestic research institutions engaged in ACM research are the Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University, the 16th Research Institute of China Electronics Technology Group Corporation, Hunan University, and other institutions.

Since 2014, domestic research work has been carried out on the ACM system, and some progress has been made. The Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University developed an airborne turbo cooler supported by gas foil bearings [50]. The structure of airborne turbo cooler is shown in Figure 5. The rotor mass is 900 g and the rated speed is 40 krpm. A series of tests have been carried out on its low temperature repeating start–stop and reliability performances. In vibration tests and extreme temperature (−55 ◦C and 70 ◦C) start-up tests, its performance was better than expectations. In the tests, the coating did not peel off after a hundred start–stop and over-speed operations. In addition, a 13-kW ACM supported by gas foil bearings was designed and developed by Hunan University, and experimental research was conducted on its rotor bearing system [51]. The critical speed and unbalanced response of the ACM rotor bearing system were analyzed, and the structural parameters of the rotor bearing system were modified according to the analyses. The vibration amplitude of the prototype rotor was finally reduced to about 10 μm. In addition, a performance test on the turbo cooler prototype with foil bearings was conducted by Deng [52]. After 10,000 start–stop life tests and vibration tests, it was found that this machine had a good working performance. So far, the total running hours of this type of turbo cooler exceeded 1000 h.

**Figure 5.** Airborne turbo cooler by Xi'an Jiaotong University [50].

#### *2.3. Application in Centrifugal Blowers*

The centrifugal blowers are important in industrial ventilation, which are widely used in many situations, such as medical, chemical, and sewage treatment. Due to the advantages of the high adaptability of the gas foil bearings, the service life and reliability of the blowers are greatly improved under extreme conditions, such as in a dusty environment. On the other hand, blower with gas foil bearings have the merits of low noise, high efficiency, and no pollution. In view of above advantages, centrifugal blowers supported by gas foil bearings have become the first choice for energy-efficient blower. Since the 1970s, foreign blower companies have tried to apply foil bearings in centrifugal blowers. By the end of the last century, highly efficient, energy-saving and highly integrated gas foil bearing blower products have come out on the market. Although domestic theoretical research on gas foil bearings has already been carried out, due to constraints such as materials and process levels, mature products are not yet available on the market. The main domestic research work in this area is the Institute of Machinery Manufacturing Technology, CAEP (China Academy of Engineering Physics) and Shijiazhuang Kingston Bearing Company. Shu X J developed a high-speed air foil bearing with a load capacity of 12.5 kg and tested

the bearing performance [53]. By assembling air foil bearings, high-speed motors, and fan impellers, a prototype centrifugal blower was designed and manufactured. The centrifugal blower is shown in Figure 6. The diameter of the rotor is 80 mm and the maximum speed of the original machine is 30 krpm.

**Figure 6.** Centrifugal blower by the Institute of Machinery Manufacturing Technology, CAEP [53].

#### *2.4. Application in Oil-Free Turbocharger*

In the vehicular industry, an engine's output power can be greatly improved by increasing the air–fuel ratio of a turbocharger. Up to now, oil lubricated bearings have mostly been used in automotive turbochargers. Due to the high working temperature, there will be oil quality deterioration and coking, which may lead to catastrophic failure. These problems can be solved effectively by foil bearings. At the same time, the weight, size and complexity of the devices are greatly reduced. In addition, fuel economy and exhaust emissions are greatly improved. At present, domestic studies on turbochargers have also been carried out.

In order to improve turbocharger performance and solve problems, such as oil leakage, Hunan University conducted theoretical and experimental research on oil-free turbochargers supported by air foil bearings, as shown in Figure 7. Subsequent acceleration and deceleration test showed that the rotors exhibited large sub-synchronous vibrations at 20–60 krpm, but the sub synchronous vibration of the rotor was significantly reduced at 68 krpm. The rotor trajectory and FFT analyses showed that the main vibration of the rotor at a steady speed of 68 krpm originated from the synchronous vibration, and the turbocharger could operate stably at 68 krpm [54,55].

**Figure 7.** Oil-free turbocharger host machine by Hunan University [54].

#### *2.5. Application in High-Speed Gas Compressor*

Recently, high-speed compressors supported by gas foil bearings have become an indispensable part of fuel cell systems. Air compressors are mainly used to provide pressurized clean air. The performance of a fuel cell system is directly related to the pressure of the oxygen provided. When the air pressure is increased, the energy density of the fuel cell system can be improved and a higher output power and power performance can be obtained. Up to now, high speed air compressors with gas foil bearings have been put into use in some renewable energy automotive fuel cell systems, for example, Toyota and Honda in Japan, Garrett in the United States, etc.

After nearly ten years of hard work, prototype compressor products can be produced in China. An air compressor with a maximum speed of 120 krpm (10 kW) was developed by the Institute of Refrigeration and Cryogenics of Xi'an Jiaotong University for automotive fuel cell systems. Subsequently, prototype trial production and preliminary tests were completed. The high-speed air compressor is shown in Figure 8. The test results showed that the air compressor operated stably, and the overall technology reached the leading level in China. In addition, the State Key Laboratory for Strength and Vibration of Mechanical Structures (Xi'an Jiaotong University) developed a 11-kW, 70-krpm high-speed experimental air compressor system supported by gas foil bearings. Experiments including loss, rotor transient temperature, and critical speed were carried out and a stable operation of the system under 70 krpm was realized [56,57]. In addition, Hunan University processed and built an oil-free high-speed air compressor experiment test rig [58]. After many test cycles, a diagram of the rotational speed and the temperature characteristics of foil bearings were drawn. Under different experimental conditions, two time-domain waveforms were used to configure the orbit of the rotor system. Meanwhile, the Dalian University of Technology designed a high-speed two-stage air compressor with gas foil bearings [59,60]. Using DyRoBes-Rotor software, the dynamic characteristics of the rotor bearing system were analyzed. The results showed that the rotor-bearing system performances had good repeatability. The maximum bearing force was 60 N, which was much lower than the bearing capacity of 100 N. For refrigerant compressors supported by foil bearings, a centrifugal compressor (630 kW) was produced by the GREE company with a rotor amplitude of less than 12 μm. The unit energy efficiency coefficient of performance (COP) reached 6.35, and the integrated part load value (IPLV) reached 10.15. Meanwhile, the R134a high-speed refrigerant compressor with a rated speed of 80,000 rpm was developed by Xi'an Jiaotong University, and is shown in Figure 8. Its pressure ratio is 2 × 2, but its size is only one third that of conventional compressors.

**Figure 8.** High-speed compressors and gas bearings by Xi'an Jiaotong University (**a**) air compressor; (**b**) refrigerant compressor.

#### **3. Challenges and Prospects**

#### *3.1. Load Capacity of Gas Foil Bearing*

At present, domestic research on gas foil bearings is mostly focused on gas foil journal bearings, and the research content is mainly on static and dynamic performance, such as bearing capacity and stability. In a rotor bearing system, the foil thrust bearing also plays an equally important role in high-speed turbo machinery. From the perspective of the working principle, there are many similarities between the gas foil journal and trust bearings. Due to the increase in the rotational speed and the power density of turbo machinery, higher requirements are imposed on the bearing capacity and stability of gas foil thrust bearings. Compared with gas foil journal bearings, gas foil thrust bearings require a greater carrying capacity. The bearing capacity and reliability of gas foil thrust bearings have become a bottleneck for turbo machinery. Due to insufficient bearing capacity of the foil thrust bearings, turbo machinery is prone to instability and failure. Research on thrust bearings mainly focuses on experimental research. It is essential to research and develop new high-precision manufacturing and high-performance gas foil thrust bearings. In the development of foil bearing, the operation reliability and the practicability of the manufacturing process also need to be considered. At present, efforts are being made to improve the bearing capacity of journal and thrust foil bearings to 0.6 MPa and 0.4 Mpa, respectively, in China.

#### *3.2. Design Criteria of Gas Foil Bearing*

Due to the diversity of industrial applications, higher requirements are put forward for the variety and performance of high-speed turbo machinery. Therefore, design criteria for gas foil bearings and appropriate foil structure parameters are crucial for their actual application; for example, the trade-off relationship between foil stiffness and Coulomb damping. In some cases, a large Coulomb damping is formed due to the high foil flexibility, which is beneficial for improving the stability of the rotor bearing system and lessening the whirl vibration amplitude. However, the large flexibility will cause a large deformation of the foil and rupture the gas film. In addition, the dynamic analysis method of the foil rotor bearing system also needs to be improved. Efforts are being made to increase the operating DN value of the foil bearing to 4.5 × <sup>10</sup><sup>6</sup> mm · r/min, in China.

#### *3.3. The Stability of Gas Foil Bearing*

The supporting force provided by the gas foil bearings is nonlinear due to the influence of geometric and structural parameters. At the same time, the high-speed rotor system also has relatively complicated dynamic characteristics, which makes stability analysis of the rotor bearing system very complicated. Larger vibration amplitudes are prone to appear in the rotor system. Subsynchronous vibration may cause serious accidents with rotating machinery. At present, a great deal of research has been carried out on subsynchronous vibration of foil bearings. However, the in-depth mechanism of the emergence of subsynchronous vibrations has not yet been found. In order to reveal the influencing mechanism of foil bearing stability, some researches could be carried out from various points of view, such as establishing a complete dynamic solution model considering gas film, foil thickness, and rotor structure, and improving experimental conditions. It is significant to analyze the formation mechanism of the non-linear bearing capacity of the gas foil bearing through experiments and theoretical calculations.

#### *3.4. Thermal Management of Gas Foil Bearing*

During high speed operation, a large amount of heat is generated in the foil bearing due to viscous dissipation or high-temperature working environments in some special cases, such as with micro gas turbines. The heat needs to be dissipated or discharged in time, otherwise it will cause a gas foil bearing failure. Therefore, various measures have been taken to transfer heat to prevent heat accumulation, such as processing multi-layer gap foils, adding holes, and adding heat dissipation channels on the bearing housing that can be used to guide cooling air. Apart from this, the development of phase coatings and supporting materials with a high heat capacity is also very beneficial for gas foil bearing thermal management. At present, a great deal of research on the thermal management of foil bearings has been carried out in China. Efforts are being made to increase the maximum operating temperature of foil bearings to 650 ◦C.

#### *3.5. The Failure Protective Measures and Applications in Extreme Condition*

During the operation of the high-speed turbo machinery supported by gas foil bearings, high temperature failure caused by the accumulation of heat and the transient load in the complex environment will cause a sudden stuck of the rotor-bearing system. If the scratches are severe, it will be necessary to replace the rotor with a new one, which would cause a certain amount of economic loss and wasted time. Therefore, failure mechanisms, failure protection, and preventive measures for gas foil bearings need to be considered, which are of great significance to a stable operation. During the operation of turbo machinery, more attention should be paid to the foil temperature and rotor vibration. The temperature of a foil surface can be monitored using temperature sensors in the prevention of bearing failure due to excessively high temperatures. In addition, rotor vibration is an important parameter reflecting the operating stability of high-speed turbomachinery. The vibration displacement, velocity, or acceleration signals can be obtained by data acquisition system. When the rotor amplitude is too large, turbo machinery should be slowed down in time to prevent sudden instability and shutdown. It is a trend to equip foil bearings in aeronautical and astronautical devices. The bearing capacity of foil bearings are greatly affected by working environment, for instance in rarefied air and extreme temperature situations. At the same time, the reliability of foil bearings is required to be higher for aeronautical equipment. The challenge is how to improve the reliability, service life, and load capacity of the foil bearing. Compared with ball bearings and rolling bearings, the load capacity of foil bearings is relatively small. Due to the influence of manufacturing materials, coating technology, and turbine mechanical thermal management, the reliability and service life of foil bearings could be greatly affected.

#### **4. Conclusions**

The development of applications for gas foil bearings in China was summarized in this paper. Related challenges during the development process were discussed, including analyses of static and dynamic characterizations, thermal characterization, and lubricating coatings. Based on the above review and discussion, the following conclusions could be drawn:

(1) Gas foil bearing technology developed rapidly in China because of application prospects in high-speed turbo machinery. The diversity of gas foil bearing applications also put forward higher requirements on bearing design. Currently, the design criteria of gas foil bearings need to be configured for different application considerations, especially for high-performance gas foil journal and thrust bearings.

(2) Advances in foil coating technology could greatly improve the reliability and service life of turbo machinery. Research on solid lubricating coating of gas foil bearings was closely related to solid lubricating technology and coating technology. Limited by the temperature range of coating materials, friction pair selection and feasibility of manufacturing and processing, further foil coating technology research is thus needed for bearing applications.

(3) One of the major failure modes for gas foil bearing was due to thermal instability. Research is needed to develop suitable predictive tools to address this problem at the design stage. A detailed non-linear thermo-elastohydrodynamic analysis is needed, considering thermal, hydrodynamic, and structural deformation. In the design and application process, different measures should be adopted to reduce heat accumulation and bearing failure during bearing operation.

**Funding:** This project was funded by the National Key R&D Program of China (2018YFB2000100), the Na-tional Natural Science Foundation of China (51976150) and the Youth Innovation Team of Shaanxi Universities.

**Institutional Review Board Statement:** This paper does not involve humans or animals.

**Informed Consent Statement:** This paper does not involve humans or animals.

**Data Availability Statement:** The study did not report any data.

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

#### **References**


## *Article* **Air-Based Contactless Wafer Precision Positioning System**

**Rico Hooijschuur, Niranjan Saikumar, S. Hassan HosseinNia and Ron A. J. van Ostayen \***

Department of Precision and Microsystems Engineering, Delft University of Technology, Mekelweg 5, 2628 CD Delft, The Netherlands; ricohooijschuur@hotmail.com (R.H.); niranjans87@gmail.com (N.S.); S.H.HosseinNiaKani@tudelft.nl (S.H.H.)

**\*** Correspondence: R.A.J.vanOstayen@tudelft.nl; Tel.: +31-15-2781647

**Featured Application: In this paper an improved contactless handling and positioning method for thin substrates is proposed based on air bearing technology. Being fully contactless, the potential application of this work can be found in the high-tech industry where thin, flexible, but at the same time fragile substrates need to be handled with minimal risk of contamination, damage, or even breakage.**

**Abstract:** This paper presents the development of a contactless sensing system and the dynamic evaluation of an air-bearing-based precision wafer positioning system. The contactless positioning stage is a response to the trend seen in the high-tech industry, where the substrates are becoming thinner and larger to reduce the cost and increase the yield. Using contactless handling it is possible to avoid damage and contamination. The system works by floating the substrate on a thin film of air. A viscous traction force is created on the substrate by steering the airflow. A cascaded control design structure has been implemented on the contactless positioning system, where the inner loop controller (ILC) controls the actuator which steers the airflow and the outer loop controller (OLC) controls the position of the substrate by controlling the reference of the ILC. The dynamics of the ILC are evaluated and optimized for the performance of the positioning of the substrate. The vibration disturbances are also handled by the ILC. The bandwidth of the system has been improved to 300 Hz. For the OLC, a linear charge-coupled device has been implemented as a contactless sensor. The performance of the sensing system has been analysed. During control in steady-state, this resulted in a position error of the substrate of 12.9 μm RMS, which is a little more than two times the resolution. The bandwidth of the OLC is approaching 10 Hz.

**Keywords:** air-bearing; precision positioning; tribotronics

#### **1. Introduction**

In the high-tech industry products such as displays, solar panels, and chips are produced. These products are made out of brittle materials where thickness is an important factor. The thickness of the raw material is a factor in the total costs. Additionally, the efficiency of the production process can be increased by increasing the size of the substrates. This is why there was a trend where the substrates increased in size while decreasing in thickness [1]. However, this trend came to a halt due to the percentage of damaged and broken wafers during the manufacturing process.

A possible solution is to handle the products without mechanical contact, floating the substrate on a layer of air. This type of air-based levitation is known from air-bearing technology, where pressurized air forms a thin air film that separates the two surfaces. This concept is an alternative to mechanical contact, where friction, stick-slip, and wear have to be minimized. Furthermore, the layer of air can support the whole surface of the substrate. This is an improvement over the point contact that is currently being used [2], since an air bearing limits particle generation and resulting contamination and forces on the substrate. This will increase the yield from a single wafer.

**Citation:** Hooijschuur, R.; Saikumar, N.; HosseinNia, S.H.; van Ostayen, R.A.J. Air-Based Contactless Wafer Precision Positioning System. *Appl. Sci.* **2021**, *11*, 7588. https://doi.org/ 10.3390/app11167588

Academic Editors: Terenziano Raparelli, Federico Colombo, Andrea Trivella and Luigi Lentini

Received: 25 July 2021 Accepted: 14 August 2021 Published: 18 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Research at the Delft University of Technology has generated multiple concepts where substrates are controlled on a controlled layer of air. The concepts vary in the way the air is controlled, i.e., by variation of the pressure inlet [3], deformation of the surface to control the airflow [4], or variation of the outlet restriction [5]. This technology has also been used to levitate and transport different kinds of materials [6,7], because the air-based levitation is not limited to specific substrate material properties.

When the substrate is floating in a positioning system it is necessary to measure and control its position. In the research cited above, the substrates were modified to measure their position. For example, gratings were added to the wafers. Research has been done on contactless positioning [6,8], but sensing systems that do not require any modifications to the substrate have been limited in either bandwidth or resolution.

This paper presents an air-bearing positioning system where no modifications are made to the substrate, in line with how such an air-bearing motion system is to be used in the high-tech industry. The implementation of a charge-coupled device (CCD) is proposed for noncontact position measurement. The sensing concept is implemented in the contactless actuation system and the results are presented.

In Section 2, the state of the art is presented. Section 3 continues with a focused and detailed description of the existing contactless actuation concept based on deformable surfaces as is relevant for the rest of the paper. The dynamics response of the actuator is evaluated and a controller is designed in Section 4. Section 5 shows the proposed sensor concept with performance analysis. The sensor is implemented in the contactless actuation system and the results are presented in Section 6.

#### **2. State of the Art**

A survey of recent developments in the field of contactless manipulation is presented in [9]. The use of airflow for manipulation and control of the localised deformation of wafers has also been studied in [10], while a pick-up head for wafers has been designed using similar concepts in [11]. Ultra-sound air film technology has also been studied for this application in [12]. However, while this technology provides a viable alternative for contactless pick-up of wafers, it is mainly foreseen to be used in a pick-up head and not for complete manipulation and/or transportation of wafers using airflow.

The start of the development of the new contactless actuation concept at TU Delft was in 2007 [3]. The concepts use the same physics to control the position of the substrate. First, the substrate is levitated on a layer of air. Subsequently, the substrate is preloaded in the out of plane direction with a controlled vacuum source. This increases the out of plane stiffness at the fly height. This height can be adjusted by changing the force the vacuum and pressure source are applying.

The in-plane actuation is achieved by creating an airflow tangential to the substrate surface. The air is flowing from the high-pressure area to the lower pressure following the path of the least resistance. This flow of air exerts a viscous traction force on the substrate. This innovative idea has lead to multiple concepts that will be discussed in this section.

#### *2.1. Variable Inlet Pressure*

In 2007, Wesseling designed a contactless positioning system based on a fixed geometry and a variable inlet pressure to control the thin film airflow in his PhD thesis "Contactless Positioning Using an Active Air-Film" [3]. Figure 1a shows a cross-section of one of the actuator cells, where a substrate is floating on top of the pressurized air.

The air will flow from the inlet *p+* with high pressure to the outlet *p-*, where the air is removed from the cell at a sub-ambient pressure. The viscous shear of the airflow creates a traction force on the substrate. The direction of the flow within a cell is controlled by changing the pressure at the different inlets. Four actuator cells have to work together to actuate the substrate in three degrees of freedom. The flow of the different cells is separated by dams. Four actuator cells are shown in Figure 1b.

The final system was used to position a 4-inch silicon wafer, a demonstrator as shown in Figure 1c. A bandwidth of 50 Hz was realized using a PID controller. The system had a maximum acceleration of 600 mm/s2, with an error of 6 nm (1*σ*), while using an external active vibration isolation system. The concept's performance was limited by the manufacturing tolerances of the system, since this placed a lower limit on the fly height of 15 μm. A lower fly height would result in a higher actuation force. The other limiting factor was the time delay caused by the pressure lines between the control valve and the pressure inlet in the cell.

**Figure 1.** Working principle of the first generation. (**a**) Shows the cross-section of an actuator cell. (**b**) Top view. (**c**) The pressure variation demonstrator. For reference, the size of 6 × 6 actuator array in the centre of the image is 60 mm × 60 mm [3,13].

#### *2.2. Variable Outlet Restriction*

In 2017, Verbruggen designed the third generation of the contactless actuator which uses an outlet restriction to create the traction force. The dams can be actuated to restrict the flow at one side, while the flow is opened up at the other side. The concept is similar to the variable inlet pressure design. The advantage of this system is that the restriction is close to the chamber. This is an improvement, because there is no delay present at the supply lines, which limited the variable inlet pressure stage.

A demonstrator was built to prove the concept. Seven circular pockets have been produced. The substrate was controlled in open-loop. Using feedforward average angular accelerations have been measured of 11 rad/s2.

#### *2.3. Transportation*

Besides positioning a substrate on a stage, this technology can also be used for transportation thus further reducing the need for mechanical contact in a typical production system. In 2016, Vagher investigated the construction of a passive contactless conveyor [6]. In this research, a cost-effective manufacturing technique for the passive conveyor surface was explored.

Vagher's conveyor concept moves substrates in 2D and measured the position of the substrate using an automated vision system. The vision system resulted in a resolution of 30 μm at 60 Hz. In continuation, Snieder investigated the transition between active and passive surfaces, which makes the transportation concept more versatile [7].

#### *2.4. Deformable Surface*

In 2011, Vuong explored a new concept regarding contactless positioning, "Air-based contactless actuation system for thin substrates" [4]. The concept is based on multiple tilted surfaces to steer the airflow, and therefore the force on the substrate. The advantage of the

tilting surface is that a change in the geometry directly generates a change of the shear force on the substrate. This concept does not suffer from the delay between the controlling action and force as seen in the first generation. In this concept, the in- and outlet pressure can be kept constant. The concept where the surfaces are tilted by bending the shafts is shown in Figure 2. This concept is used in this paper for the precision contactless positioning of the wafer and is explained in detail in the next section.

**Figure 2.** The tilting surface concept, the surfaces are titled by bending the shafts [4].

#### **3. Flowerbed**

A demonstrator was built based on the bendable shaft concept. However, instead of bendable shafts, two membranes are used which allows the shafts to tilt. This reduces the force requirement of the actuator. The concept is shown in Figure 3. The top membrane connects the flowers to the world. The bottom membrane is attached to the movable plate, which is connected to an actuator. The top membrane also functions as a seal for the upper chamber, which is connected to the vacuum pressure outlet. The inlet pressure comes in through the hollow centres of the shafts.

**Figure 3.** Schematic of the second generation contactless actuator, showing the flow of air through the system, for a displacement *d* of the movable plate [13].

Figure 4 shows a cut-through view of the Flowerbed. The Flowerbed has 61 hexagonal heads, with an in-radius of 14 mm. It is nicknamed the Flowerbed because the surface looks like flowers that tilt towards the sun. The hexagon shape of the flowers was selected because it provides a full packing of the surface in combination with an optimal circumference to surface area ratio, and thus a maximum force for a given airflow. Furthermore, the figure also shows the fixed plate and movable plate, where both consist of two plates that are clamped together with a 50 μm thick membrane in between. The membrane is also clamped in the flowers. To better align the flowers, the movable plate is preloaded downwards with three springs, with a force of approximately 100 N. After assembly, the final step was the spark erosion of the surface with wire-cut EDM technology to further level the surface.

**Figure 4.** A cross-section of the Flowerbed showing the fixed plate, flowers, and movable plate. Source: [13].

Figure 5 shows a single flower in both a neutral and actuated position. The thin film airflow is visualized in the area between the effective surface of the flower and the substrate. The combination of the positive and negative pressures provides stiffness in the out-of-plane direction at the fly-height *h*.

**Figure 5.** Schematic of a single flower of the Flowerbed in (**a**) neutral position and (**b**) actuated position. The force on the wafer is visualized in the actuated position [13].

The flowers are tilted with an angle *α*, through a displacement of the movable plate, as shown in Equation (1). Where *dp* is the displacement of the movable plate and L is the distance between the movable plate and the fixed plate.

$$
\sin(\alpha) = \frac{d\_p}{L} \approx \alpha \tag{1}
$$

when tilted the majority of the airflow will follow the path of least resistance. This flow will, in turn, exert a viscous traction force *Fw* on the wafer in addition to the out of plane stiffness.

The tilting angle of the flower has been measured to be proportional to the displacement of the movable plate up to 1000 Hz, without delay. The tilt response between the individual flowers is similar, there is only a small deviation in the gain. The values of gain vary between 71 rad/m and 80 rad/m, where the predicted value is 77 rad/m [4]. Furthermore, Vuong experimentally validated that the pneumatics can be modelled by a proportional gain, up to at least 400 Hz. This means that if the movable plate position is known, the force exerted on the wafer is known.

The stiffness of the Flowerbed is different for the translation and rotation of the movable plate. Each flower has an elastic hinge at the top and bottom, which provide stiffness against rotation. Figure 6, shows the behaviour of the Flowerbed during the translation and rotation of the movable plate. The stiffness and inertia values for translation and rotation result in two suspension eigenvalues. These values have been calculated and experimentally validated by Jansen [13]. The important properties are summarized in Table 1.

**Figure 6.** Schematic of the flowerbed during (**a**) translation and (**b**) totation [13].



#### *3.1. Flowerbed Actuation*

In 2016, Krijnen designed a reluctance actuator to actuate the movable plate of the Flowerbed [14,15]. Reluctance actuators can provide a high force, but they also have a high mass. With the use of fractional order control, this resulted in a control bandwidth for the Flowerbed of 100 Hz. The bandwidth was limited by the eigenmodes caused by the actuator as shown in Figure 7. Figure 8 shows the flapping mode of the actuators that was limiting the bandwidth. The position of the wafer is measured using optical sensors. This required that gratings are attached to the wafer. The position of the wafer is controlled with a bandwidth of 60 Hz.

**Figure 7.** The simulated and measured frequency response function of the Flowerbed with the reluctance actuator [14].

**Figure 8.** Flapping modeshape of the actuators at 1080 Hz (measured) [14].

To increase the bandwidth of the system, Jansen redesigned the actuator in 2018. The new actuator was designed with special care for dynamic performance. Instead of traditional hinges, a design with notch flexures is chosen. The flexures provide a high stiffness in the out of the plane direction and work without backlash. The reluctance actuators have been replaced with Lorentz actuators, because they have a linear force current relation. Furthermore, the permanent magnets were placed on the stationary frame, and the coils on the mover, resulting in a minimal moving mass.

Figure 9 shows the chosen guidance for the movable plate. In the figure, the application point of the force and sensor are shown, as well as the notch flexure and guide linkages. The sensors that measure the actuator displacement are placed close to the actuator.

**Figure 9.** Compliant mechanism design used to house the Lorentz actuators. Adapted from [13].

Figure 10 shows the Flowerbed including the realized actuator. From top to bottom it shows the connection between the actuator and movable plate. The coil is used to generate the force on the movable plate. The sensor is placed directly behind the actuator. The actuators have a relatively small motion range thus the wires running from the coil to the base plane are no problem. Furthermore, the stiffness is small compared to the stiffness of the movable plate [13].

**Figure 10.** Photograph of the Flowerbed with the compliant mechanism and Lorentz actuator design [13].

The position of the actuators is measured using Micro-Epsilon CapaNCDT CS08 capacitive sensors, which have a range of 800 μm. With the 16 bits ADC this resulted in a resolution of 24.4 nm. The noise is measured at 0.6 μm peak-to-peak, with an RMS value of 72 nm.

#### **4. Flowerbed Dynamics and Control**

In the previous sections, the state of the art has been presented. The bandwidth of the compliant actuator design was placed on 233 Hz. No sensor was implemented to measure the position of the wafer. In this section, the dynamic response of the system is presented. A controller is designed and the closed-loop results are discussed. At last, the possibilities for improvement are discussed.

#### *4.1. Flowerbed Dynamics*

Figures 11 and 12 show the results from the identification where all the straight transfers and cross transfers are plotted. The directions of force and position sensor are inverted and account for the 180◦ phase seen. Both transfers show a quick drop in phase which is the limiting factor for the bandwidth. Jansen suggested this could be the effect of the capacitive sensing system [13]. This large phase drop was investigated as part of this work. After investigation of the delay sources in the system— amplifier, capacitive sensors, sampling time, and coil dynamics—the main sources of the delay have been identified. The phase lag is due to a combination of the 20 kHz sampling which produces a phase lag and the magnetic field which is lagging the current. The lagging field is due to the eddy currents in the mover which will create a magnetic field opposing the movement. When the magnetic field of the coil and eddy current are summed the field lags that of the coil alone. This provides a time delay seen in the bode diagram as a lagging phase. Furthermore, the eddy current damping provides damping of magnitude proportional to the velocity of the mover.

**Figure 11.** Frequency response function from actuator force to actuator displacement, the phase is shifted upwards by 180 degrees since the position is measured in opposite direction to applied force.

**Figure 12.** The cross-coupled frequency response functions of the Flowerbed, the phase is shifted upwards by 180 degrees since the position is measured in opposite direction to applied force.

Figure 13 shows the mover in the magnetic field of the permanent magnets. The mover consists of an aluminium plate with spaces cut out for the coils and a top plate to hold them in place. A solution to limit the intensity of the generated eddy currents is to laminate the mover. Another method is to look at the possibility of flux feedback for the amplifier of the system, to compensate for the effects of the eddy current [16]. These solutions have not been implemented in this work and are left for future work.

**Figure 13.** Simulation showing the magnetic flux density (in T on the colorbar) and the eddy current density (arrows) [13].

#### *4.2. Flowerbed Control*

The controller of the inner loop will control the position of the actuator to the corresponding sensor. The measured position will be controlled using a feedforward controller and feedback controller. The diagram is shown in Figure 14.

**Figure 14.** Movable plate control strategy (ILC).

For the feedback controller, a PID controller is designed for a bandwidth of 300 Hz. The parameters are shown in Table 2, the gain for actuator 2 is different to compensate for the slightly different stiffness at the bandwidth. The controller creates a phase margin of 27◦, gain margin of 32.7 dB, and a modulus margin of 11.3 dB.

**Table 2.** ILC-PID parameters.


The feedforward controller has been designed using a simplification of the plant. The feedforward controller is the inverse of this plant made proper with two second-order low pass filters at the frequency of the bandwidth. The closed-loop response of the system with the PID controller and feedforward is shown in Figures 15 and 16. The feedforward controller is given in Equation (2).

$$FF = \frac{\left(\frac{s}{2\pi 300} + 1\right)\left(\frac{s^2}{(2\pi 49)^2} + \frac{0.1s}{2\pi 49} + 1\right)}{0.37\left(\frac{s^2}{2\pi 300} + \frac{2s}{2\pi 300} + 1\right)^2} \cdot \frac{\left(\frac{s^2}{(2\pi 75)^2} + \frac{0.17s}{2\pi 75} + 1\right)}{\left(\frac{s^2}{(2\pi 60)^2} + \frac{2s}{2\pi 60} + 1\right)}\tag{2}$$

**Figure 15.** Closed loop frequency response of the actuators, from actuator reference position to measured position.

**Figure 16.** Closed loop cross coupled frequency response of the actuators, from actuator reference position to measured position.

The closed-loop response of the system is similar for all three actuators. The response provides a good basis for the outer loop controller up to, at least 100 Hz. The system shows strong coupling at 470 Hz, this could be improved using a decoupling filter or by prefiltering the input. Because the bandwidth of the outer loop is limited, it is chosen to use the outer loop controller to prefilter the reference signal.

#### **5. Wafer Sensing**

The wafer floating on air is a floating mass with a viscous damping system, thus to gain control over the position of the wafer it is necessary to have a collocated sensor. For this sensor in the outer loop, there are some restrictions, that influence the sensor choice. For the high-tech industry, it is unappealing to use a wafer modified to accommodate the sensing solution. Thus, it is preferred that the wafer is not altered.

As discussed in the state of the art, different sensor techniques have been used in the past, such as techniques based on capacitive sensors, edge detection using charge-coupled devices (CCD), and vision and optical encoders. The best speed and resolution was achieved using optical encoders with gratings on the wafer. However, this is an alteration to the wafer. Without altering the wafer edge detection using the linear CCD sensor array is the best choice. Edge detection can reach sampling speeds up to multiple kHz while having a sub-micrometre resolution that can be achieved using sub-pixel interpolation. This will be discussed later in this section.

#### *5.1. Linear Sensor Array*

A linear CCD has been selected to measure the edge position of the wafer. It consists of an array of pixels that measure the intensity of the light. The change from light to dark is used to find the edge of the wafer. The concept is shown in Figure 17. As discussed in the state of the art, Wesselingh had implemented CCD sensors but used a somewhat different concept where noise from stray light limited the performance [3].

The selection of the CCD sensor is based on a trade-off between the resolution, sampling time, and range. The sampling time is limited by the available modules for the National Instruments CompactRio. The NI9201 is used which supports 500 kS/s. The TCD1103GFG has 1500 pixels of 5.5 μm by 64 μm on a 5.5 μm center. It supports a data rate of 2 MHz. With 1500 pixels and a low integration time this results in a sampling rate of almost 1.3 kHz, and with the limitation of the 500 kS/s this is lowered to 320 Hz. High power cool white LEDs are used as the light source. The intensity of the light is tuned to maximize the dynamic range.

**Figure 17.** Linear CCD sensor concept to measure the edge of a wafer.

In this paper, the sensor output is compared to a threshold. This threshold is used to detect the transition from light to dark. Instead, a threshold a sub-pixel interpolation method could also be used. Using this method it is possible to increase the resolution of the sensors to better than 0.05 of the pixel size, while at the same time reducing the effects of noise [17–20]. Using the sub-pixel interpolation method the theoretic resolution of the chosen CCD could be as low as 0.11 μm. However, to limit the calculation time on the FPGA, it was chosen to implement the threshold model. Figure 18 shows the output of the CCD sensor. The threshold is placed at 0.6. When the threshold is reached the pixel number corresponds to the position on the CCD where the shadow ends.

**Figure 18.** Output of the CCD sensor. It shows the edge of the wafer. The edge is a curve over 10 pixels.

When using a point light source, the measured position should be adjusted to get the real position. Figure 19 shows that the distance between the shadow and the wafer position is dependent on the position and height of the wafer. The wafer position is calculated using Equation (3a). Because the height of the wafer is not fixed a disturbance could enter the system through a variation in the height. The size of the disturbance can be estimated by Equation (3a).

The difference in fly-height is estimated to be in the order of micrometers. When minimizing the height between the wafer and CCD the disturbance is estimated to be within 0.1 of the pixel size for large movements, which is acceptable when using the threshold model, if sub-pixel interpolation is used it is better to remove the height by using a collimated light source or by adding another LED. With two LEDs, another data point is known which can be used to eliminate the height from the equation as shown in Equation (3b) [21].

**Figure 19.** A CCD sensor with 2 LEDs, showing the actual position of the wafer *xw*, the shadow created by the first led *xs*<sup>1</sup> and second led *xs*2.

$$X\_w = \frac{h\_{wfer}}{h\_{lead}}(x\_{L1} - xs1) + s1\tag{3a}$$

$$X\_w = \frac{\mathbf{x}\_{L2}\mathbf{x}\_{s1} - \mathbf{x}\_{L1}\mathbf{x}\_{s2}}{-\mathbf{x}\_{L1} + \mathbf{x}\_{L2} + \mathbf{x}\_{s1} - \mathbf{x}\_{s2}}\tag{3b}$$

#### *5.2. Performance Analysis*

Before the implementation of the CCD sensor inside the contactless actuator, a performance analysis experiment was done. In these tests, the linearity of the sensor was tested, as well as the steady-state behaviour. The performance test was performed in an environment where disturbance sources from the external lighting, vibrations, and the electrical side were not removed. The steady-state response shows no drift and high precision, since the position returns to the same position. Figure 20 shows the performance analysis test where the position of the CCD is compared to a stage that has a resolution of 100 nm. The figure shows no drift and follows the position of the stage with an error of 3.0 μm RMS.

**Figure 20.** Performance of CCD sensor analysed for linearity. RMS error value 3.0 μm.

#### *5.3. Metrology*

The kinematics are used to convert the measured position from the sensors into the position of the wafer. Figure 21, shows the placement of the sensors. The kinematic equations are shown in Equation (4), where *r* is the radius of the wafer and 41.5 mm the distance between the sensors. The input of the kinematic equations is in millimetres and the output is in millimetres and radians. The photograph of the Flowerbed system with the integrated contactless sensing is shown in Figure 22.

**Figure 21.** The sensor placement in reference to the wafer to be controlled.

$$y = \frac{d\_2 + d\_3}{2} \tag{4a}$$

$$\mathbf{x} = d\_1 - \left( r - \sqrt{r^2 - y^2} \right) \tag{4b}$$

$$Rz = \frac{d\_3 - d\_2}{41.5} \tag{4c}$$

The validation of the outer loop control was performed using a 200 mm dummy wafer. The dummy wafer was cut out of a 750 μm thick aluminium sheet according to the JEIDA standard. This means a circle with a radius of 100 mm with a flat of 57.5 mm. The aluminium wafer had a mass of 63 g, which is 8 percent more than a comparable 8 inch silicon wafer.

**Figure 22.** Photograph of the Flowerbed system with the implemented contactless sensing system without the wafer. The position of the sensors has been annotated according to Figure 21, and for reference, the diameter for the inscribed circle in each of the 61 hexagonal flowers is 14 mm.

#### **6. Wafer Control**

To control the position of the wafer, a second controller needs to be added to the inner loop Flowerbed control. This leads to a cascaded control structure as shown in Figure 23. The Flowerbed controller will be called the inner loop controller (ILC) from now on. The ILC will provide a stiff connection between the Flowerbed and actuator. The setpoint of the ILC will be controlled by the wafer controller, which is also called the outer loop controller (OLC). The kinematics and inverse kinematics are placed in the loop of the OLC, because the OLC has a lower bandwidth and thus more time for the calculations.

**Figure 23.** The control strategy for the flowerbed.

#### *6.1. Pressure Control*

For the operation of the Flowerbed, constant pressure and vacuum are crucial, even more so because of the non-uniform behaviour of the flowers on the surface for a given supply and vacuum pressure. As a result, setting a higher pressure will cause an extra force in the negative y-direction, whereas a lower vacuum will create an extra force in the positive y-direction.

To control the pressure, sensors and valves are placed in the supply lines. The dynamic behaviour of the pressure and vacuum supply has been identified around the operation pressure of 240 kPa supply and 6 kPa vacuum. The controller design is a proportional gain together with a second-order low pass filter with the values provided in Table 3. The lower bandwidth has been chosen because of the limitations of the solenoid valves, increasing the bandwidth further does not improve the performance. The closed-loop frequency responses for the pressure controls are shown in Figure 24. During control the steady-state error is measured, the error for the vacuum controller is 0.01 kPa(*σ*) and for the pressure controller 0.3 kPa(*σ*).

**Table 3.** Pressure controller properties.

**Figure 24.** The closed loop frequency response of the Flowerbed pressure control.

#### *6.2. Pneumatics*

Previous experiments by Vuong showed that the air film dynamics can be described as a pure proportional gain up to at least 400 Hz, after which the measurement became unreliable [4]. The proportional gain of the system was determined at 1.4 mN/μm. In short, this means that if the position of the movable plate is known, the viscous traction force on the wafer can be calculated up to at least 400 Hz. This value of the proportional gain has been verified during the experiments of Krijnen [14]. Furthermore, Krijnen provided the displacement force graph for different supply pressures while maintaining the vacuum supply at ambient pressure minus 7 kPa.

Multiple experiments were done as part of this work with the values provided by Krijnen, however, the previously measured pneumatic gain could not be replicated again without mechanical contact. Multiple factors play a role here, the Flowerbed has been taken apart and reassembled multiple times to replace the actuators. Upon close inspection, it was noticed that multiple flowers were out of alignment. Furthermore, the supply pressure was leaking between the pressure sensor and the Flowerbed. This made the readings from the pressure sensor unreliable. The combination of these factors resulted in the fact that the noncontact criteria could not be guaranteed for a given pressure. Therefore, in this work, the pressure was manually tuned to obtain the best performance while maintaining the no-contact condition for a full movement of 60 μm for the movable plate. The optimum was found at a vacuum pressure at minus 6 kPa and a supply pressure of 245 kPa. As a result, the value of the pneumatic gain was less than optimal and was approximated at 0.3 mN/μm.

#### *6.3. Outer Loop Controller*

In the outer loop, two separate controllers have been designed, one for controlling the rotation at lower bandwidth and one for positioning the wafer in *x* and *y* at a higher bandwidth to maximize the performance of the total system. The parameters of the controller will be discussed later in this section.

The outer loop could not be identified using the open-loop response because the system consists of a floating mass with viscous damping. The identification of the system could be done by adding stiffness to the mass or through simulation. The simulated openloop transfer function could be described by the transmissibility of the inner loop controller (*TILC*), the pneumatic gain of the system and the wafer, which is a mass system. This is shown in Equation (5), where *COLC* is the outer loop controller.

$$PC = \frac{T\_{ILC}C\_{OLC}}{m\_{water}s^2} \tag{5}$$

For the controller, it is decided to use a PD controller with a low pass filter. The low pass filter must be placed before 300 Hz to cancel out the coupling in the inner loop. Because of the high gain of the mass-line, it was chosen to leave out the integrator.

$$K\_p = \frac{m\_{water}\omega\_{bw}^2}{k\_{pw} \ast 4};\tag{6}$$

The gain of the system was calculated using Equation (6). The factor 4 accounts for the frequency range of the derivative part of the controller as given in Table 4, which adds a gain of the same value. To start, the pneumatic gain *kpn* values calculated by Krijnen were used. This gives an initial place to start after which the controller has been tuned iteratively. After the initial measurement, the new controller gain could be calculated to compensate for the lower pneumatic gain.


**Table 4.** Outer loop controller properties.

#### *6.4. Performance*

After tuning the controller for the inner loop, outer loop, and pressure, the performance of the positioning of the wafer has been measured. Figure 25, shows the closed-loop response of the system. The measurements for the cross-coupling between the x and y-axis are insignificant and thus not presented in the frequency domain. The error of the x and y-axis during tracking is presented in Figure 26. Furthermore, when tracking a sine wave with an amplitude of 0.33 mm at 0.5 Hz, the RMS error of x and y remained under 9 μm.

The error at steady-state is shown in Figure 27 with the RMS error values provided in Table 5. The periodic movement is caused by the correlation between the rotation and translational degrees of freedom. This is expected because the dummy wafer is not perfectly round. It must be noted that less rotation was observed during tracking than under steady-state.

**Figure 25.** The closed loop frequency response of the Flowerbed—T from the reference position to the measured position for the translations.

**Figure 26.** Tracking response of the wafer.

**Figure 27.** The steady state error during control of the wafer in x and y.



#### **7. Conclusions**

The push for thinner and larger wafers has substantially increased the need for contactless handling of the substrates. The Flowerbed and other concepts were developed for this purpose using an air-bearing concept. This paper has presented an air-bearing contactless wafer handling system where the wafer's position is sensed without any alterations to the wafer.

The performance of the inner loop has been improved from 233 Hz to 300 Hz by improving the controller. This is an improvement of almost 30 percent compared to previously published results [13]. The inner loop creates a stable platform for the outer loop controller which controls the position of the wafer. External disturbances are handled in the inner loop, since the air does not transmit the vibrations from the floor.

With the use of three charge-coupled devices (CCD) the position of the wafer could be measured by detecting the edge of the wafer. This is done without any alteration of the wafer. The CCD sensors have a resolution of 5.5 μm.

Due to the unforeseen lowering of the pneumatic gain than attained in previous studies, the bandwidth for control has not been improved. However, the wafer has been controlled using the new sensor's design, resulting in a bandwidth of 10 Hz. The wafer can be positioned with a high accuracy of 13 μm RMS, which is under 3 resolution counts of the sensor. This showcases the high precision positioning capability of the Flowerbed.

Further improvement of the inner loop could be achieved without redesign of the actuator. This improvement can be made by laminating the mover such that the effect from the eddy currents is removed. If lamination is not an option switching to flux feedback could remove the effect from the eddy currents, which are lowering the phase and thus the bandwidth. Additionally, decoupling filters could be used to improve the closed-loop sensitivity function of the inner loop. However, even in the current state the inner loop is no longer limiting the performance of the outer loop.

The resolution of the CCD could be improved to sub-micrometre resolution using subpixel interpolation techniques. Additionally, using improved data acquisition hardware it is possible to increase the sampling frequency of the current CCD from 320 Hz to 1.2 kHz.

**Author Contributions:** Conceptualization, all; methodology, R.H., N.S. and S.H.H.; software, R.H. and N.S.; validation, R.H., N.S., S.H.H. and R.A.J.v.O.; investigation, R.H.; writing—original draft preparation, R.H. and N.S.; writing—review and editing, all; supervision, R.A.J.v.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


## *Article* **Rotation Accuracy Analysis of Aerostatic Spindle Considering Shaft's Roundness and Cylindricity**

**Guoqing Zhang 1,2, Jianming Zheng 1,\*, Hechun Yu 2, Renfeng Zhao 1, Weichao Shi <sup>1</sup> and Jin Wang <sup>2</sup>**


**Abstract:** The rotation accuracy of the aerostatic spindle can easily be affected by shaft shape errors due to the small gas film clearance. Thus, the main shaft shape errors with the largest scale—that is, the roundness and cylindricity errors—are studied in this paper, and a dynamic mathematical model is established with the consideration of the roundness, cylindricity errors, and spindle speed. In order to construct the shaft model, the discrete coefficient index of the shaft radius based on roundness measurement data are proposed. Then, the simulation calculations are conducted based on the measured cylindricity data and the constructed shaft model. The calculation results are compared with the spindle rotation accuracy measured using the spindle error analyzer. The results show that the shaft with a low discrete coefficient is subjected to less unbalanced force and smaller rotation errors, as obtained by the experiment.

**Keywords:** rotation accuracy; aerostatic spindle; roundness; cylindricity; dispersion coefficient

#### **1. Introduction**

In contrast with traditional spindles, an aerostatic spindle uses compressed gas to form a gas film, which has the characteristics of high precision, high speed, and low noise. The film thickness of an aerostatic spindle is usually about 10 micrometers, which meets the high requirements regarding manufacturing errors relating to the shaft [1]. Manufacturing errors can be classified as dimensional errors and shape errors, which include roundness, cylindricity, waviness, and roughness. Shape errors of the shaft will cause the uneven distribution of the gas film in the bearing. When the shaft rotates, the thickness of the gas film in the bearing will change continuously, which will cause pressure fluctuation, and then produce an unbalanced force on the shaft. The unbalanced force will affect the stability and rotation accuracy of the spindle.

Pande et al. [2] pointed out that the shape errors of journal bearings can generally reduce the load-carrying ability and increase the flow rate. Song et al. [3] calculated the half frequency whirl phenomena and the instability threshold speeds for the bearings with shape errors. They found that the precision of the shaft is much more important than that of the bearing. Sun et al. [4] studied the effects of shape errors of a shaft-bearing system on dynamic characteristics, and results show that the rotating speed at which the fluid whip occurs increased when shape errors existed. In order to reduce the radial error motion, Cappa et al. [5] analyzed the influence of several manufacturing errors, bearing parameters and feeding geometries of an aerostatic journal bearing and suggested that the radial error motion is mainly influenced by the shape errors of the shaft. Cui et al. [6] presented research on the manufacturing errors of aerostatic porous journal bearings, and found that the circumferential waviness errors caused the obvious inhomogeneity of the flow field and the transformation of the morphology of the high-pressure region. Zhang et al. [7] established a new approximate model to predict the radial error motion of

**Citation:** Zhang, G.; Zheng, J.; Yu, H.; Zhao, R.; Shi, W.; Wang, J. Rotation Accuracy Analysis of Aerostatic Spindle Considering Shaft's Roundness and Cylindricity. *Appl. Sci.* **2021**, *11*, 7912. https://doi.org/ 10.3390/app11177912

Academic Editors: Terenziano Raparelli, Federico Colombo, Andrea Trivella and Luigi Lentini

Received: 27 July 2021 Accepted: 25 August 2021 Published: 27 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

hydrostatic journal bearings, and pointed out that the main working harmonic components in roundness errors had a major influence over rotation accuracy. Wang et al. [8] proposed a numerical model using linear perturbation theory for the research on the influence of journal rotation and bearing surface waviness on the dynamic performance of aerostatic journal bearings. Lee et al. [9] examined the influence of the waviness errors of a hydrostatic journal bearing by considering the rotational and bearing installation angles, and revealed that the load-carrying capacity varied according to the amplitude and phase angle of the waviness functions. Li et al. [10] investigated the effects of different forms of surface waviness on the stability of the hydrodynamic journal bearing system by using non-linear dynamic analysis and calculation method, and the results show that the existence of the surface waviness causes the oil film of the system to wave and increases the energy loss of the system. In order to study the surface roughness effects in hydrodynamic bearings, Quiñonez et al. [11] developed an analytical solution method based on the Reynolds equation using perturbation techniques under the assumption of waviness coupled with linear superposition that can be used in wide exponential land slider bearings. Lin [12] analyzed the effect of surface roughness on the dynamic stiffness and damping characteristics of hydrostatic thrust bearings; according to the results, the mean stiffness and damping behaviours are significantly affected by the roughness pattern and the height of the roughness. Kumar et al. [13] studied the effect of stochastic roughness on performance of a Rayleigh step bearing operating under thermo-elastohydrodynamic lubrication. The results show that the directional orientation of the surface roughness could affect the performance of a bearing significantly. Based on the stochastic method and the Ng–Pan turbulent model, Zhu et al. [14] calculated the turbulent lubrication performance of a journal bearing with an isotropic rough surface. Kim et al. [15] proposed an approach with both stochastic and contact characteristics to evaluate the effects of surface roughness for a slider bearing operating under partial or boundary lubrication, and found that the influence of roughness parameters on friction and load capacity increased rapidly upon the application of asperities contact. Maharshi et al. [16] proposed a stochastic analysis of hydrodynamic journal bearings including the effect of surface roughness and development of the efficient radial basis function based uncertainty quantification algorithm. By using Galerkin's technique to solve the Reynolds equation, Rajput [17] indicated that the performance of a journal bearing system is significantly deteriorated due to the geometric imperfections.

In the above references, the waviness and roughness in the shaft shape errors are mostly studied. However, the roundness and cylindricity are the largest errors among the shaft shape errors, usually between 1 and 4 micrometers. These errors have the most obvious impact on the aerostatic spindle's rotation accuracy. Although Zhang et al. [7] studied the radial error motion of the bearing in a quasi-static state considering roundness errors of the shaft, the rotation accuracy of the aerostatic spindle will change with rotation speed, which is not considered in the present study. Therefore, it is necessary to study the dynamic rotation accuracy of the aerostatic spindle under the effects of roundness and cylindricity errors.

Therefore, a dynamic mathematical model considering roundness, cylindricity errors and spindle speed is established by using the finite difference method based on the Reynolds equation in this paper. The roundness and cylindricity errors of the shaft are measured by Taylor Hobson 585 LT cylindricity meter. The measured data are resampled and filtered and then brought into the model for calculation. Through the comparison of the simulation and the experiment, it is proved that the shaft roundness and cylindricity errors will produce an unbalanced gas film force on the shaft and affect the spindle's rotation accuracy. Based on the measurement data of roundness error, the evaluation index of the dispersion coefficient is proposed. The research shows that this discrete coefficient is a more effective index to predict the spindle's rotation accuracy compared with the roundness and cylindricity.

#### **2. Mathematical Model**

*2.1. Modeling of the Aerostatic Spindle*

The schematic illustration of the aerostatic journal bearing with shaft shape errors is shown in Figure 1. This bearing has two column orifice restrictors at the front and rear, and each column includes eight orifices. The compressed gas flows into the clearance between the bearing and the shaft through the orifice, and then flows out from both ends of the bearing.

Due to the shape errors of shaft surface, the film thickness of different positions in the circumferential and axial directions of the bearing is different. In Figure 1b, *O*<sup>0</sup> is the bearing center, *O*<sup>1</sup> is the shaft center, and *R* is the bearing radius. *h* is the film thickness at point *A* on the shaft surface that can be denoted as

$$h = R - \left| \overrightarrow{O\_0O\_1} + \overrightarrow{O\_1A} \right| \tag{1}$$

When *<sup>θ</sup>* is the angle between −−→*O*0*<sup>A</sup>* and the *<sup>x</sup>*-axis, −−−→ *<sup>O</sup>*0*O*<sup>1</sup> = (*ex*,*ey*), and −−→*O*1*<sup>A</sup>* <sup>=</sup> (*cxθ*, *cy<sup>θ</sup>* ), the film thickness at *θ* is

$$\ln \theta = R - \sqrt{\left(c\_{x\theta} + c\_x\right)^2 + \left(c\_{y\theta} + c\_y\right)^2} \tag{2}$$

If the shaft is rotating clockwise at a constant angular velocity *ω*, then at *t* time,

$$\begin{cases} \begin{aligned} \theta' &= \theta + \omega t + \alpha - 2N\pi \\ c\_{x\theta} &= x\_{\theta'} \\ c\_{y\theta} &= y\_{\theta'} \end{aligned} \end{cases} \tag{3}$$

where *α* is the initial angle of the shaft, *N* is the integer number of rotations of the shaft, and the values of *x<sup>θ</sup>* and *y<sup>θ</sup>* come from the measurement results of the cylindricity meter.

**Figure 1.** Illustration of the aerostatic journal bearing with shaft shape errors: (**a**) axial direction; (**b**) circumferential direction.

By combining Equations (2) and (3), the film thickness at any angle at any time can be obtained.

In order to study the pressure distribution of the gas film in the bearing, the circumferential cylindrical coordinates system is used to reconstruct the journal bearing. Figure 2 is the partial schematic illustration of aerostatic journal bearing in circumference cylindrical coordinates. The circumferential direction of bearing's inner surface is the *x*-axis and the *y*-axis points and is perpendicular to the bearing axis, and the *z*-axis is parallel to the bearing axis.

This spindle system incorporates the following assumptions:


**Figure 2.** Circumference cylindrical coordinates.

The modified Reynolds equation for this type of aerostatic journal bearing [18] is

$$\frac{\partial}{\partial \mathbf{x}} \left( ph^3 \frac{\partial p}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{z}} \left( ph^3 \frac{\partial p}{\partial \mathbf{z}} \right) = 6 \eta u \frac{\partial (ph)}{\partial \mathbf{x}} + 12 \eta \frac{\partial (ph)}{\partial t} + 12 \eta p V\_{\text{in}} \tag{4}$$

where *p* is the gas film pressure, *η* is the dynamic viscosity of the gas, *u* is the velocity of the gas in *x* direction, and *Vin* is the injection velocity at the orifice entrance and can be expressed by Equation (5) [19]

$$V\_{\rm in} = -\frac{P\_{\rm s} - p}{4\eta l} [\frac{d^2}{4} - (\mathbf{x} - \mathbf{x}\_i)^2 - (z - z\_i)^2] \cdot \delta\_i \tag{5}$$

where *δ<sup>i</sup>* = 1 at the orifice entrance, whereas it is zero elsewhere.

By employing the dimensionless parameters

$$\mathbf{x} = \mathbf{x}\_0 \mathbf{X}, \mathbf{z} = \mathbf{z}\_0 \mathbf{Z}, \mathbf{h} = h\_0 \mathbf{H}, \mathbf{t} = \tau \frac{\mathbf{x}\_0}{\mathbf{u}}, \mathbf{p} = p\_a \mathbf{p}, \mathbf{A} = \frac{12 \eta \mu \mathbf{x}\_0}{h\_0^2 p\_a}, \mathbf{Q} = \frac{24 \eta \mathbf{x}\_0^2}{h\_0^3 p\_a^2} \mathbf{V}\_{\text{in}} \tag{6}$$

Equation (4) can be written as:

$$\frac{\partial}{\partial X}\left(H^3 \frac{\partial P^2}{\partial X}\right) + \frac{x\_0^2}{z\_0^2} \frac{\partial}{\partial Z}\left(H^3 \frac{\partial P^2}{\partial Z}\right) = \Lambda \frac{\partial PH}{\partial X} + 2\Lambda \frac{\partial PH}{\partial \tau} + QP \tag{7}$$

#### *2.2. Numerical Analysis*

With reference to paper [20], some items in Equation (7) are rewritten as follows by using the finite difference method:

$$\frac{\partial}{\partial X}\left(H^3 \frac{\partial P^2}{\partial X}\right) = \frac{H\_{i+0.5,j}^3}{\Delta X^2} \left(P\_{i+1,j}^2 - P\_{i,j}^2\right) - \frac{H\_{i-0.5,j}^3}{\Delta X^2} \left(P\_{i,j}^2 - P\_{i-1,j}^2\right) \tag{8}$$

$$\frac{\partial}{\partial Z} \left( H^3 \frac{\partial P^2}{\partial Z} \right) = \frac{H\_{i,j+0.5}^3}{\Delta Z^2} \left( P\_{i,j+1}^2 - P\_{i,j}^2 \right) - \frac{H\_{i,j-0.5}^3}{\Delta Z^2} \left( P\_{i,j}^2 - P\_{i,j-1}^2 \right) \tag{9}$$

$$\frac{\partial PH}{\partial X} = \frac{H\_{i+0.5,j} - H\_{i-0.5,j}}{2\Delta X} P\_{i,j} + \frac{H\_{i+0.5,j}}{2\Delta X} P\_{i+1,j} - \frac{H\_{i-0.5,j}}{2\Delta X} P\_{i-1,j} \tag{10}$$

$$\frac{\partial PH}{\partial \tau} = \frac{H\_{i,j}}{\Delta \tau} P\_{i,j} - \frac{H\_{(\tau - \Delta \tau)(i,j)}}{\Delta \tau} P\_{(\tau - \Delta \tau)(i,j)} \tag{11}$$

where *<sup>H</sup>*(*τ*−Δ*τ*)(*i*,*j*) and *<sup>P</sup>*(*τ*−Δ*τ*)(*i*,*j*) in Equation (11) are *Hi*,*<sup>j</sup>* and *Pi*,*<sup>j</sup>* in (*<sup>τ</sup>* − <sup>Δ</sup>*τ*) time, respectively.

By substituting Equations (8)–(11) into Equation (7), the equation is simplified as

$$E\_{i,j}P\_{i,j}^2 + (F\_{i,j} + I\_{i,j})P\_{i,j} - K\_{i,j} = 0\tag{12}$$

with

$$A\_{i,j} = \frac{H\_{i+0.5,j}^3}{\Delta X^2} \tag{13}$$

$$B\_{i,j} = \frac{H\_{i-0.5,j}^{3}}{\Delta X^2} \tag{14}$$

$$\mathbf{C}\_{i,j} = \frac{x\_0^2}{z\_0^2} \frac{H\_{i,j+0.5}^3}{\Delta Z^2} \tag{15}$$

$$D\_{i,j} = \frac{\varkappa\_0^2}{z\_0^2} \frac{H\_{i,j-0.5}^3}{\Delta Z^2} \tag{16}$$

$$E\_{i,j} = A\_{i,j} + B\_{i,j} + C\_{i,j} + D\_{i,j} \tag{17}$$

$$F\_{i,j} = \Lambda \frac{H\_{i+0.5,j} - H\_{i-0.5,j}}{2\Delta X} \tag{18}$$

$$G\_{i,j} = \Lambda \frac{H\_{i+0.5,j}}{2\Delta X} \tag{19}$$

$$H\_{i,j} = -\Lambda \frac{H\_{i-0.5,j}}{2\Delta X} \tag{20}$$

$$I\_{i,j} = 2\Lambda \frac{H\_{i,j}}{\Delta \tau} + Q \tag{21}$$

$$J\_{i,j} = -2\Lambda \frac{H\_{(t-\Delta \tau)(i,j)}}{\Delta \tau} \tag{22}$$

$$K\_{i\bar{j}} = A\_{i\bar{j}}P\_{i+1,\bar{j}}^2 + B\_{i\bar{j}}P\_{i-1,\bar{j}}^2 + C\_{i\bar{j}}P\_{i,\bar{j}+1}^2 + D\_{i,\bar{j}}P\_{i,\bar{j}-1}^2 - G\_{i,\bar{j}}P\_{i+1,\bar{j}} - H\_{i,\bar{j}}P\_{i-1,\bar{j}} - f\_{i,\bar{j}}P\_{(t-\Lambda\tau)(i,\bar{j})} \tag{23}$$

The dimensionless gas film pressure in *τ* time is given by

$$P\_{i,j} = \frac{\sqrt{\left(F\_{i,j} + I\_{i,j}\right)^2 + 4E\_{i,j}K\_{i,j}} - F\_{i,j} - I\_{i,j}}{2E\_{i,j}} \tag{24}$$

#### **3. Data Acquisition and Processing**

In addition to the basic parameters of the spindle, it is necessary to know the initial distribution of the gas film to solve the pressure distribution using Equation (24). That is, the *H* of any node at *τ* = 0 time is needed, which is very critical. It can be seen from Equations (2), (3), and (6) that in order to calculate *H*, the shape errors of the shaft should be measured. A cylindricity meter was used to measure the roundness and cylindricity errors for study in this paper.

#### *3.1. Roundness and Cylindricity Errors Measurement of Shaft*

In this paper, we processed 4 shafts with the same parameters and technique. The roundness and cylindricity errors of the working area of these shafts corresponding to the journal bearing was measured using the Taylor Hobson 585 LT cylindricity meter, and 21 sections of each shaft were collected. The section interval was 5 mm, and the number of sampling points was 18,000 points per section. Figure 3 shows the roundness and cylindricity error measurement process, and the cylindricity error results of shaft 1 are shown in Figure 4. The data on the left-hand side represent the runout of each section, which is equal to the roundness errors under the least square method. The data on the right side represent the measured height. The cylindricity errors of four shafts are 4.00, 4.50, 4.25, and 4.15 μm, respectively.

**Figure 3.** The roundness and cylindricity errors measurement process.

**Figure 4.** The cylindricity error measurement results of shaft 1.

#### *3.2. Data Processing*

Since the data collected by the cylindricity meter were based on relative coordinates, it was necessary to calibrate the measured data with the actual radius of the shaft. The calculation method is shown in Equation (25)

$$\begin{cases} \begin{array}{c} \mathbf{x}\_{\theta'} = \mathbf{x}\_{c\theta'} \frac{\mathbf{R}\_{\varsigma}}{\sqrt{{x\_{c\theta'}^2 + y\_{c\theta'}^2}}}\\\\ \mathbf{y}\_{\theta'} = \mathbf{y}\_{c\theta'} \frac{\mathbf{R}\_{\varsigma}}{\sqrt{{x\_{c\theta'}^2 + y\_{c\theta'}^2}}} \end{array} \end{cases} \tag{25}$$

where *xc<sup>θ</sup>* and *yc<sup>θ</sup>* are the data collected by cylindricity meter, and *Rc* is the shaft radius measured by micrometer.

In the data measured by cylindricity meter, each section had 18,000 data points. The average value of each adjacent 50 data points was taken to save the computing time, and the number of data points was reduced to 360. Gaussian filter was used for roundness filtering with 150 UPR (undulations per revolution) cut-off [21]. The number of axial nodes increased from 21 to 41 by using Equation (26)

$$h\_{i,k+0.5} = \frac{h\_{i,k} + h\_{i,k+1}}{2} \qquad (k = 1, 2, \dots, 20) \tag{26}$$

#### *3.3. Calculation Settings*

Programs were compiled according to Equations (7)–(24); Table 1 shows the parameters of the spindle. The measured data of shaft 1, shaft 2, shaft 3, and shaft 4 were used for calculation, respectively. The initial eccentricity was (*ex*,*ey*)=(0, 0), and the speed range was [500, 10,000] r/min.

**Table 1.** Spindle parameters.


#### **4. Results and Discussions**

#### *4.1. Simulation Results*

First of all, an ideal shaft without shape errors is used for the calculation to verify the correctness of the numerical model and the simulation program; Figure 5 shows the result of dimensionless pressure distribution when the shaft rotates 180◦ at 6000 r/min. It can be seen that the pressure distribution on the shaft surface is very regular, which is consistent with the calculation results in reference [22]. In this ideal case, the gas film force acting on the shaft is zero.

**Figure 5.** Dimensionless pressure distribution of ideal shaft.

Figure 6a–d show the dimensionless pressure distribution of shaft 1 at 6000 r/min when it rotates 90◦, 180◦, 270◦ and 360◦ with the same coordinate scales, respectively. Due to the roundness and cylindricity errors, the gas pressure on the shaft surface varies greatly. As the shaft rotates, the pressure changes obviously. Thus, the shaft will produce an unbalanced gas film force, which changes both in size and direction, and then affects the rotary accuracy of the spindle.

**Figure 6.** Dimensionless pressure distribution of shaft 1 at 6000 r/min at different angles: (**a**) 90◦; (**b**) 180◦; (**c**) 270◦; (**d**) 360◦.

Figure 7a–d are box charts of the surface pressure distribution of shaft 1, shaft 2, shaft 3, and shaft 4 at different speeds with the same coordinate scales. With the increase in rotating speed, the distribution of pressure value becomes more and more divergent. By comparing the interquartile range (IQR) and the range within 1.5IQR of the dimensionless pressure in Figure 7a–d, it can be observed that the pressure fluctuation of shaft 2 is significantly greater than that of the other three shafts.

**Figure 7.** Box charts of surface pressure distribution: (**a**) shaft 1; (**b**) shaft 2; (**c**) shaft 3; (**d**) shaft 4.

If the pressure is known, the gas film force acting on the shaft can be calculated using Equation (27) [23]:

$$
\begin{pmatrix} F\_x \\ F\_y \end{pmatrix} = \int\_0^L \int\_0^{2\pi} (p - p\_{atm}) \begin{pmatrix} \cos\theta \\ \sin\theta \end{pmatrix} r d\theta dz \tag{27}
$$

Figure 8a shows the gas film force (with 5 revolutions) of shaft 1 at 6000 r/min. With the rotation of the shaft, the direction and size of the force change, and the shaft adopts an unsteady state. Plotting the average force of the shaft at different speeds in Figure 8b, it can be seen that with the increase in speed, the force on the shaft increases basically. The force of shaft 2 at each speed is significantly higher than that of the other three shafts. According to the force from large to small, the order is as follows: shaft 2 > shaft 1 > shaft 4 > shaft 3.

**Figure 8.** Gas film force of shaft: (**a**) gas film force of shaft 1 (5 revolutions) at 6000 r/min; (**b**) gas film force of shafts in different speeds.

#### *4.2. Comparison with Experimental Results*

The aerostatic spindle used in this paper to validate the simulation results is mainly used for precision milling. For ensuring the stability of continuous running, the spindle is water-cooled. The four shafts are fitted into the spindle in turn, and the bearing, motor and other parts are installed. The motor is a Kollmorgen KBMS-10X01 frameless motor with the rated speed of 18,600 r/min. The driver belongs to Sieb and Meyer's SD2S series. After two stages of drying and filtering, the gas supply pressure is set to 0.4 Mpa (same as the numerical calculation). This system is mounted on a natural granite cubic base with a side length of 200 mm. For the consistency of the measurement results, other parts and related configurations remain unchanged throughout the test, except for the replacement of the shaft.

Using sensors to measure the standard gauge installed on the spindle, the rotation errors are usually measured without load, and the track of the ideal axis of the rotation of the spindle is fitted as the basis of the analysis. Due to the small size of the spindle, there is not enough space to install the standard gauge and there is a section of finish-machined cylinder in the front of the shaft as an alternative.

The dynamic radial rotation errors were measured using the Lion Precision CPL290 spindle errors analysis instrument. The bandwidth of the sensor is 15 kHz and the resolution is 0.01 μm. The spindle speed range is set to 500 to 10,000 r/min, and the data are collected every 500 r/min. Figure 9 is the test rig of the spindle rotation error measurement, and the results of shaft 1 at 1000 r/min are shown in Figure 10. Following measurement, the data of the four shafts are plotted in Figure 11.

**Figure 9.** Test rig of spindle rotation error measurement.

**Figure 10.** Result diagram of spindle radial rotation errors.

**Figure 11.** Spindle radial rotation errors of four shafts.

As can be seen from Figure 11, in the rotation speed range of [500, 10,000] r/min, the rotation radial errors of shaft 1, shaft 3 and shaft 4 are always close to each other, except for individual speeds. When the speed exceeds 2000 r/min, the errors of shaft 2 are obviously larger than those of the other three shafts, and the "particularity" of shaft 2 in Figure 11 is consistent with that in Figures 7 and 8b. In the analysis results of Figures 7 and 8, the original data are obtained from the measurement of the roundness and cylindricity errors of the shaft. It can be seen that there must be a correlation between the measurement results of the roundness and cylindricity errors of the shaft and the rotation errors.

The roundness errors of all shaft sections measured using the cylindricity meter are plotted, as shown in Figure 12. It can be seen that the roundness error lines of the four shafts overlap with each other and the difference is not obvious. For example, the roundness error values of section 11 of shaft 1 and section 15 of shaft 2 are both 2.69 μm, and the roundness error values of section 19 of shaft 2 and section 10 of shaft 4 are both 2.25 μm. Figure 13 is a composite drawing of roundness errors for section 19 of shaft 2 and section 10 of shaft 4. It can be seen that although the roundness errors values are the same (2.25 μm), their actual

shapes are quite different. Therefore, the roundness error values cannot reflect the shape of the shaft.

**Figure 12.** Roundness errors of four shafts.

**Figure 13.** Composite drawing of roundness error for section 19 of shaft 2 and section 10 of shaft 4.

#### *4.3. Evaluation with Dispersion Coefficient*

The shape of the shaft determines the gas film thickness between the shaft and the matched bearing. The more uniform the film thickness, the smaller the pressure fluctuation in the shaft; therefore, the discrete coefficient index of shaft radius can be used to quantify the section shape of the shaft. The calculation is conducted by means of Equation (28):

$$V\_S = \frac{1}{r\_0} \sqrt{\frac{1}{n-1} \sum\_{i=1}^{n} \left(r\_i - r\_0\right)^2} \tag{28}$$

where *n* is the node number of single section, *ri* is the radius of the shaft at the *i*-th node, and *r*<sup>0</sup> is the average radius of the shaft in the current section.

The radius dispersion coefficients of all sections are calculated using Equation (28) and the results are plotted in Figure 14. It can be seen that the dispersion coefficients of all sections of shaft 2 are significantly larger than those of the other three shafts, which indicates that there are great fluctuations in radius for each section of shaft 2. The particularity results here are consistent with those shown in Figures 7, 8 and 11. Although the roundness errors of section 19 of shaft 2 and section 10 of shaft 4 are identical, the differences in the dispersion coefficients shown in Figure 14 are significant. Therefore, it is not reliable to infer the shaft's rotational accuracy based on only its roundness and cylindricity errors.

To further study this problem, sections with similar roundness errors (see Table 2) and similar discrete coefficients (see Table 3) are found from all cross sections of four shafts to construct the ideal shafts (Figure 15 is a schematic diagram), and then the unbalanced gas film forces are analyzed in Figure 16. Because the shaft is constructed with a single section, the roundness errors are equal to the cylindricity errors. Figure 16a,b have the same coordinate scales to ensure the comparability of data.

**Figure 14.** Dispersion coefficients of four shafts.

**Figure 15.** Schematic diagram of an ideal shaft constructed by single section.

**Table 2.** Sections with similar roundness.


**Table 3.** Sections with similar discrete coefficients.


**Figure 16.** Unbalanced gas film force of ideal shaft: (**a**) ideal shaft with similar roundness; (**b**) ideal shaft with similar discrete coefficients.

It can be observed from Figure 16a that, although these shafts have the same roundness, the unbalanced film forces of these four shafts are quite different. The order of the force from large to small is section 19 of shaft 2 > section 7 of shaft 1 > section 10 of shaft 4 > section 1 of shaft 3, which is consistent with the order of the discrete coefficients in Table 2. As presented in Figure 16b, for the similar dispersion coefficient, the forces of the four shafts are generally similar, the lowest of which is for section 2 of shaft 1. In Table 3, the roundness errors of section 2 of shaft 1 is the largest among the four ideal shafts, but its discrete coefficient is the smallest. It can be proved that, compared with roundness errors, the discrete coefficient can be used more reliably to predict the unbalanced film force generated by the shaft due to the shape errors. In addition, when the results are combined with the data in Tables 2 and 3, it can be seen that if the roundness is similar, the discrete coefficient may be very different, but if the discrete coefficient is similar, the roundness difference is not large, so the discrete coefficient can be used to evaluate the processing quality of the shaft.

#### **5. Conclusions**

In this paper, the roundness and cylinder error data of the shaft are collected by the Taylor Hobson cylinder meter as the data source of FDM simulation. The influence of shaft shape errors on spindle rotation accuracy is verified by comparing the simulation results with experimental results. The following conclusions have been drawn.


**Author Contributions:** Conceptualization, G.Z.; methodology, J.Z.; software, G.Z. and H.Y.; validation, R.Z. and W.S.; formal analysis, J.Z. and H.Y.; investigation, J.Z. and H.Y.; data curation, G.Z. and J.W.; writing—original draft preparation, G.Z.; writing—review and editing, J.Z. and J.W.; visualization, R.Z. and W.S.; funding acquisition, J.Z. and H.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51875586; training plan for young backbone teachers in colleges and universities of Henan province, grant number 2018GGJS105.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on reasonable request from the corresponding author.

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

#### **References**


## *Article* **Porous Gas Journal Bearings: An Exact Solution Revisited and Force Coefficients for Stable Rotordynamic Performance**

**Luis San Andrés 1, Jing Yang 1,\* and Andrew Devitt <sup>2</sup>**

<sup>2</sup> New Way Air Bearings, Aston, PA 19014, USA; DDevitt@newwayairbearings.com

**\*** Correspondence: yangjing@tamu.edu

**Abstract:** Having come of age, gas film bearings enable high-speed oil-free (micro) rotating machinery with gains in efficiency and reliability, longer maintenance intervals, and a reduction in contaminants released to the atmosphere. Among gas bearing types, porous surface gas bearings (PGBs) have proven successful for 50+ years and presently are off-the-shelf mechanical elements. This paper reviews the literature on PGBs since the 1970s and reproduces an exact solution for the performance of cylindrical PGBs. Both the analytical model and an accompanying finite-element (FE) model predict the performance for two PGBs, a commercially available 76 mm diameter bearing and a smaller 25 mm diameter laboratory unit whose experimental performance is available. As expected, the FE model results reproduce the analytical predictions obtained in a minuscule computing time. For a set external supply pressure, as the radial clearance increases, the flow rate through the bearing grows until reaching a peak magnitude. The PGB load capacity is a fraction of the product of the set pressure difference (*pS* − *pa*) and the bearing projected area with a significantly large centering static stiffness evolving over a narrow region of clearances. Operation with shaft speed enhances the bearing load capacity; however, at sufficiently high speeds, significant magnitude cross-coupled forces limit the stable operation of a PGB. At constant operating shaft speed, as the whirl frequency grows, the bearing effective stiffness (*Keff*) increases, while the effective damping (*Ceff*) becomes positive for whirl frequencies greater than 50% shaft speed. Similar to a plain hydrodynamic journal bearing, the PGB is prone to a half-frequency whirl, albeit the system natural frequency can be high, mainly depending on the external supply pressure. In essence, for the cases considered, PGBs are linear mechanical elements whose load capacity is proportional to the journal eccentricity.

**Keywords:** gas bearings; porous materials; rotordynamics; force coefficients

#### **1. Introduction**

Gas bearings support oil-free high-speed turbomachinery with significant savings in drag power losses and a reduction in system footprint [1]. Bearings constructed with a layer of porous material offer an alternative to other bearing types, such as orifice compensated bearings and foil bearings [2]. Porous surface gas bearings (PGBs) are commercially available for industrial use, in particular in linear guide systems and in the health imaging industry [3].

The contemporary approach to the modeling of PGBs relies on the numerical solution of the flow equations, using to advantage desktop computational resources [4]. However, the expedience in obtaining numerical predictions often dispenses the flow physics, avoids dimensional analysis, and disregards the effects of physical parameters that carry information on limiting solutions. Hence, this paper carefully reviews the abundant published literature on PGBs and recreates elegant theoretical analyses developed in the 1970s. The predictions from the analytical model are compared against those of a computational finite element (FE) model as well as some published experimental results.

**Citation:** Andrés, L.S.; Yang, J.; Devitt, A. Porous Gas Journal Bearings: An Exact Solution Revisited and Force Coefficients for Stable Rotordynamic Performance. *Appl. Sci.* **2021**, *11*, 7949. https://doi.org/ 10.3390/app11177949

Academic Editors: Terenziano Raparelli, Andrea Trivella, Luigi Lentini and Federico Colombo

Received: 5 July 2021 Accepted: 26 August 2021 Published: 28 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> J. Mike Walker '66 Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843, USA; lsanandres@tamu.edu

#### **2. The Governing Equations for a PGB and Derived Operating Parameters**

The performance of PGBs is determined by their geometry (radial clearance *c*, length *L*, and diameter *D =* 2*R*), the gas physical properties (density *ρ* and viscosity *μ*), the magnitudes of supply pressure (*pS*) and ambient pressure (*pa*), the porous material permeability coefficient (*κ*) and the liner radial thickness (*tp*), and the operating conditions of shaft angular speed (Ω and precessional frequency (*ω*).

Figure 1 displays a schematic diagram of a cylindrical PGB with the journal spinning with speed Ω. The film thickness is *h =* (*c+eX* cos*θ + eY* sin*θ*), where (*eX*, *eY*) are the components of the journal static eccentricity. Consider an ideal gas with density *ρ* = *p*/(*RgT*), whose temperature *T* is constant. Then, the film pressure (*p*) is governed by (Gargiulo, 1979) [5]. A detailed investigation on the influence of non-zero slip flow condition and tangential flow within the porous material is out of scope in the current analysis.

$$\vec{\nabla} \cdot \left( \frac{p \ h^3}{12 \mu} \vec{\nabla} p \right) = \frac{\Omega R}{2} \frac{\partial (ph)}{\partial \mathbf{x}} + \frac{\partial (ph)}{\partial t} + \frac{\kappa}{2 \mu} \frac{}{t\_p} (p^2 - p\_S^2) \tag{1}$$

where <sup>→</sup> ∇ is the vector gradient operator, *μ* is the gas viscosity, and *pS* is the pressure of the pressurized gas supplied through the porous layer. The field *p* is periodic around the bearing circumference,*p*(*θ*, *<sup>z</sup>*,*t*) = *p*(*θ*+2*π*,*z*,*t*), and on the axial ends of the bearing *p*(*θ*, <sup>1</sup> <sup>2</sup> *<sup>L</sup>*,*t*) <sup>=</sup> *<sup>p</sup>*(*θ*,<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>L</sup>*,*t*) <sup>=</sup> *pa*.

**Figure 1.** Schematic diagram of a cylindrical journal and bearing with porous material layer (not to scale). Coordinate system noted for reference.

Note the supplied mass flow rate through the porous material follows Darcy's law [6–8].

$$
\dot{m} = \frac{\rho}{\mu} \frac{\kappa}{t\_p} (p\_S - p) \sim \frac{1}{2} \frac{1}{\mu} \frac{1}{R\_\text{gT} T} \frac{\kappa}{t\_p} \left(p\_S^{-2} - p^2\right) \tag{2}
$$

Flow unsteadiness, circumferential flow through the porous media, and inertial effects due to the volume in the porous material are not particularly important in modern commercial products such as carbon graphite [3,9,10].

The PGB geometry, porous material properties, and operating conditions combine to produce three fundamental (dimensionless) parameters,

$$
\Lambda\_{\mathbf{x}} = \frac{12\kappa}{t\_p} \left(\frac{R}{c}\right)^2; \; \Lambda\_{\Omega} = 6 \frac{\mu \,\Omega}{p\_S} \left(\frac{R}{c}\right)^2; \; \Lambda\_{\omega} = 12 \frac{\mu \,\omega}{p\_S} \left(\frac{R}{c}\right)^2. \tag{3}
$$

known as the flow feeding number, the speed number, and the frequency number, respectively. Clearly *pS* > *pa* for the pressurized gas to flow through both the porous layer and then through the film clearance to exit at the bearing sides. Most PGBs have a large slenderness ratio, *L/D* > 1, hence giving enough surface area for the pressure field to produce a significant bearing load.

The ΛΩ aerodynamic effects due to shaft rotation dominate over aerostatic effects (~*pS*). Similarly, squeeze film effects (high frequency) dominate for Λ*<sup>ω</sup>* > 1. Incidentally, the porosity parameter Λ*<sup>κ</sup>* relates the flow conductance across the porous layer (*κ/tp*) to the conductance through the film, which is proportional to (*c*3/*R*2).

#### **3. An Appraisal of the Past Literature**

The archival literature on PGBs is abundant with a sound beginning in the early 1960s and continuous developments through the 1980s. Only in the 2000s did commercial applications abound, in particular in the flat panel manufacturing industry, semiconductor manufacturing industry, coordinate measuring machines, computed tomography imaging machines, and precision machine tools and spindles [3].

From 1965 to 1967, Sneck et al. publish three seminal papers [6–8] to introduce the theoretical framework for the analysis of gas lubricated PGBs. The authors derived analytical solutions for the flow rate and load capacity of PGBs under aerostatic conditions and produced experimental verification for both the load capacity as well as the flow coefficients. The authors noted the importance of surface roughness that modifies the PGB nominal clearance to produce agreement with the measured load and flow. The later reference [8] noted that shaft surface speed aids to increase the load capacity of PGBs and presents ultimate solutions for operation under very large speed numbers (ΛΩ >> 1). Sneck et al. [6–8] did not study the effect of frequency on the bearing stiffness coefficients, and altogether ignored damping coefficients.

At about the same time, in 1968, Mori et al. [11,12] developed a similar solution, albeit for an incompressible fluid, and produced predictive formulas for load capacity and gas flow that agreed well with experimental results for finite length PGBs (*L/D ~* 1/2). Both Sneck et al. and Mori et al. found that the PGB load capacity was a fraction of *W\* =* ((*pS* − *pa*) × (*L D*)), i.e., the product of the pressure difference (*pS* − *pa*) times the bearing projected area. Having introduced an equivalent clearance for the layer of porous material, *c<sup>κ</sup>* = (12 *κ tp*) 1/3, Mori et al. reported a maximum aerostatic load capacity ~(0.7 × *W\**) for the range (*cκ*/*c*) ~ 0.6–1.0; hence, the practical flow feeding number Λ*<sup>κ</sup>* ~ (*R/tp*) <sup>2</sup> >> 1.

A decade later (1978), Rao and Majumdar [13] presented a perturbation analysis and numerical solution to calculate the stiffness (*K*) and damping (*C*) force coefficients of aerostatic PGBs (ΛΩ = 0). The unique predictions showed *C* < 0 for small frequency numbers (Λ*ω*). Interestingly enough, the analysis did not include the volume of the porous material, which, if large enough and under dynamic conditions, traps the gas to produce a pneumatic-hammer effect, i.e., a self-excited instability due to the absence (even negative) of effective damping [14].

Just a year later, in 1979, Gargiulo [5] presented a comprehensive analysis that included both steady-state performance parameters and the dynamic force coefficients via a perturbation analysis that produces analytical expressions as a function of ΛΩ, Λ*ω*, and Λ*κ*, and includes the effect of the material porosity volume. For very small porous volumes, Cargiulo [1] reported a hardening bearing stiffness as frequency grows while the damping coefficient decreases steadily. The findings are similar to those for orifice-compensated aerostatic gas bearings, as published by Lund [15] in 1968. Large porous volumes could lead to a negative direct stiffness but not negative damping, in opposition to Rao and Majumdar's findings [13]. Under hybrid operations, i.e., operation with shaft speed (ΛΩ > 0), the analysis reports cross coupled stiffness and damping coefficients, growing in proportion to the shaft speed until ΛΩ → 10, to then drop for operation at higher shaft speeds. Cargiulo [5] left for the future (the 21st century) the dynamic stability of PGBs operating in a hybrid mode.

In a companion paper, Cargiulo [16] described the outcome of experiments that for the most part validated the theoretical predictions, in particular the PGB load capacity and the direct force coefficients. The test bearing has *L=D=* 50.4 mm and a porous liner with thickness *tp* = 6.4 mm and permeability *<sup>κ</sup>* = 2.3 × <sup>10</sup>−<sup>9</sup> mm2 (*c<sup>κ</sup>* = 5.6 <sup>μ</sup>m). The experiments also revealed the onset of unstable behavior characterized by vibrations of the test system at its natural frequency. The tests conducted with various radial clearances, *c =* 3 to 17 μm, demonstrated that the test aerostatic PGB produced a journal displacement directly proportional to the applied load, even up to the point of journal contact with the porous liner surface. Hybrid operation, i.e., with shaft rotation (ΛΩ > 0), increases the PGB load capacity and its stiffnesses, direct and cross-coupled. Incidentally, the higher the supply pressure, the larger the bearing centering stiffness (and natural frequency), and hence the more likelihood of a PGB to be tailored for dynamically stable operation.

After a long hiatus, in 2006, Miyatake et al. [17] assessed the stability of a rotor supported on PGBs coated with a surface layer that was ~300 times more restrictive than that of the porous media itself. The manufacturing of surface loaded porous materials is a major breakthrough in PGBs, as they practically eliminate the potential for pneumatic hammer; see Otsu et al. [18].

In practice, PGBs are mounted in series with elastomeric supports (o-rings) to add more effective damping, albeit the soft mounts also reduce the system natural frequency. See for example the synopsis by Hwang [14]. Most recently, various publications detail the application of tilting pad PGBs to enable oil-free turbomachinery operating at high speeds; see San Andrés et al. [2,4,19,20] and Feng et al. [21–23]. Tilting pad journal bearings effectively eliminate the cross-coupled stiffness coefficient that could easily excite a hydrodynamic instability and thus extend the maximum operating shaft speed of rotor-bearing systems.

At Texas A&M University, since 2015, San Andrés et al. have been conducting measurements of the rotordynamic response of solid rotors supported on carbon-graphite PGBs: one rotor of small diameter (29 mm) and operating at a high shaft angular speed (55 krpm) [19], and another with large diameter (100 mm) and turning at 18 krpm maximum [4,20]. The rotors' surface speeds (*Us =* Ω*R*) equal 82 m/s and 94 m/s, respectively. Rotordynamic tests show stable responses, free of sub synchronous whirl frequencies. Imbalance mass induced synchronous whirl speed rotor motion amplitudes show peak magnitudes while crossing critical speeds largely determined by the pressure (*pS*) supplied to the bearings. Derived from recorded amplification factors, the system damping ratio is larger than 10%, uncharacteristically high for gas bearing supported rotors. The computational analysis in [4] solves the Reynolds equation and Darcy's diffusion equation that couple the flow in the PGB film land to that through the porous layer. The model produces force coefficients that agree well with experimentally-derived effective damping and stiffness coefficients and shows that the magnitude of supply pressure (*pS*), the permeability coefficient (*κ*), the pad-pivot stiffness, and the assembly clearance affect significantly the performance, static and dynamic, of the test bearings. In addition, just published in 2021, [2] describes an experimental campaign and the results that quantify the force coefficients of a large size, four tilting pad PGB. From dynamic load experiments conducted with shaft speeds at 6 and 9 krpm (32 and 64 m/s surface speed), the identified PGB's direct stiffness and damping coefficients are practically invariant with excitation frequency (to max 200 Hz), The pads' pivot compliance largely determines the bearing direct stiffness coefficients. The experiments also reveal that the load capacity of the bearing is mainly aerostatic, mostly determined by the supply pressure and the bearing preload assembly.

Feng at Hunan University leads an extensive research program whose objective is to develop oil-free bearing technologies applicable to high speed rotating machinery. Feng's group has effectively tackled bump-type foil GBs, metal mesh foil GBs, and most recently, PGBs since 2018. Feng et al. develop computational analyses for tilting pad PGBs [21,22] and cylindrical PGBs [23] and present model validations against their own experimental results as well as those in [20].

Other minor recent computational developments by Bohle [24] and Lawrence and Kemple [25] provide numerical steps toward modeling aerostatic PGBs while ignoring the vast prior literature reviewed above. Li et al.'s recent work on cylindrical PGBs [23] is particularly of notice, since the recorded static load performance of the test PGB agrees with numerical predictions. However, the identified test bearing force coefficients are rather unique since they follow a trend opposite to prior knowledge, both experimental and analytical [5,16]. That is, the PGB direct stiffness coefficient drops with frequency, while the direct damping coefficient increases.

The recent literature from 2018 through 2021 [2,21,24,25] focuses on building complex comprehensive computational tools with model verifications specific to a tailored configuration for oil-free machinery. The literature misses careful dimensional analysis, physical scaling, and dimensionless number representation that could help improve the selection of the numerical method and its validation.

This paper reworks the classical analyses of Gargiulo [5], Sneck [6], and Mori [11] to derive an exact solution for the ultimate flow rate and flow coefficient, static load capacity, and dynamic force coefficients of a generic cylindrical PGB. In the last stage of the analysis, the evaluation of integrals, the current development relies on modern mathematical software capable of symbolic algebra processing. The authors wish the engineering student and neophyte researcher appreciate the art of fundamental physical analysis, while staying for a short time away from the mundane routine of number crunching.

#### **4. A Close Form Solution to the Flow and Dynamic Force Coefficients in a PGB**

In reference to Figure 1, consider journal center motions with small amplitude (*δeX*, *δeY*) and frequency (*ω*) about an equilibrium position (*eX*0, *eY*0). The film thickness (*h*) and pressure (*p*) fields equal the sum of zeroth-order fields (*h*0, *p*0) and first-order or perturbed fields (*h*σ, *p*σ)|<sup>σ</sup> <sup>=</sup> *<sup>X</sup>*, *<sup>Y</sup>*,

$$h = h\_0 + \delta c\_X \ln\_X e^{i\omega t} + \delta c\_Y \ln\_Y e^{i\omega t} \to p = p\_0 + \delta c\_\sigma \ p\_\sigma \ e^{i\omega t}, \sigma = \text{X}, \text{Y} \tag{4}$$

where *hX <sup>=</sup>* cos(*θ*) and *hY <sup>=</sup>* sin(*θ*). Above, <sup>σ</sup> <sup>=</sup> *<sup>X</sup>*, *<sup>Y</sup>*, *<sup>i</sup>* <sup>=</sup> √−1.

Substituting Equation (4) into the Reynolds Equation (1) produces the zero-order equation for *p*<sup>0</sup>

$$\overrightarrow{\nabla} \cdot \left(\frac{p\_0 h\_0^3}{12\mu} \overrightarrow{\nabla} p\_0\right) = \frac{\Omega R}{2} \frac{\partial (p\_0 h\_0)}{\partial \mathbf{x}} + \frac{\kappa}{2\mu} \left(p\_0^2 - p\_S^2\right) \tag{5}$$

and the first-order equation

$$\begin{aligned} \stackrel{\rightarrow}{\nabla} \cdot \Big( \begin{array}{c} \frac{\mu\_0}{12\mu} \Big[ \begin{array}{c} \stackrel{\rightarrow}{\nabla} p\_\sigma \stackrel{\rightarrow}{\nabla} p\_0 + p\_0 \stackrel{\rightarrow}{\nabla} p\_\sigma \end{array} \Big] \Big) \\ \stackrel{\rightarrow}{\partial} \Big( \begin{array}{c} \frac{\mu\_0}{4\mu} \stackrel{\rightarrow}{\partial} \mathbf{x} \end{array} \Big) \Big. \\ \stackrel{\rightarrow}{\nabla} \cdot \Big( \begin{array}{c} \frac{p\_0 \mu\_0}{4\mu} \stackrel{\rightarrow}{\nabla} p\_\sigma \stackrel{\rightarrow}{\nabla} p\_0 \end{array} \Big) + \frac{\kappa}{\mu \frac{\Gamma}{\Gamma\_P}} p\_0 \, p\_\sigma \end{aligned} \Big( \begin{array}{c} \sigma \end{array} = X \text{, } Y \end{aligned} \Big) \begin{aligned} \stackrel{\rightarrow}{\nabla} \begin{aligned} \begin{array}{c} \frac{\Gamma}{\Gamma\_P} \stackrel{\rightarrow}{\nabla} p\_\sigma \end{aligned} \Big) \Big. \end{aligned} \Big) \begin{aligned} \stackrel{\rightarrow}{\nabla} \begin{aligned} \frac{\Gamma}{\Gamma\_P} \stackrel{\rightarrow}{\nabla} p\_\sigma \end{aligned} \Big) \end{aligned}$$

Following San Andrés [26], the equations above are integrated in a typical finite element (FE) and assembled over the flow domain to produce sets of algebraic (nonlinear) equations for solution of the pressure fields; and subsequent calculation of the PGB reaction force, drag torque and power, and dynamic force coefficients. Details of the FE model are omitted for brevity.

For operation at a centered condition, *h*<sup>0</sup> = *c*, Equation (5) reduces to

$$\stackrel{\rightarrow}{\nabla} \cdot \left( \stackrel{\rightarrow}{\nabla} p\_0 \,^2 \right) = \frac{12 \,\text{\kappa}}{c^3 \, t\_p} \left( p\_0 \,^2 - p\_S^2 \right) \tag{7}$$

That dispenses with the effects of gas viscosity (*μ*) and shaft speed (Ω). Let <sup>ψ</sup><sup>0</sup> <sup>=</sup> *p*2 <sup>0</sup> − *<sup>p</sup>*<sup>2</sup> *S* , then Equation (7) reduces to the simple ordinary differential equation,

$$\frac{d^2}{d\mathbb{Z}^2}(\psi\_0) - \Lambda\_k \,\psi\_0 = 0. \tag{8}$$

Above, *z* = *<sup>z</sup> R* with *γ*<sup>2</sup> = Λ*<sup>κ</sup>* = <sup>3</sup> *<sup>κ</sup> tp c <sup>D</sup> c* 2 as the feed flow parameter. The solution of Equation (8) with ambient pressure (*pa*) at the bearing sides *z* = ± *<sup>L</sup> D* gives

$$\psi\_{0(\mathsf{T})} = \psi\_0|\_{\mathsf{P} = p\_a} \frac{\cosh(\gamma \overline{\mathsf{z}})}{\cosh\left(\gamma \frac{\mathsf{L}}{\mathsf{D}}\right)} \rightarrow p\_{0(\mathsf{T})} = \sqrt{p\_S^2 + \left(p\_a^2 - p\_S^2\right) \frac{\cosh\left(\gamma \overline{\mathsf{z}}\right)}{\cosh\left(\gamma \frac{\mathsf{L}}{\mathsf{D}}\right)}}.\tag{9}$$

The mass flow rate ( . *M*) supplied to the bearing through the porous layer equals

$$\dot{M} = \iiint \dot{m} \, R \, d\theta \, dz = \gamma \left( \frac{\pi}{6} \frac{c^3}{\mu} \right) \tanh \left( \gamma \frac{L}{D} \right) \, \frac{\left( p\_S^{-2} - p\_d^{-2} \right)}{R\_\% T}. \tag{10}$$

Note that . *<sup>M</sup>* is proportional to #*<sup>κ</sup> <sup>c</sup>*<sup>3</sup> *tp* , the average gas density *<sup>ρ</sup>*<sup>∗</sup> <sup>=</sup> <sup>1</sup> 2 (*pS*+*pa* ) *RgT* , and the pressure difference (*pS* − *pa*). If the bearing clearance (*c*) is very large (as when there is no shaft inserted in the bearing), then *p=pa* over the flow domain, and the supplied flow rate is the largest and equal to

$$\dot{M}\_{\text{max}} = \frac{1}{2} \frac{\kappa}{\mu} \frac{\left(p\_S^2 - p\_a^2\right)}{R\_\text{\textdegree C}} (\pi \,\, D \,\, L) \tag{11}$$

$$\frac{\dot{M}}{\dot{M}\_{\text{max}}} = \frac{\tanh\left(\gamma \frac{L}{D}\right)}{\left(\gamma \frac{L}{D}\right)}.\tag{12}$$

The ratio of flows is only a function of the bearing geometry and the permeability coefficient (*κ*), since *γ <sup>L</sup> D* = # 3*κ tpc L <sup>c</sup>* . Equation (11), known as the flow coefficient for PGBs [6,7], serves to estimate the permeability coefficient (*κ*) from measurements of the supplied flow ( . *M*max) for various supply pressures, as will be shown later.

Incidentally, the shear drag torque (*To*) and power loss (*Ploss*) for the laminar flow in a centered journal bearing equal

$$T\_o = \mu \,\,\Omega \,\, \text{(R/c)} \,\, (\pi D \,\, L), \\
 P\_{loss} = T\_o \,\, \Omega. \tag{13}$$

Note that for an air film, the drag friction coefficient *f=To/(W\* R)* << 1.

As the analysis considers small amplitude motions about the centered condition (*e =* 0), then from *p* = *p*<sup>0</sup> + (*δeX pX* + *δeY pY*)*e<sup>i</sup> <sup>ω</sup> <sup>t</sup>* , deduce

$$p^2 = p\_0^2 + 2\ p\_0 \left(\delta \varepsilon\_X p\_X + \delta \varepsilon\_Y p\_Y\right) \; \varepsilon^{i\ \omega \ t} ; \; \psi = \left(p^2 - p\_S^2\right) = \psi\_0 + \left(\delta \varepsilon\_X \psi\_X + \delta \varepsilon\_Y \psi\_Y\right) \; \varepsilon^{i\ \omega \ t} . \tag{14}$$

Hence, ψ*<sup>X</sup>* = 2*p*<sup>0</sup> *pX*; ψ*<sup>Y</sup>* = 2*p*<sup>0</sup> *pY*. Presently, for small amplitude motions about *e =* 0, let

$$
\psi\_X = \psi\_{X\varepsilon(z)}\cos\theta \,\, +\psi\_{X\varepsilon(z)}\sin\theta \,\, \psi\_Y = \psi\_{Y\varepsilon(z)}\cos\theta \,\, +\psi\_{Y\varepsilon(z)}\sin\theta \,. \tag{15}
$$

Then, write the first order Equation (6) for the perturbed pressure fields as

$$\begin{aligned} \frac{d^2}{d\tau^2} \begin{bmatrix} \boldsymbol{\Psi}\_{X\boldsymbol{\varepsilon}} \\ \boldsymbol{\Psi}\_{X\boldsymbol{\varepsilon}} \end{bmatrix} - \mathbf{A} \begin{bmatrix} \boldsymbol{\Psi}\_{X\boldsymbol{\varepsilon}} \\ \boldsymbol{\Psi}\_{X\boldsymbol{\varepsilon}} \end{bmatrix} &= \mathbf{C}\_1 \cosh(\gamma \overline{\boldsymbol{\varepsilon}}) + \mathbf{C}\_2; \\\ \frac{d^2}{d\tau^2} \begin{bmatrix} \boldsymbol{\Psi}\_{Y\boldsymbol{\varepsilon}} \\ -\boldsymbol{\Psi}\_{Y\boldsymbol{\varepsilon}} \end{bmatrix} - \mathbf{A}^T \begin{bmatrix} \boldsymbol{\Psi}\_{Y\boldsymbol{\varepsilon}} \\ -\boldsymbol{\Psi}\_{Y\boldsymbol{\varepsilon}} \end{bmatrix} &= \mathbf{C}\_1 \cosh(\gamma \overline{\boldsymbol{\varepsilon}}) + \mathbf{C}\_2 \end{aligned} \tag{16}$$

$$\begin{aligned} \text{where } \mathbf{A} &= \begin{bmatrix} \begin{pmatrix} 1 + \gamma^2 \\ -\overline{\Lambda}\_{\Omega} \end{pmatrix} + i\overline{\Lambda}\_{\omega} & \begin{matrix} \overline{\Lambda}\_{\Omega} \\ 1 + \gamma^2 \end{matrix} + i\overline{\Lambda}\_{\omega} \\ &= \begin{bmatrix} i 2\overline{\Lambda}\_{\omega} - 3\gamma^2 \\ -\overline{\Lambda}\_{\omega} \end{bmatrix} \underline{\underline{\Psi}\_{\text{eff}}} \xrightarrow{\overline{\underline{1}}} \text{and } \mathbf{C}\_{\omega} = \begin{bmatrix} i 2\overline{\Lambda}\_{\omega} & \frac{\underline{\nu}^2}{2} \end{bmatrix} \end{aligned} \tag{17}$$

$$\mathbf{C}\_{1} = \begin{bmatrix} \,^1\Sigma\Lambda\_{\omega} - \,^3\gamma \,^4 \\ -2\overline{\Lambda}\_{\Omega} \end{bmatrix} \begin{array}{c} \frac{\overline{\Psi\_{\alpha}}}{c} \frac{1}{\cosh\left(\gamma\frac{\pi}{6}\right)} \text{ and } \mathbf{C}\_{2} = \begin{bmatrix} \,^1\Sigma\Lambda\_{\omega} \\ -2\overline{\Lambda}\_{\Omega} \end{bmatrix} \begin{array}{c} \frac{p\_{S}^{\omega}}{c} \\ \frac{\overline{\Psi\_{\alpha}}}{c} \end{array}$$

$$\text{with } \psi\_{as} = \left(p\_a^2 - p\_S^2\right), \overline{\Lambda}\_{\Omega} = \frac{12\mu}{2\overline{p}\_0} \left(\frac{R^2}{c^2}\right); \ \overline{\Lambda}\_{\omega} = \frac{24\mu}{2\overline{p}\_0} \left(\frac{R^2}{c^2}\right) \tag{18}$$

using average pressure *<sup>p</sup>*<sup>0</sup> = <sup>1</sup> *L* ´ *L* <sup>0</sup> *p*<sup>0</sup> *dz*. The boundary conditions for the perturbed fields are homogeneous, i.e.,ψ*X*(*z*=<sup>±</sup> *<sup>L</sup> <sup>D</sup>* ) <sup>=</sup> <sup>ψ</sup>*Y*(*z*=<sup>±</sup> *<sup>L</sup> <sup>D</sup>* ) <sup>=</sup> 0. Then,

$$
\psi\_{X\varepsilon} = \psi\_{Y\varepsilon} = \,\,\psi\_{X\mathfrak{s}} = \psi\_{Y\mathfrak{c}} = 0 \,\,\text{at}\,\,\overline{z} = \pm \frac{L}{D} \,. \tag{19}
$$

Note that Equation (16) reveals that ψ*Xc* = ψ*Ys*, ψ*Xs* = −ψ*Yc*.

Once found, integration of the pressure field (*p*) acting on the rotor surface gives the bearing reaction force (**F**) with components

$$-\left[\begin{array}{c} F\_X\\F\_Y \end{array}\right] = \oint \left(P\_0 - P\_a + \left(\delta\mathfrak{c}\_X \, P\_X + \,\delta\mathfrak{c}\_Y \, P\_Y\right) \, e^{i\,\omega \cdot t}\right) \begin{bmatrix} \cos\theta\\\sin\theta \end{bmatrix} \, Rd\theta dz. \tag{20}$$

As the bearing is centered, the zeroth order pressure field does not produce a static force, i.e., **F** = 0. The first order pressure fields produce the matrix of complex dynamic stiffnesses **H**, where *Hij* = *Kij* + *i ω Cij* i,j = X,.

$$\begin{bmatrix} H\_{XX} \\ \end{bmatrix} = \begin{bmatrix} H\_{YY} \\ -H\_{XY} \end{bmatrix} = \oint -p\_X \begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix} \mathbb{R} \, d\theta dz = -\oint \begin{pmatrix} \Psi\_X \\ 2 \ p\_0 \end{pmatrix} \begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix} \mathbb{R} \, d\theta dz. \tag{21}$$

Then substituting Equation (15) above leads to

$$H\_{XX} = H\_{YY} = -\pi R^2 \int\_0^{L/D} \frac{\Psi\_{Xc\left(\Xi\right)}}{p\_{0\left(\Xi\right)}} \, d\Xi;\,\, H\_{XY} = -H\_{YX} = -\pi R^2 \int\_0^{L/D} \frac{\Psi\_{Yc\left(\Xi\right)}}{p\_{0\left(\Xi\right)}} \, d\Xi.\tag{22}$$

The exact solution for the first-order fields is as follows: Let <sup>ψ</sup>*<sup>X</sup>* <sup>=</sup> <sup>ψ</sup>*Xc* <sup>ψ</sup>*Xs <sup>T</sup>* , then Equation (16) has the general solution

$$
\Psi\_X = \psi\_{XP1} \cosh(\gamma \overline{z}) + \psi\_{XP2} + \,\,\Psi\_{XH} \tag{23}
$$

where the coefficients of the particular solution are

$$\boldsymbol{\Psi}\_{X\boldsymbol{\Psi}1} = -\left[\mathbf{A} - \mathbf{I}\boldsymbol{\gamma}^2\right]^{-1}\mathbf{C}\_1 = -\frac{1}{\Delta} \begin{bmatrix} 1 + i\overline{\Lambda}\_{\omega} & -\overline{\Lambda}\_{\Omega} \\ \overline{\Lambda}\_{\Omega} & 1 + i\overline{\Lambda}\_{\omega} \end{bmatrix} \begin{bmatrix} i\,2\overline{\Lambda}\_{\omega} - 3\,\boldsymbol{\gamma}^2 \\ -2\overline{\Lambda}\_{\Omega} \end{bmatrix} \frac{1}{\cosh\left(\gamma\frac{\boldsymbol{I}}{D}\right)} \frac{\psi\_{as}}{c} \tag{24}$$

$$\Psi\_{X}p\_{2} = -\mathbf{A}^{-1}\mathbf{C}\_{2} = -\frac{1}{\Delta\_{\gamma}} \begin{bmatrix} \left(1+\gamma^{2}\right) + i\overline{\Lambda}\_{\omega} & -\overline{\Lambda}\_{\Omega} \\ \Lambda\_{\Omega} & \left(1+\gamma^{2}\right) + i\overline{\Lambda}\_{\omega} \end{bmatrix} \begin{bmatrix} i\,2\overline{\Lambda}\_{\omega} \\ -2\overline{\Lambda}\_{\Omega} \end{bmatrix} \frac{p\_{S}^{2}}{c} \tag{25}$$

with Δ = 1 + *i*Λ*<sup>ω</sup>* <sup>2</sup> <sup>+</sup> ΛΩ<sup>2</sup> and <sup>Δ</sup>*<sup>γ</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> *<sup>γ</sup>*<sup>2</sup> + *i*Λ*<sup>ω</sup>* <sup>2</sup> <sup>+</sup> ΛΩ<sup>2</sup> as determinants of the system of equations. <sup>ψ</sup>*XH* <sup>=</sup> - ψ*Hesz*, the homogenous solution of Equation (16), satisfies

$$\left[\mathbf{I}\,\mathbf{s}^2 - \mathbf{A}\right] \widehat{\boldsymbol{\Psi}}\_H = \mathbf{0} \tag{26}$$

whose roots are *x*1,2 = (*s*1,2) <sup>2</sup> = 1 + *γ*<sup>2</sup> + *i* Λ*<sup>ω</sup>* ± ΛΩ , and with corresponding eigenvectors - <sup>ψ</sup>*H*<sup>1</sup> <sup>=</sup> <sup>1</sup> *<sup>i</sup> T* , - <sup>ψ</sup>*H*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup>*<sup>i</sup> T* . Thus,

$$
\Psi\_{XH} = \stackrel{\frown}{\Psi\_H} \epsilon^{\mathfrak{F}} = d\_{1c} \stackrel{\frown}{\Psi\_{H\_{1c}}} \cosh(\sqrt{\mathfrak{x}\_1}\overline{\mathfrak{z}}) + d\_{2c} \stackrel{\frown}{\Psi\_{H\_{2c}}} \cosh(\sqrt{\mathfrak{x}\_2}\overline{\mathfrak{z}}).\tag{27}
$$

Satisfying the end condition ψ*X*(*z*<sup>=</sup> *<sup>L</sup> <sup>D</sup>* ) <sup>=</sup> **<sup>0</sup>** leads to

$$
\begin{bmatrix} d\_{1\varepsilon} \\ d\_{2\varepsilon} \end{bmatrix} = - \begin{bmatrix} \cosh\left(\sqrt{\mathbf{x}\_{1}} \frac{l}{D}\right) & \cosh\left(\sqrt{\mathbf{x}\_{2}} \frac{l}{D}\right) \\\ i \cosh\left(\sqrt{\mathbf{x}\_{1}} \frac{L}{D}\right) & -i \cosh\left(\sqrt{\mathbf{x}\_{2}} \frac{L}{D}\right) \end{bmatrix}^{-1} \begin{bmatrix} \Psi\_{XP\_{1}} \cosh\left(\gamma \frac{L}{D}\right) + \Psi\_{XP2} \end{bmatrix} . \tag{28}
$$

Finally, the exact solution is

$$\Psi\_X = \begin{bmatrix} \ \Psi\_{X\varepsilon} & \Psi\_{X\varepsilon} \end{bmatrix}^T = \Psi\_{X\mathbb{P}\_1} \cosh(\gamma \overline{z}) + \Psi\_{X\mathbb{P}\_2} + d\_{1\varepsilon} \begin{bmatrix} 1 \\ i \end{bmatrix} \cosh(\sqrt{x\_1}\overline{z}) + d\_{2\varepsilon} \begin{bmatrix} 1 \\ -i \end{bmatrix} \cosh(\sqrt{x\_2}\overline{z}). \tag{29}$$

Equation (22) implements the functions above, and using symbolic mathematical software, produces the force coefficient in a very short time (fractions of a second). In 1979, Gargiulo [5] produced a nearly identical analysis, albeit he used an influence coefficient method to obtain the solution of the ordinary differential equations (the authors were unware of the original solution in [5] until later in their work when searching for examples to validate their FE and analytical models. There is no excuse for the ignorance of past work).

#### **5. PGB Aerostatic Operation**

It is of interest to quantify the bearing operation under aerostatic conditions (without shaft speed), in particular to predict the PGB static stiffness (*KS*). Let Ω = *ω* = 0, then ΛΩ = Λ*<sup>ω</sup>* = 0, and Equation (16) reduces to

$$\frac{d^2}{d\,\overline{z}^2} \begin{bmatrix} \Psi\_{Xc} \\ \Psi\_{Xs} \end{bmatrix} - \left(1 + \gamma^2\right) \begin{bmatrix} \Psi\_{Xc} \\ \Psi\_{Xs} \end{bmatrix} = \begin{bmatrix} -3\,\gamma^2 \\ 0 \end{bmatrix} \frac{\Psi\_{as}}{c} \frac{\cosh(\gamma\overline{z})}{\cosh\left(\gamma\frac{L}{D}\right)}\tag{30}$$

whose solution is

$$\Psi\_{\rm Xc(\underline{r})} = \frac{3}{\varepsilon} \gamma^2 \left( p\_a^2 - p\_S^2 \right) \left\{ \frac{\cosh(\gamma \overline{z})}{\cosh\left(\gamma \frac{L}{D}\right)} - \frac{\cosh\left[ \left( 1 + \gamma^2 \right)^{0.5} \overline{z} \right]}{\cosh\left[ \left( 1 + \gamma^2 \right)^{0.5} \frac{L}{D} \right]} \right\}, \quad \Psi\_{\rm Xs} = 0. \tag{31}$$

Since ψ*<sup>X</sup>* = 2*p*<sup>0</sup> *pX*; ψ*<sup>Y</sup>* = 2*p*<sup>0</sup> *pY*, then

$$p\_X = \frac{1}{2p\_0} \left[ \psi\_{X\varepsilon(z)} \cos \theta \right], \ p\_Y = \frac{1}{2p\_0} \left[ \psi\_{Y\varepsilon(z)} \sin \theta \right]. \tag{32}$$

Substitute into Equation (21) to obtain the static stiffness*KS* = *KXXS* = *KYYS* , and *KXYS* = *KYXS* = 0.

#### **6. Experimental Estimation of Porous Material Permeability Coefficient (***κ***)**

Figure 2 displays a photograph of a cylindrical PGB with physical dimensions listed in Table 1. The table also details the conditions for measurement of the supplied flow rate and estimation of *κ*. A product specification sheet is available in [27].

**Figure 2.** Photograph of a porous material cylindrical journal bearing.


**Table 1.** Dimensions of porous journal bearing and operating conditions.

Without a shaft installed in the bearing, measurements of the air supply pressure (*pS*) and ensuing mass flow rate ( . *M*max) deliver an estimation of the permeability coefficient (*κ*) [20]. From Equation (11), .

$$\kappa = \frac{M\_{\text{max}}}{\pi \, D \, L} \frac{2 \, t\_p \, \mu \, R\_{\text{g}} \, T}{(p\_{\text{S}}^2 - p\_{\text{a}}^2)}. \tag{33}$$

In the measurements, the supply pressure (*pS*) increases from 2.4 bar to 6.5 bar (absolute), and an accurate turbine type meter records the flow. Figure <sup>3</sup> shows . *M*max vs. (*pS* <sup>2</sup> − *pa* 2), and the line fitting the data has a correlation factor *R*<sup>2</sup> = 1.0. From the experimental data, *<sup>κ</sup>* = 8.2 × <sup>10</sup>−<sup>16</sup> m2, which is typical for the product.

**Figure 3.** PGB measured mass flow rate . *M*max vs. (*pS* <sup>2</sup> <sup>−</sup> *pa* 2) and a line fit (no shaft installed).

#### **7. PGB Mass Flow Rate, Peak Pressure, and Aerostatic Stiffness vs. Clearance**

From Equation (9), the pressure at the bearing mid-plane (*z =* 0) is *p*(*z*=0) = \* *p*2 *<sup>S</sup>* <sup>+</sup> *p*2 *<sup>a</sup>* − *<sup>p</sup>*<sup>2</sup> *S* ( cosh *γ <sup>L</sup> D* )−<sup>1</sup> . For the PGB with geometry listed in Table 1 and *<sup>e</sup>* = 0, Figure <sup>4</sup> depicts the mass flow rate ( . *M*),*p*(*z*=0)and the aerostatic stiffness *KS* (= *KXXs* = *KYYs*) vs. clearance for the PGB supplied with air at *pS =* 6 bar. Note the equivalent clearance for the porous layer *c<sup>κ</sup>* = (12 *κ tp*) 1/3~3 μm.

The graphs include results from the exact solution and the FE computational model (1872 elements = 72 circumferential × 26 axial). Note that shaft speed has no effect on both parameters; see Reynolds Equation (7). The FE model flowrate is almost identical to the exact solution (difference less than 1%). For small clearances (*<sup>c</sup>* < 0.040 mm), . *M* decreases quickly, *p*(*z*=0) → *pS*, *KS* reaches a maximum (184 MN/m) at *c =* 0.010 mm, and then sharply decreases as *<sup>c</sup>* <sup>→</sup> 0. On the other hand, for *<sup>c</sup>* > 0.06 mm, . *<sup>M</sup>* <sup>→</sup> . *<sup>M</sup>*max *<sup>c</sup>*→<sup>∞</sup> = 0.73 g/s, *p*(*z*=0) → *pa*, and *KS* → 0, i.e., a complete loss of load carrying capacity. Note *c* = 0.010 mm (*D/c* = 7610) ~ 3.3 *cκ*, though appearing to be tight (too small bearing clearances increase cost and make installation difficult), is quite appropriate for a gas bearing. Hence, the selection of the PGB clearance is most important to both keep a low flow rate while also providing enough load capacity and centering ability.

Having established the validity of the analytical solution versus a numerical solution, Figure 5 depicts *p*(*<sup>z</sup>* = 0)/*pa*, mass flow rate, and the (centering) aerostatic stiffness (*KS*) vs. clearance (*c*) for operation with supply pressure (*pS*) = 2, 4, 6, and 8 bar. The clearance ranges from a thin *c* = 0.002 mm (*D/c* = 38,105) to *c* = 0.1 mm (*D/c* = 762). Clearly, the film peak pressure, mass flow, and aerostatic stiffness increase with *pS.* As noted earlier, too large clearances (*c* > 0.040 mm → *D/c* < 1,900) produce no film pressure (*p* → *pa*); hence, the PGB leaks too much (<sup>→</sup> . *<sup>M</sup>*max *<sup>c</sup>* <sup>→</sup><sup>∞</sup> ) and is devoid of load capacity (*KS* <sup>→</sup> 0). For a very tight

clearance (*<sup>c</sup>* < 0.005 mm), the mid-plane pressure *pz=*<sup>0</sup> <sup>→</sup> *pS*, the flow ( . *M*) is small, and the direct stiffness is along a decreasing path. Independent of the supply pressure, the peak aerostatic stiffness appears at the same (narrow) clearance, *c ~* 0.010 mm = 3.3 *cκ*. At this *c*, *pz =* <sup>0</sup>*/pa* = 0.92 and . *<sup>M</sup>*/ . *M*max = 0.37.

**Figure 4.** PGB mass flow, pressure at bearing mid-plane, and aerostatic stiffness vs. clearance (*c*); (**a**) Mass flow rate; (**b**) Pressure at mid-plane *p*(*<sup>z</sup>* = 0)/*pa*; (**c**) Static stiffness *KS*. Analytical and FE predictions. Air supply pressure *pS* = 6 bar. Aerostatic operation at *e =* 0.

**Figure 5.** PGB static pressure at mid-plane, mass flow rate, and aerostatic stiffness vs. clearance (*c*). Air supply pressure *ps* = 2, 4, 6, and 8 bar; (**a**) Mass flow rate; (**b**) Pressure at mid-plane *p*(*<sup>z</sup>* = 0)/*pa*; (**c**) Static stiffness *KS*. Aerostatic operation (Ω = 0) and null eccentricity (*e* = 0). Vertical lines denote clearance *c =* 0.01 mm for largest stiffness.

Although not shown for brevity, a dimensionless aerostatic stiffness scaled as *K =* (*KS c*/*W*∗) clusters the data in Figure 5c. At *c* = 0.010 mm, *K ~* 0.6 for all *pS.* Being of order (1), the parametrization of *K* allows one to quickly estimate the (maximum) centering stiffness of any PGB having the same (*κ*/*tp*). The finding is nearly the same as that reported by Mori et al. [11,12] in 1968, *K* ~ 0.7.

The PGB product specification sheet [27] recommends a PGB with clearance *c =* 0.010 mm and, for operation with pressure supply *pS =* 5.15 bar [60 psig], quotes a static stiffness *KS =* 159 MN/m and flow of 13.2 L per minute (LPM) at standard conditions. For the same conditions, the current analysis predicts *KS =* 155.5 MN/m and a flow rate equaling to 9.86 LPM. The agreement in stiffness is remarkable, while the difference in flows indicates differing clearances. For example, for *c =* 0.012 mm (~20% larger), the analysis predicts 13.2 LPM flow and *KS =* 149 MN/m. Since the porous layer is affixed to the bearing housing using an adhesive, the manufacturer [27] specifies a maximum operating supply pressure of 7.9 bar (100 psig).

#### **8. PGB Force Coefficients vs. Rotor Speed (Synchronous Frequency Condition)**

The analysis continues for the PGB with dimensions in Table 1 and for clearance *c\* =* 0.010 mm, the one providing the peak aerostatic stiffness (*KS*). The operating speed range covers 0 to 30 krpm, with a mean operating speed (MOS) Ω\* = 25 krpm (417 Hz); hence the shaft surface speed Ω*R* = 120 m/s. The pressure supply *pS* = 2 bar → 8 bar, and the ambient pressure *pa* = 1 bar. Note the feed flow parameter Λ*<sup>κ</sup>* = *γ*<sup>2</sup> = 5.3; and at the lowest supply pressure and MOS, ΛΩ<sup>∗</sup> <sup>=</sup> <sup>6</sup> *<sup>μ</sup>* <sup>Ω</sup><sup>∗</sup> *pS <sup>R</sup>*<sup>2</sup> *c*2 ~21, thus indicating a moderate aerodynamic effect will assist with the generation of load capacity.

Figure 6 depicts the PGB force coefficients vs. shaft speed (Ω) and four *pS.* The coefficients are evaluated at a whirl frequency coinciding with shaft speed (*ω* = Ω), namely a synchronous speed condition. Note that at the centered condition, *KXX* = *KYY*, *KXY* = −*KYX*, *CXX* = *CYY*, and *CXY* = −*CYX*; see Equation (22). Effective force coefficients are

$$K\_{\rm eff} = K\_{XX} + C\_{XY} \,\omega \,\, \mathcal{C}\_{\rm eff} = C\_{XX} - K\_{XY} / \omega \,\,. \tag{34}$$

In general, the magnitude of the PGB force coefficients increases as the supply pressure grows. For hybrid operation (aerostatic plus aerodynamic effects), the direct stiffness (*KXX*) increases with rotor speed (Ω), whereas the direct damping *CXX* drops rapidly. At the top speed (30 krpm), *KXX* is approximately two times larger than at the low speed (aerostatic) condition. On the other hand, *CXX* reduces to ~1/5 in magnitude at low shaft speed. The cross-coupled stiffness *KXY* (<<*KXX*) is peculiar, as it reverses sign at a certain shaft speed. The larger *pS* is, the higher the speed at which *KYX* changes from positive to negative. Note −*CXY* << *CXX* with a peak magnitude at a certain shaft speed. Since *CXY* < 0, then *Keff* < *KXX* for all shaft speeds. At a low shaft speed condition, the effective damping is ~50% of *CXX*. On the other hand, at the high end of shaft speeds, *Ceff* is slightly larger than *CXX* since *KXY* < 0.

**Figure 6.** PGB (synchronous frequency) force coefficients vs. shaft speed: (**a**,**b**) stiffnesses (*KXX*, *KXY*), (**c**,**d**) damping coefficients (*CXX*, *CXY*), (**e**,**f**) effective stiffness (*Keff*) and damping (*Ceff*) coefficients. Excitation frequency *ω* = Ω (synchronous with speed). Supply pressure *pS* = 2 to 8 bar; MOS = 25 krpm.

#### **9. PGB Force Coefficients vs. Excitation Frequency**

Figure 7 displays the PGB force coefficients vs. frequency ratio (*ω/*Ω\*) for the bearing supplied with air at increasing pressures. In general, *KXX* grows (hardens) as the whirl frequency increases, whereas the direct damping (*CXX*) quickly drops. The magnitude of all force coefficients increases with an increase in gas supply pressure. In addition, at the cross-over frequency, *ω =* 1/2 Ω, half-frequency whirl, the effective stiffness (*Keff*) has a dip, and the effective damping turns positive, *Ceff* > 0. Hence, the PGB has the identical aerodynamic stability characteristics as a plain journal bearing. The predictions reveal that the cross-coupled stiffness (*KXY*) reverses sign for *ω* > 1/2 Ω, this frequency increasing as the supply pressure grows. Incidentally, for a low-pressure supply, *pS* < 4 bar, the direct damping is negative, *CXX* < 0, at low whirl frequencies (*ω <<* Ω\*). This effect is remarkable and a precursor to a pneumatic hammer; see [13].

**Figure 7.** PGB force coefficients vs. frequency ratio (*ω/*Ω). Mean operating speed Ω\* = 25 krpm (417 Hz); (**a**,**b**) stiffnesses (*KXX*, *KXY*), (**c**,**d**) damping coefficients (*CXX*, *CXY*), (**e**,**f**) effective stiffness (*Keff*) and damping (*Ceff*) coefficients. Supply pressure *pS* = 2 to 8 bar.

In sum, shaft whirl motions at super synchronous motions (*ω* > 1/2 Ω), produce stiffness hardening and a marked reduction in effective damping. On the other hand, SSVs are likely to occur at 1/2 frequency whirl since *Ceff* = 0.

#### **10. Stability of PGB**

\$ The stability of a point mass rotor (*Mr*) supported on one PGB is derived from \$ \$ **<sup>H</sup>**(*ω*) − *<sup>ω</sup>*2*Mr* **<sup>I</sup>** \$ \$ \$**<sup>=</sup>** 0 , where **<sup>I</sup>** is the 2 <sup>×</sup> 2 identity matrix. Following San Andrés [1], the instability threshold occurs at frequency <sup>s</sup> defined when the equivalent complex stiffness *He* has its imaginary part Im (*He*) **=** 0 and its real part > zero. The equivalent complex stiffness (*He*) is defined as

$$H\_{\mathfrak{E}\_{\{\varphi\}}} = \frac{1}{2} \left( H\_{XX\_{\{\varphi\}}} + H\_{YY\_{\{\varphi\}}} \right) - \left[ \frac{1}{4} \left( H\_{XX\_{\{\varphi\}}} - H\_{YY\_{\{\varphi\}}} \right)^2 + H\_{XY\_{\{\varphi\}}} H\_{YX\_{\{\varphi\}}} \right]^{1/2}.\tag{35}$$

Recall that at the centered condition (*e =* 0), *HXX* = *HYY*, *HXY* = −*HYX*. Hence, *H*<sup>e</sup> becomes

$$H\_{\mathfrak{C}\_{\langle \varphi \rangle}} = H\_{\text{XX}\_{\langle \varphi \rangle}} - i \, H\_{\text{XY}\_{\langle \varphi \rangle}} = (\mathcal{K}\_{\text{XX}} + \mathcal{a} \mathcal{C}\_{\text{XY}})\_{\langle \varphi \rangle} + i (\mathcal{a} \mathcal{C}\_{\text{XX}} - \mathcal{K}\_{\text{XY}})\_{\langle \varphi \rangle} = \mathcal{K}\_{\text{tf}/\langle \varphi \rangle} + i \, \mathcal{a} \mathcal{C}\_{\text{tf}/\langle \varphi \rangle} \tag{36}$$

$$\operatorname{Im}\left(H\_{\mathfrak{E}\_{\left(\varphi\right)}}\right) = 0 \to \mathbb{C}\_{\operatorname{eff}\_{\left(\varphi\right)}} = 0 \, ; \, \operatorname{Re}\left(H\_{\mathfrak{E}\_{\left(\varphi\right)}}\right) > 0 \to M\_{\operatorname{Tr}} = \frac{K\_{\operatorname{eff}\_{\left(\varphi\right)}}}{\mathfrak{Q}^2}.\tag{37}$$

When = <sup>1</sup> <sup>2</sup><sup>Ω</sup> <sup>→</sup> *Ceff*() <sup>=</sup> 0, and taking the operating speed as being the threshold speed of instability (Ω*<sup>T</sup>* = Ω\*=2), then the largest mass the rotating system can hold is *Mcr* = *Keff*( <sup>1</sup> <sup>2</sup> <sup>Ω</sup> *<sup>T</sup>*) ( 1 <sup>2</sup> <sup>Ω</sup> *<sup>T</sup>*) <sup>2</sup> . Table 2 lists the critical rotor mass for operation at various supply pressures and at MOS = 25 krpm. *Mcr* increases from 26.2 kg to 146.1 kg as *pS* = 2 bar → 8 bar. Note the weight of *Mcr* is much smaller than the PGB load capacity of ~(*K*<sup>Ω</sup> × *c*). Here*K*<sup>Ω</sup> = # *K*2 *XX* + *<sup>K</sup>*<sup>2</sup> *XY* is a hybrid mode operation stiffness evaluated at the MOS (Ω\*) and a zero whirl frequency (*ω* = 0).

**Table 2.** Aerostatic and hybrid stiffnesses for PGB at the centered condition and rotor speed Ω = 0 and 25 krpm. Air supply pressure *pS* = 2, 4, 6, and 8 bar. Clearance *c* = 0.010 mm (Λ*<sup>κ</sup>* = 5.3).


**11. PGB Load Capacity and Attitude Angle**

Predictions of the load carrying ability of the PGB for off-centered operation (*e* > 0) follow. Figure 8 depicts the load (*W*) vs. eccentricity ratio (*e*/*c*) for aerostatic operation (0 rpm) and at the mean MOS, Ω\* = 25 krpm, and increasing magnitudes of air supply pressure. The results from the FE computational model (1872 elements = 72 circumferential × 26 axial), shown with dark symbols, reveal *W* increases with *pS* and is proportional to journal eccentricity (*e*). The graphs include light-colored lines that represent an approximate load derived from the product of the (exact) static stiffness at *e =* 0 times the journal eccentricity, i.e., *Wapprox =* (*KS* × *e*) when Ω = 0 (left graph), and *Wapprox=* (*K*<sup>Ω</sup> × *e*) at the MOS. Table 2 lists the magnitudes of the analytical force coefficients for the aerostatic (Ω = 0) and hybrid (Ω\*) conditions.

**Figure 8.** FE model of PGB load (*W*) vs. journal eccentricity ratio (*e*/*c*) at 0 rpm and 25 krpm shaft speed. Air supply pressure *pS* = 2, 4, 6, and 8 bar; (**a**) Load at 0 rpm; (**b**) Load at 25 krpm. PGB clearance *c* = 0.010 mm. (Light colored) straight lines represent approximate load capacity derived from the exact solution at *e =* 0.

Note the remarkable agreement between the analytical solution and the FE numerical solution for operation with journal eccentricity (*e*) as large as 70% of the bearing clearance. In short, the load capacity of the cylindrical PGB is proportional to the journal displacement (*e*). Furthermore, for the hybrid mode operating condition, the attitude angle (*β*) between the eccentricity vector and the applied load vector remains constant. As seen in Table 2, *β* increases from 14.6◦ to 25.7◦ as the supply pressure increases from two bar to eight bar.

The analysis and predictions demonstrate that the PGB is a linear mechanical element. Presently, based on the results shown in Figure 8, an estimation for the bearing load capacity equals *<sup>W</sup> <sup>W</sup>*<sup>∗</sup> <sup>=</sup> *<sup>W</sup>* (*pS*−*pa* )*L D* ∼ ( 0.5 + 0.3 *pa pS* ) (1 + 0.13ΛΩ).

#### **12. An Example of Validation for the Static Performance of a PGB**

Li et al. [23] detailed an investigation of the static and dynamic performance of an (in-house constructed) cylindrical PGB. As in the current paper, the authors of [23] built a computational finite difference (FD) model for prediction of performance of simple PGBs and to make pertinent distinctions for operation as purely aerostatic (Ω = 0), aerodynamic (*pS = pa*), and hybrid (Ω > 0, *pS* > *pa*). Authors in [6] also detailed a handful of experimental results for the measurement of applied load vs. journal eccentricity for a PGB with two distinct clearances, small (16 μm) and large (31 μm). Table 3 details the geometry, porous layer physical properties, and the operating conditions of the test PGB in [23]. Note that the reference has many clerical errors including dubious captions in several of the figure captions. In addition, the direct force coefficients presented in [23] follow trends not in accordance with the theory, the direct stiffness not agreeing with one derived from the static load (*ω* → 0). The work also gives a cursory review of the cross-coupled force coefficients.

Figure 9 presents comparisons of the current analytical solution and FE model predictions vis-a-vis those in [23], experimental and numerical. As Figure 9a shows, the analytical solution and FE prediction of journal eccentricity (*e*) agrees well with the FD prediction [23] for *pS* = 4.7 bar, *c* = 31 μm and three static loads, 5 N to 15 N. For the top load of 15 N, the specific load *W*/(*LD*) = 0.1 bar; hence *W*/(*LD*)/(*pS* − *pa*) *=* 0.029, i.e., a very small load. The model predictions are slightly larger than the experimental results (ref. [23] does not report uncertainty or variability for the measured parameters nor variability for the force coefficients).

**Figure 9.** Predictions from current models (analytical and FE) compared to experimental results and predictions reported in [23]. (**a**) Journal eccentricity (*e*) vs. static load (*W*); (**b**) eccentricity (*e*) vs. shaft speed (Ω); (**c**) eccentricity (*e*) vs. supply pressure (*pS*).


**Table 3.** Operating conditions and dimensions of a cylindrical PGB reported in Ref. [23].

(\*) Assumed since [23] has many clerical mistakes, as confirmed by one of the authors [28]. Ref. [23] also uses gauge pressure to report *ps*.

As per the journal eccentricity (*e*) vs. shaft speed and a fixed load (*W =* 15 N) (see Figure 9b), the current analytical and FE model predictions agree well with the experimental result at the top speed of 24 krpm. For pure aerostatic operation (Ω = 0), the current model predictions are ~32% smaller than the experimental result. Figure 9c displays eccentricity (*e*) vs. supply pressure (*pS*) for operation at 24 krpm (ΛΩ = 0.097 at *pS* = 4.7 bar). For *pS* = 5.5 and 6.3 bar, the current predictions match the experimental values, the difference <2%. For *pS* = 4.7 bar, the models produce a slightly larger (*e*) than the measured one.

Comparisons of predictions of the experimental force coefficients are omitted. As frequency grows, the test coefficients follow differing paths and have lower physical magnitudes than the model predictions, current and those in [23]. Feng [28] attributes the difference to the absence of a surface restrictive layer in the constructed PGB. Otsu et al. [18] explain difference in the performance of PGBs with and without the said restrictive layer. Interestingly enough, the current predictive model as well as that in [23] do not apply to PGBs with restrictive layers.

#### **13. Conclusions**

Porous surface gas bearings (PGBs) have come of age to enable high-speed, near friction free rotating machinery with improved reliability and availability. In the last decade (2012 and onward), numerous computational analyses for calculation of PGB forced performance, static and dynamic, have appeared. Most analyses, however, tackled the flow physical equations by computer while ignoring the vast body of archival literature on the subject.

The paper reviews the archival literature on PGBs and reproduces, using modest analytical means, an exact solution to the flow field and forced performance of a cylindrical PGB. This solution has been available since 1979 [5].

Predictions from the analysis serve to validate the results of a finite element computational model. Further comparisons to recent experimental results validate the analysis results and serve to elucidate the effect of external pressure and shaft speed on PGB performance. The learning from this work is as follows:


achieve when also considering manufacturing costs and devising procedures for easy installation.


The authors hope other analysists, engineering students in particular, appreciate the effort to obtain the solution, the one that exactly predicts bearing performance and allows quick assessment of trends and selection of physical parameters. In a world ruled by complex numerical models attacking physics with computers, fundamental mathematics does bring meaning to an engineering endeavor as well as personal solace.

**Author Contributions:** Conceptualization, L.S.A. and A.D.; Formal analysis, L.S.A.; Funding acquisition, L.S.A.; Investigation, J.Y.; Project administration, L.S.A.; Resources, L.S.A. and A.D.; Software, J.Y.; Supervision, L.S.A.; Validation, J.Y.; Visualization, L.S.A. and J.Y.; Writing—original draft, L.S.A.; Writing—review & editing, L.S.A. and J.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was funded by Texas A&M University-CONACYT Collaborative Research Program, and Texas A&M University J. Mike Walker '66 Department of Mechanical Engineering 2020 Seed Grant.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Thanks to the Texas A&M University-CONACYT Collaborative Research Program, Texas A&M University J. Mike Walker '66 Department of Mechanical Engineering 2020 Seed Grant, for partial financial support of the research. Thanks to NewWay® Air Bearings for donating the bearing.

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

#### **Nomenclature**


#### **References**


## *Article* **Parameter Sensitivity Analysis on Dynamic Coefficients of Partial Arc Annular-Thrust Aerostatic Porous Journal Bearings**

**Pyung Hwang 1,\*, Polina Khan <sup>2</sup> and Seok-Won Kang 1,\***


**Abstract:** Aerostatic bearings are widely used in high-precision devices. Partial arc annular-thrust aerostatic porous journal bearings are a prominent type of aerostatic bearings, which carry both radial and axial loads and provide high load-carrying capacity, low air consumption, and relatively low cost. Spindle shaft tilting is a resource-demanding challenge in numerical modeling because it involves a 3D air flow. In this study, the air flow problem was solved using a COMSOL software, and the dynamic coefficients for tilting degrees of freedom were obtained using finite differences. The obtained results exhibit significant coupling between the tilting motion in the x-and y-directions: cross-coupled coefficients can achieve 20% of the direct coefficient for stiffness and 50% for damping. In addition, a nonlinear behavior can be expected, because the tilting motion within 3◦, tilting velocities within 0.0012◦/s, and relative eccentricity of 0.2 have effects as large as 20% for direct stiffness and 100% for cross-coupled stiffness and damping. All dynamic coefficients were fitted with a polynomial of eccentricity, tilting, and tilting velocities in two directions, with a total of six parameters. The resulting fitting coefficient tables can be employed for the fast dynamic simulation of the rotor shaft carried on the proposed bearing type.

**Keywords:** aerostatic journal bearings; porous bearing; annular-thrust bearing; numerical simulation; finite element method; dynamic coefficients; eccentricity; tilting; tilting velocity; partial arc

#### **1. Introduction**

Porous aerostatic bearings are widely used in machinery, owing to their relatively high load-carrying capacity and low air consumption [1–3]. Similar to other types of aerostatic bearings, owing to their high accuracy and zero contamination, they are suitable for high-precision devices. In addition, the friction drag of a porous aerostatic bearing is low [4,5], its motion accuracy is higher than that of a conventional orifice-type aerostatic bearing [5,6], and its load thresholds at high rotation speeds are the highest among other types of aerostatic bearings [5,7]. The production of porous materials with desirable permeabilities is relatively complicated. However, the amount of porous material required can be reduced via design optimization [8].

The optimal governing equations required for porous air bearing analysis remain under discussion. Zhong et al. [9] experimentally determined Ergun's equation coefficients for the pressure drop in sintered metal porous media for air bearings. Belforte et al. [10] experimentally verified Forchheimer's law for porous air bearing applications. Zhong et al. [11] verified the accuracy of Forchheimer's law for air flow through sintered metal porous media and Darcy's law, under slight pressure drops. Recently, it has been reported that a deep learning-based approach can be employed for high-precision and low-cost analysis of the 3D flow state inside a porous material [12].

Various numerical methods and engineering packages have been adopted for air pressure calculations in porous bearings. Huang et al. [13] analyzed the pressure in a porous conveyor air bearing using the FLUENT software. Cui et al. [14] used FLUENT

**Citation:** Hwang, P.; Khan, P.; Kang, S.-W. Parameter Sensitivity Analysis on Dynamic Coefficients of Partial Arc Annular-Thrust Aerostatic Porous Journal Bearings. *Appl. Sci.* **2021**, *11*, 10791. https://doi.org/ 10.3390/app112210791

Academic Editors: Terenziano Raparelli, Federico Colombo, Andrea Trivella and Luigi Lentini

Received: 30 September 2021 Accepted: 11 November 2021 Published: 15 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

software to numerically analyze the effect of manufacturing errors on the load-carrying capacity and stiffness of an annular-thrust porous aerostatic bearing. In addition, the deformation of porous thrust bearings was numerically investigated based on the fluid– solid interaction method by Wang et al. [15]. Van Ostayen et al. [16] used the COMSOL software for active aerostatic bearing analysis. Hwang and Khan [17] conducted 3D finite element analysis, while Plante et al. [18] proposed the shed 1D generalized flow theory. Otsu et al. [19] and Hosokawa et al. [8] used the FDM.

Several papers [19–22] deal with the numerical calculation of stiffness and damping coefficients corresponding to eccentricity or parallel shaft displacement (Figure 1a). However, insufficient data are available on stiffness and damping for shaft tilting (Figure 1b) in porous air bushing. Cui et al. [14,23] analyzed the non-flatness effect of the air-bearing surface on the stiffness and damping of an aerostatic porous bearing, which is computationally similar to modeling the tilting.

**Figure 1.** Two spindle motion modes: (**a**) eccentricity; (**b**) tilting. sb and sr denote the bearing and rotor axes, respectively.

The rotor in ultraprecision devices is usually supported by journal and thrust bearings, which share a common air supply system but separate porous pads [14,23]. Annularthrust bearings share the same porous pad carrying loads in both the annular and thrust directions (Figure 2). This dual role decreases the manufacturing cost while triggering flow coupling [23–27]. A partial arc bearing has a higher load-carrying capacity but lower stiffness than a full-arc bearing [27].

**Figure 2.** Annular-thrust porous aerostatic bearing.

The numerical modeling of the rotor motion, including the effect of shaft tilting, as well as the modeling of a partial arc annular-thrust porous aerostatic bearing is substantially time-consuming, because it requires a 3D air flow solution. Dynamic modeling can be performed in a less resource-consuming manner if the full bearing model is replaced with stiffness and damping coefficients. The nonlinear dynamic motion can be modeled, provided the influences of eccentricity, tilting, and velocities on the dynamic coefficients are included. In this study, the air flow in the porous pad and air lubrication film between the porous pad and solid surface of the shaft were numerically modeled using the COMSOL Multiphysics® software. The pressure distribution and total forces and moments were calculated for low eccentricity, low tilting angles, and tilting speeds. The stiffness and damping coefficients were determined using finite differences.

#### **2. Materials and Methods**

The porous air bearing was modeled using the COMSOL Multiphysics software, "Darcy's law", and "Thin-Film Flow, Shell" modules. Here, the static solution is sufficient in determining the response of the displacement and velocity of the spindle. In addition, the porous bushing is a hollow cylinder parallel to the z-axis with a cut sector *α*. Figure 3 illustrates the different boundary conditions.

**Figure 3.** Porous pad boundaries: (**a**) airbearing surface carrying both radial and axial loads; (**b**) air supply surface; (**c**) outflow to atmospheric pressure; (**d**) outflow to atmospheric pressure in *xy*-plane with axes and bearing partial arc.

The pressure in the porous pad is modeled with the "Darcy's law" module:

$$-\nabla \left(\rho \frac{\kappa}{\mu} \nabla p\right) = 0\tag{1}$$

where *p*, *ρ*, *κ*, and *μ* denote the air pressure, air density, permeability coefficient, and dynamic viscosity of air, respectively. ∇ is the gradient operator in three-dimensional space. The pressure at the air-bearing surface (see Figure 3a) is modeled with the "thin-film flow, shell" module:

$$\nabla\_t \left( \rho h \left( \frac{1}{2} v\_t - \frac{h^3}{12\mu} \nabla\_t p\_f \right) \right) = -\rho \frac{\partial}{\partial t} h + Q \tag{2}$$

where *pf* , *h*, *vt*, and *Q* represent the air film pressure, air film thickness, tangential spindle surface velocity, and air flux, respectively. ∇<sup>t</sup> is the gradient operator in the twodimensional space of the air film surface and ∇ is the partial derivative operator. At the interface between the porous pad and the air-bearing surface, the pressure and air flux must be continuous. In the "Darcy's law" module, the "Pressure" boundary condition is adopted at the surfaces that open to the air lubrication film, which is expressed as:

$$p = p\_f.\tag{3}$$

In the "Thin-Film Flow, Shell", the continuity of the air flux is expressed by admittance Y in the "Perforation" interface, which is expressed as:

$$\mathcal{Y} = \frac{\mathcal{Q} \cdot \mathfrak{n}}{\rho \ p\_f} \tag{4}$$

where *n* = *nx*, *ny*, *nz* is the air-bearing surface normal vector. The air film thickness in the "Thin-Film Flow, Shell" module depends on the spindle eccentricity (*εx*, *εy*), tilting (*tx*, *ty*), and tilting velocities (. *tx*, . *ty*) (refer to Figure 1):

$$
\hbar = \mathfrak{c} + \mathfrak{u} \cdot \mathfrak{n} - \mathfrak{u} \cdot \nabla\_t \mathfrak{h} \tag{5}
$$

$$\frac{\partial}{\partial t}h = \boldsymbol{\sigma} \cdot \boldsymbol{n} - \boldsymbol{\sigma} \cdot \nabla\_t h \tag{6}$$

$$\mathbf{u} = \left( c\varepsilon\_x + t\_x(z - L/2); c\varepsilon\_y + t\_y(z - L/2); -t\_x\mathbf{x} - t\_y\mathbf{y} \right) \tag{7}$$

$$\mathbf{w} = \left( -y\omega + \dot{t}\_x (z - L/2); \mathbf{x}\omega + \dot{t}\_y (z - L/2); -\dot{t}\_x \mathbf{x} - \dot{t}\_y y \right) \tag{8}$$

where *u* is the spindle surface displacement, and *v* is the spindle surface velocity.

The isothermal compressibility of the air was modeled using the user-defined density in both modules. Here, the density is related to the air pressure in the "Variables" block as follows:

$$
\rho = \frac{\rho\_0 \, p}{p\_0},
\tag{9}
$$

where the subscript 0 refers to the atmospheric conditions, and *p*<sup>0</sup> = 1 bar. The air supply is modeled by another "Pressure" interface in the "Darcy's law" module with a constant pressure of *ps* (see Figure 3b). The outflow at the air-bearing edges is given by the "Border" interface in the "Thin-Film Flow, Shell" module with a fixed pressure of *p*<sup>0</sup> (see Figure 3c). The other boundaries were sealed.

After the static problem was solved and the air pressure was calculated, the moments in the *x* and *y* tilt directions were calculated, using the "Surface Integration" in the "Derived Values" block, according to the following expressions:

$$M\_x = \iint p(n\_x(z - L/2) - n\_z x) dx \, dy \, \tag{10}$$

$$M\_{\mathcal{Y}} = \iint p \left( n\_{\mathcal{Y}} (z - L/2) - n\_{z} \mathcal{y} \right) dx \prime dy \prime \tag{11}$$

where *x* and *y* are the local coordinates of the air film surface.

Stiffness and damping for tilting degrees of freedom were calculated by using finite differences:

$$K\_{\rm tx,tx} = -\left(M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x} + \Delta t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right) - M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right)\right) / \Delta t\_{\rm x} \tag{12}$$

$$K\_{\rm tx,ty} = -\left(M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y} + \Delta t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right) - M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right)\right) / \Delta t\_{\rm y} \tag{13}$$

$$K\_{\rm ty,tx} = -\left(M\_y\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x} + \Delta t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right) - M\_y\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right)\right) / \Delta t\_{\rm x} \tag{14}$$

$$K\_{\rm ty.ty} = -\left(M\_y\left(\varepsilon\_{x\prime}\varepsilon\_{y\prime}, t\_{x\prime}t\_y + \Delta t\_{y\prime}\dot{t}\_{x\prime}\dot{t}\_y\right) - M\_y\left(\varepsilon\_{x\prime}\varepsilon\_{y\prime}t\_{x\prime}t\_y, \dot{t}\_{y\prime}\dot{t}\_{x\prime}\dot{t}\_y\right)\right) / \Delta t\_y \tag{15}$$

$$D\_{\rm tx,tx} = -\left(M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x} + \Delta \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right) - M\_{\rm x}\left(\varepsilon\_{\rm x}, \varepsilon\_{\rm y}, t\_{\rm x}, t\_{\rm y}, \dot{t}\_{\rm x}, \dot{t}\_{\rm y}\right)\right) / \Delta \dot{t}\_{\rm x} \tag{16}$$

$$D\_{tx,ty} = -\left(M\_{\mathbf{x}}\left(\varepsilon\_{x\_{\prime}}\varepsilon\_{y\_{\prime}}t\_{\mathbf{x}\_{\prime}}t\_{\mathbf{y}\_{\prime}}\dot{t}\_{\mathbf{x}\_{\prime}}\dot{t}\_{\mathbf{y}} + \Delta\dot{t}\_{\mathbf{y}}\right) - M\_{\mathbf{x}}\left(\varepsilon\_{x\_{\prime}}\varepsilon\_{y\_{\prime}}t\_{\mathbf{x}\_{\prime}}t\_{\mathbf{y}\_{\prime}}\dot{t}\_{\mathbf{x}\_{\prime}}\dot{t}\_{\mathbf{y}}\right)\right) / \Delta\dot{t}\_{\mathbf{y}}\tag{17}$$

$$D\_{\mathcal{Y},tx} = -\left(M\_{\mathcal{Y}}\left(\varepsilon\_{x\_{\prime}}\varepsilon\_{y\_{\prime}}t\_{x\_{\prime}}t\_{y\_{\prime}}\dot{t}\_{x} + \Delta\dot{t}\_{x\_{\prime}}\dot{t}\_{y}\right) - M\_{\mathcal{Y}}\left(\varepsilon\_{x\_{\prime}}\varepsilon\_{y\_{\prime}}t\_{x\_{\prime}}t\_{y\_{\prime}}\dot{t}\_{x\_{\prime}}\dot{t}\_{y}\right)\right) / \Delta\dot{t}\_{x} \tag{18}$$

$$D\_{ty,ty} = -\left(M\_y\left(\varepsilon\_{x\cdot}\varepsilon\_{y\cdot}t\_{x\cdot}t\_y, \dot{t}\_{x\cdot}\dot{t}\_{y} + \Delta\dot{t}\_y\right) - M\_y\left(\varepsilon\_{x\cdot}\varepsilon\_{y\cdot}t\_{x\cdot}t\_{y\cdot}\dot{t}\_{x\cdot}\dot{t}\_y\right)\right) / \Delta\dot{t}\_y\tag{19}$$

#### **3. Results**

The airbearing moments, tilting stiffness, and damping were calculated for *ε<sup>x</sup>* = 0.1, 0.2, 0.3 , *ε<sup>y</sup>* = 0.0, 0.1, *tx* = −0.003◦, −0.00275◦, ... , 0.003◦, and *ty* = 0.0◦, 0.00025◦,. . . , 0.003◦ , respectively; with *ps* = 8 bar, bearing length *L* = 0.11 m, shaft diameter *D* = 0.1 m, porous cylinder thickness *b* = 0.01 m, clearance *c* = 10 μm, rotation speed = 1000 rpm, bearing arc <sup>α</sup> = 240◦, porosity = 0.4, and permeability = 8.33 × <sup>10</sup>−<sup>17</sup> <sup>m</sup>2. The film variation due to tilting in the given range does not exceed 57% for simultaneous change of *tx* and *ty*. Together with 30% variation due to eccentricity, it gives a minimal film thickness of 13% of clearance, or 1.3 μm. The air-bearing pressure, forces, and moments were calculated for all combinations of parameters *tx*, *ty*, *εx*, and *εy*, and their combinations with the following pairs of tilting velocity values: . *tx*, . *ty* ∈ , (0, 0), Δ . *tx*, 0 , 2Δ . *tx*, 0 , 0, Δ . *ty* , 0, 2Δ . *ty* -. The total number of cases was 13000. The steps in the differences (12–19) were Δ*tx* = Δ*ty* = 0.00025◦, and Δ . *tx* = Δ . *ty* = 0.00573◦/s. The maximum and average absolute errors of the finite differences are presented in Table 1.

**Table 1.** Finite difference errors.


#### *3.1. Preliminary Analysis*

As illustrated in Figure 4, the direct stiffness coefficient *Ktx*,*tx* has values in the range of 80–140 kN·m/rad. Within the range under consideration, *Ktx*,*tx* depends linearly on *ε<sup>x</sup>* and nonlinearly on *tx*, *ty*, . *tx*, . *ty*. The parameters *tx*, . *tx*, and *ε<sup>x</sup>* exhibit significant effects (more than 10%). *Ktx*,*tx* is almost insensitive to *tx* at low negative tilting angles; however, it decreases linearly with *tx* at positive and high negative tilting angles (see Figure 4a). The effect of the x-tilting velocity . *tx* is the diminishing of *Ktx*,*tx*, which is more significant at negative tilting angles. *Ktx*,*tx* increases with decreasing *ty*, increasing . *ty* and *εx*, and is insensible to *ε<sup>y</sup>* (see Figure 4b,c).

**Figure 4.** Parameter sensitivity of the direct stiffness coefficient *Ktx*,*tx*: (**a**) Effects of the tilting angle and tilting velocity in the x-direction; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

The cross-coupled stiffness coefficient *Ktx*,*ty* is approximately 10 times smaller than the direct stiffness coefficient *Ktx*,*tx*, and lies in the range of 6–15 kN·m/rad (see Figure 5). It is highly sensitive to all parameters, except for . *tx*. *Ktx*,*ty* increases with decreasing *tx* and . *ty*, and increasing *ty* (see Figure 5a,b). Eccentricity *ε<sup>x</sup>* can increase or decrease *Ktx*,*ty*, depending on *ε<sup>y</sup>* (see Figure 5c).

**Figure 5.** Parameter sensitivity of the cross-coupled stiffness coefficient *Ktx*,*ty*: (**a**) Effects of tilting angle and tilting velocity in the x-direction; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

*Kty*,*tx*, which is another cross-coupled stiffness coefficient, is of the same order of magnitude; however, it is negative, and falls in the range of −18–0 kN·m/rad (refer to Figure 6). It is not sensitive to . *tx*, and the . *ty* effect on it is also negligible (see Figure 6a,b). The magnitude of *Kty*,*tx* is higher at lower *tx* and *ε<sup>y</sup>* and at higher *ty* and *εx*.

**Figure 6.** Parameter sensitivity of the cross-coupled stiffness coefficient *Kty*,*tx*: (**a**) Effects of the tilting angle and tilting velocity in the x-direction; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

The direct stiffness coefficient *Kty*,*ty* is higher than that of *Ktx*,*tx*. *Kty*,*ty* falls in the range of 110–140 kN·m/rad (refer to Figure 7). This is triggered by the cut arc, which increases force and decreases the stiffness in the x-direction [27]. The y-tilting angle *ty* exerts a strong nonlinear effect on *Kty*,*ty*. The stiffness coefficient *Kty*,*ty* is almost constant for a low positive *ty*, and then increases with increasing *ty*. In contrast to *Ktx*,*tx*, which is primarily sensitive to "its own" parameter *tx* and significantly less sensitive to *ty*. *Kty*,*ty* is equally sensitive to both *tx* and *ty*. In addition, it is equally sensitive to both higher *<sup>ε</sup><sup>x</sup>* and *<sup>ε</sup>y*, with a negligible . *ty* effect, and no . *tx* effect. Furthermore, it is higher for lower *tx* and higher *ty*, *εx*, and *εy*.

**Figure 7.** Parameter sensitivity of the direct stiffness coefficient *Kty*,*ty*: (**a**) Effects of the tilting angle and tilting velocity in the x-direction; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

Because the simultaneous variation of tilting velocities was not studied, the effect of . *tx* was solely studied for *Dtx*,*tx* and *Dty*,*tx*, while the effect of . *ty* was solely studied for *Dtx*,*ty* and *Dty*,*ty*. Instead of the missing parameter, the coupling between *tx* and *ty* is comprehensively presented in Figures 8–11.

The direct damping coefficient *Dtx*,*tx* is positive and falls in the range 0.2–2.2 kN·m s/rad (refer to Figure 8). It is higher for lower *tx*, and higher *ty* and higher . *tx* (see Figure 8a,b). The effect of eccentricity is also substantial, nonlinear, and strongly coupled between *ε<sup>x</sup>* and *ε<sup>y</sup>* (see Figure 8c).

**Figure 8.** Parameter sensitivity of the direct damping coefficient *Dtx*,*tx*: (**a**) Effects of the tilting angle and tilting velocity in the x-direction; (**b**) Effects of the x and y tilting angles; (**c**) Effects of the relative eccentricities. All unreferenced parameters have zero values.

The cross-coupled damping coefficient *Dtx*,*ty* does not exceed 1 kN m s/rad in magnitude, being mostly negative (see Figure 9). *Dtx*,*ty* rapidly decreases with decreasing *tx* or increasing *ty* in the negative *tx* range, while in the positive *tx* range, it is much less sensible both to *tx* and *ty* (see Figure 9a). *Dtx*,*ty* decreases with increasing *ty* and is almost insensible to . *ty* (see Figure 9b). The dependence of *Dtx*,*ty* on eccentricity parameters is highly nonlinear. *Dtx*,*ty* can take relatively high positive values for zero *εx*. and positive *εy*, but in most of the range, the sensitivity to *ε<sup>x</sup>* and *ε<sup>y</sup>* is low (see Figure 9c).

**Figure 9.** Parameter sensitivity of the cross-coupled damping coefficient *Dtx*,*ty*: (**a**) Effects of x and y tilting angles; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

The cross-coupled damping coefficient *Dty*,*tx* falls in the range −0.6–0.13 kN m s/rad (see Figure 10). It is negative for positive *ty* and negative *tx*, otherwise it is positive. Similar to *Dtx*,*ty*, *Dty*,*tx* rapidly decreases with decreasing *tx* or increasing *ty* in the negative *tx* range, while in the positive *tx* range, it is much less sensible both to *tx* and *ty*. However, for *tx* > 0.002 ◦, *Dty*,*tx* increases with increasing *ty* (see Figure 10b), while *Dtx*,*ty* decreases with increasing *ty* for any *tx* (see Figure 9a). For zero *ty*, *Dty*,*tx* decreases with increasing *tx* and decreasing . *tx*. Similar to all previously considered damping coefficients, the dependence of *Dty*,*tx* on the eccentricity parameters is highly nonlinear.

**Figure 10.** Parameter sensitivity of the cross-coupled damping coefficient *Dty*,*tx*: (**a**) Effects of the tilting angle and tilting velocity in the x-direction; (**b**) Effects of x and y tilting angles; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

In contrast with all other damping coefficients, the direct damping coefficient *Dty*,*ty* does not significantly depend on *ε<sup>x</sup>* and *ε<sup>y</sup>* (refer to Figure 11). It falls in the range 0.6–2.0 kN·m s/rad. It is higher for lower *tx*, and higher *ty*, similar to *Dtx*,*tx*. *Dty*,*ty* is almost insensible to the y-tilt velocity . *ty*. It slightly increases with increasing *ε<sup>x</sup>* and increasing *εy*.

**Figure 11.** Parameter sensitivity of the direct damping coefficient *Dty*,*ty*: (**a**) Effects of x and y tilting angles; (**b**) Effects of the tilting angle and tilting velocity in the y-direction; (**c**) Effects of relative eccentricities. All unreferenced parameters have zero values.

In summary, the direct stiffness and damping coefficients were positive. The crosscoupled stiffness coefficients are approximately 10 times smaller than the direct stiffness coefficients, and the cross-coupled damping coefficients are approximately 2–3 times smaller than the direct damping coefficients. Hence, a stable tilting motion was expected. However, information on the rotor inertia and the bearing at the other end is required to make precise conclusions on the dynamic behavior. Moreover, the results above exhibit a strong coupling between the eccentricity and tilting motion, including the nonlinearity of the moments relative to the position parameters. This can make the rotor motion severely complicated.

#### *3.2. Regression*

Although it is possible to obtain the direct solution of the time-dependent problem in the model described above, the process is significantly time consuming because the air flow problem must be solved at every step. Conversely, the air-bearing forces and moments can be calculated using the interpolated stiffness and damping coefficients. In this study, each dynamic coefficient for each combination of *εx*, *εy*, . *tx*, and . *ty* was first fitted with a quadratic polynomial of *tx* and *ty*, and then each coefficient of the polynomial was fitted with a linear function of *εx*, *εy*, . *tx*, and . *ty*. The accuracy of the fitting is expressed in the form of determination coefficients, as presented in Table 2. The accuracy is relatively good, although it can be improved for *Ktx*,*ty*, *Dtx*,*tx*, *Dtx*,*ty*, and *Dty*,*tx* by providing more points and increasing the fitting order in *ε<sup>x</sup>* and *ε<sup>y</sup>* degrees of freedom.

**Table 2.** Coefficients of determination for dimensionless dynamic coefficient fitting.


The fitting results are presented in Tables 3–10. An example of the table interpretation is given below for *Ktx*,*tx*:

$$\begin{aligned} \check{K}\_{tx,tx} &= \left(1.09533 + 1.49129\epsilon\_x + 0.17719\epsilon\_y - 0.02161\dot{T}\_x\right) + \\ &\quad \left(-0.30745 - 3.94666\epsilon\_x + 0.10915\dot{T}\_x + 0.00478\dot{T}\_y\right)T\_x + \\ &\quad \left(-0.15003 + 0.09324\epsilon\_x - 1.16823\epsilon\_y + 0.01582\dot{T}\_x + 0.0114\dot{T}\_y\right)T\_y + \\ &\quad \left(-0.3232 + 2.45267\epsilon\_x + 2.42515\epsilon\_y + 0.16265\dot{T}\_x - 0.28434\dot{T}\_y\right)T\_x T\_y + \\ &\quad \left(7.99445\epsilon\_x - 0.26854\epsilon\_y - 0.35236\dot{T}\_x + 0.14364\dot{T}\_y\right)T\_x^2 + \\ &\quad \left(0.34256 + 0.41188\epsilon\_x + 1.09915\epsilon\_y - 0.13117\dot{T}\_x + 0.17535\dot{T}\_y\right)T\_y^2 \end{aligned}$$

**Table 3.** Fitting results for *Ktx*,*tx*.


**Table 4.** Fitting results for *Ktx*,*ty*.




**Table 6.** Fitting results for *Kty*,*ty*.


**Table 7.** Fitting results for *Dtx*,*tx*.


**Table 8.** Fitting results for *Dtx*,*ty*.


**Table 9.** Fitting results for *Dty*,*tx*.




Here, the dimensionless variables are adopted:

$$\frac{K\_{i,j}}{\overline{K}\_{i,j}} = \frac{p\_0 D L^3}{16 \text{ c}} = 83,187 \frac{\text{N m}}{\text{rad}} \tag{20}$$

$$\frac{D\_{i,j}}{\bar{D}\_{i,j}} = \frac{p\_0 D L^2}{2 \left(0.0001 \frac{\text{rad}}{\text{s}}\right)} = 1.51 \times 10^5 \frac{\text{N m s}}{\text{rad}}\tag{21}$$

$$\frac{t\_x}{T\_x} = \frac{t\_y}{T\_y} = \frac{c}{L} = 0.000182 \text{ rad} = 0.0104^\circ \tag{22}$$

$$\frac{\dot{t}\_x}{\dot{T}\_x} = \frac{\dot{t}\_y}{\dot{T}\_y} = 0.0001 \frac{\text{rad}}{\text{s}} = 0.00573^\circ \frac{1}{\text{s}} \tag{23}$$

#### **4. Discussion**

The tilting stiffness and damping analyzed in this study could be significant for the spindle shaft tilting motion, provided the shaft length is not significantly longer than the bearing length. Otherwise, the shaft tilting will be controlled by the air-bearing force applied at a long distance from the shaft center, while the tilting stiffness and damping will influence the shaft bending. In addition, bending can also significantly influence the air-bearing forces, moments, and dynamic coefficients, because of the convexity and concavity of the bearing surface [14,23]. Based on the obtained results, the tilting motion is independently stable because the direct stiffness and damping coefficients are positive and significantly higher in amplitude than the cross-coupled coefficients. However, the tilting dynamic coefficients were substantially influenced by shaft eccentricity. Hence, the stability of the shaft parallel and tilting motions must be analyzed together. The increased rotor speed would result in a more pronounced hydrodynamic component in the bearing forces, moments, stiffness, and damping components, as well as in more significant x-y coupling, that is typical for the journal bearings.

#### **5. Conclusions**

The problem of tilting stiffness and damping analysis is hardly ever addressed in the literature since the fundamental books on rotodynamic and tribology where the theoretical analysis was presented for idealized types of bearings. This ignorance can be explained by the minor role of individual bearing moments, compared to the total shaft-bearing systems moments in most applications. However, the tilting moments exist, and they become more and more important in high-precision machinery. Besides the bearing moments, the tilting motion greatly affects the bearing forces, as well as other types of air film distortions, which are widely considered in recent publications.

The purpose of this study was to predict the behavioral characteristics according to the design and operating parameters of porous aerostatic bearings that can be utilized in applications requiring high-precision positioning, such as semiconductor manufacturing equipment. Most importantly, the highly efficient and reliable numerical analysis approach proposed in this study can contribute to increasing system stability and reducing costs and risks at the product level prior to the actual manufacturing process.

**Author Contributions:** Conceptualization, P.H. and P.K.; methodology, P.K.; software, P.K.; validation, P.H. and P.K.; formal analysis, P.K.; investigation, P.H. and P.K.; resources, P.H. and P.K.; data curation, P.H.; writing—original draft preparation. P.K.; writing—review and editing. P.H. and S.-W.K.; visualization, P.K.; supervision, P.H.; project administration, P.H.; funding acquisition, P.H. and S.-W.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partly supported by Korea Institute for Advancement of Technology (KIAT) grant funded by the Korea Government (MOTIE) (P0002092, The Competency Development Program for Industry Specialist) and Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea Government (MOTIE) (20214000000010, Gyeongbuk Regional Wind Energy Cluster Human Resources Development Project) and State Assignment Project (No. FWEU-2021-0005) of the Fundamental Research Program of the Russian Federation 2021–2030.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available by request.

**Acknowledgments:** The research was performed with use of the equipment of High-Temperature Circuit Multi-Access Research Center (Project No. 13.ЦKΠ.21.0038).

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

#### **References**


## *Article* **Experimental and Numerical Dynamic Identification of an Aerostatic Electro-Spindle**

**Federico Colombo \*, Luigi Lentini, Andrea Trivella, Terenziano Raparelli and Vladimir Viktorov**

Politecnico di Torino, Department of Mechanical and Aerospace Engineering, Corso Duca degli Abruzzi 24, 10129 Turin, Italy; luigi.lentini@polito.it (L.L.); andrea.trivella@polito.it (A.T.); terenziano.raparelli@polito.it (T.R.); vladimir.viktorov@formerfaculty.polito.it (V.V.)

**\*** Correspondence: federico.colombo@polito.it

**Featured Application: Aerostatic spindles are used in high-speed micromachining applications. The main goal of the work is to validate the developed non-linear numerical model through the proposed identification procedure and the performed experimental tests.**

**Abstract:** This paper proposes a method to experimentally identify the main modal parameters, i.e., natural frequencies and damping ratios, of an aerostatic spindle for printed board circuit drilling. A variety of methods is applied to the impulse-response function of the spindle in the presence of zero rotational speed and different supply pressures. Moreover, the paper describes the nonlinear numerical model of the spindle, which consists of a four-degree-of-freedom (DOF) rigid and unsymmetrical rotor supported by two aerostatic bearings. The main goal of the work is to validate the developed non-linear numerical model through the proposed identification procedure and the performed experimental tests. The comparison proves satisfactory, and the possible sources of uncertainty are conjectured.

**Keywords:** aerostatic journal bearings; PCB spindles; dynamic identification; logarithmic decrement method; complex exponential method; half-power method; damping factor identification; rigid rotor; modal analysis

#### **1. Introduction**

Gas bearings, owing to their properties, play an essential role in turbochargers, gyroscopes, turbo blowers, gas turbines and micromachining spindles. In particular, micromachining spindles are widely used in surface finishing, printed circuit board (PCB) drilling, micro-processing, micro-milling, micro-grinding and, in general, for micromachining materials with low shear resistance. This is because compared to oil and rolling bearings, gas bearings can reach higher rotation speeds while simultaneously reducing the amount of machine maintenance required.

The aim of research on spindles supported by gas bearings is to develop systems with excellent dynamic stability and stiffness, along with very low runout of the tool, even at very high rotation speeds. Poor damping is one of the main disadvantages of journal gas bearings, which can compromise the accuracy of high-precision machining. Moreover, friction power losses in motors and bearings can be a non-negligible source of heat at high speeds. The contribution of thermally induced deformations of a machine tool can exceed 50% of the total machining error [1], which is very significant in precision milling. To reduce these problems, an efficient and precise refrigeration system is required. Recent developments in manufacturing technologies have made it possible to realize innovative high-precision devices that can even be miniaturized. Research has therefore been enriched by numerous experimental contributions, some of which are discussed below. Machine tools of 4 mm and 1/8" diameters used for micro-milling can now be supported on aerostatic bearings up to 400 krpm and 450 krpm [2,3]. Experiments on

**Citation:** Colombo, F.; Lentini, L.; Trivella, A.; Raparelli, T.; Viktorov, V. Experimental and Numerical Dynamic Identification of an Aerostatic Electro-Spindle. *Appl. Sci.* **2021**, *11*, 11462. https://doi.org/ 10.3390/app112311462

Academic Editor: Homer Rahnejat

Received: 8 October 2021 Accepted: 29 November 2021 Published: 3 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

dynamic characterization of air-bearing spindles are described in [4,5], where the frequencyresponse function (FRF) of the spindle is measured. Impact tests have been carried out with an impact hammer [4] or with a custom-designed impact excitation system able to provide excitations up to 20 kHz [5]. Ref. [6] describes a test rig used to measure the dynamic stiffness and damping coefficients of a hemisphere spiral groove hybrid gas bearing. In this rig, the excitation is provided by means of two electromagnetic exciters.

Mathematical models have long been used as a valuable tool to investigate the dynamic performance and stability of high-speed rotors. Early attempts at modelling this kind of system were aimed at gaining a better insight into the half-whirl instability, known today as asynchronous whirl instability [7]. At the time, this was one of the most troublesome problems in journal bearings. The first attempts to actually solve the equations of motion and the gas-lubricated equation simultaneously were reported by Sternlicht [8], Poritsky and Arwas [9]. In their solution, they linearized the equations of motion of a 2 DOF rotor, along with the Reynolds equation, to obtain lubricant pressure distribution. However, in the Reynolds equation, they omitted the time derivative of pressure, and therefore, the resulting quasi-static solutions do not adequately represent the squeeze-film effects. In the same period, the literature presented many other models based on the "*p*-*h*" linearization of the time-dependent Reynolds equation and the equations of motion of a 4 DOF rigid rotor [10,11]. However, in most cases, these kinds of simplified models cannot replace more costly trial-and-error experimentation since they do not represent a reliable theory.

The results obtained from these models have become ever more accurate with the use of finite-difference and finite-element models, where the perturbation method can be used to compute stiffness and damping coefficients for rotordynamic models. Ref. [12] investigates the impulse response of an ultra-precision raster milling spindle under the simplified hypothesis of constant coefficients. In [13], critical speeds and modal shapes of a micro-spindle machine tools were evaluated, but this analysis is limited because it only considers bearings with constant stiffnesses and neglects their damping properties. A rotordynamic model of an air spindle with a diameter of 1/8" was considered [14], which involves coefficients that depend on the rotational speed. Here, the critical rotational speeds were calculated as a result of the rotordynamic analysis. A test rig to measure the rotordynamic response of a spindle supported on porous gas bearings is described in [15]. In [16–19], the dependence of stiffness and damping properties on both rotational speed and perturbation frequency is expressed with analytical formulations. With this technique, it is possible to calculate the critical speeds and verify the stability of rigid or flexible rotors. During the design process of a spindle supported on gas bearings, it is of great importance to estimate the conical and cylindrical natural frequencies. A clear design guideline is illustrated in [20] to avoid critical speeds within the operating speed range. The influence of temperature on the dynamic behavior of the spindle is also discussed. Conical and cylindrical modes are investigated in [21] for a non-uniform slot-clearance journal bearing, together with the onset speed of whirl instability.

The orbit method is an alternative to the linearized coefficient analysis employed in the aforementioned papers. It includes the complete nonlinear equations of the rotor and the bearings, which are integrated numerically to obtain the shaft center orbits corresponding to any set of geometrical, running, and initial conditions. This technique essentially uses the computer as an accurate experimental rig; it operates exactly in accordance with the assumed governing equations [22]. This method is currently used to calculate the unbalance response, the onset speed and the non-linear rotor behavior in case of large displacements [23]. Ref. [24] makes use of this method for a spindle with herringbonegrooved journal bearings. Ref. [25] highlights the differences between linear and non-linear analyses of a rotor supported by plain journal bearing. Past research by the authors concerns an electro-spindle of 50 mm shaft dia. supported on aerostatic bearings [26], a textile spindle [27] and a mesoscopic spindle of 10 mm diameter [28]. However, the experimental campaigns to validate the models are not always time- and cost-effective.

This paper proposes a cheap and simple experimental identification procedure for aerostatic spindles. Impulse and static tests are performed to identify the natural frequencies, damping ratio and static stiffness evaluated at the nose of an electro-spindle for PCB drilling. The natural frequencies and damping factors are measured by considering the impulse-response function of the system, whereas the stiffness at the nose spindle is detected through the application of static forces. These spindle features are evaluated at zero rotational speed and with different supply pressures. Moreover, the paper describes the non-linear numerical model of the spindle, which consists of a 4 DOF rigid and unsymmetrical rotor supported by two aerostatic bearings. The equations of motion of the rigid rotor are solved simultaneously with the isothermal time-dependent Reynolds equation through Euler's explicit method to obtain the lubricant pressure distribution. The main goal of the work is to validate the developed non-linear numerical model through the proposed identification procedure and the performed experimental tests.

#### **2. Experimental Setup**

#### *2.1. The Prototype*

Figure 1a,b show the investigated electro-spindle. During tests, the system was kept in horizontal position by means of a nylon block and a clamp. The spindle is driven by a permanent magnet (PM) synchronous motor that is placed between the rear and front aerostatic journal bearings. A sketch of the rotor-bearings system is shown in Figure 2. The bearings have slightly different diameters: 25 mm (front bearing) and 22 mm (rear bearing). Their axial length is 30 mm. They are supplied by means of two series of supply orifices of diameter *ds* = 0.15 mm, located at 7 mm from the borders of the bearings. Each series includes 12 holes equally spaced along the circumferential direction. Downstream each supply hole, a small pocket of dia. 1.1 mm and depth 0.2 mm is present. Considering the manufacturing tolerances, the radial air-film thicknesses are *hf* = 17.5 ± 1.5 μm and *hr* = 21.5 ± 1.5 μm on the front and rear bearings, respectively. Table 1 lists the mass properties of the rotor. The spindle is equipped with a water-cooling system to ensure temperature stability during operation.

**Table 1.** Mass properties of rotor.


**Figure 1.** (**a**) Photo of the electro-spindle in horizontal attitude; (**b**) photo of the test setup.

**Figure 2.** Sketch of the rotor-bearings system.

#### *2.2. Shaft-Displacement Detection*

Figure 2 shows the cartesian reference system, *Oxyz*, used to define the rotor center of mass and its position with respect to the bearings. The origin, *O*, is located at the front edge of the front journal bearing. Axis *z* is along the axial direction and directed from the front to the rear bearing. The axial coordinates of journal centers are indicated with *z <sup>f</sup>* and *zr* for front and rear bearings, respectively.

The rotor displacement is detected by means of two couples of capacitance probes located at *z*<sup>1</sup> = −21 mm (front plane) and *z*<sup>2</sup> = 91.5 mm (rear plane). The center of mass of the shaft is located at *zG* = 41 mm. The shaft center-of-mass displacements and the rotations about *x* and *y* axes are obtained from the sensor readings on front (*x*1, *y*1) and rear planes (*x*2, *y*2), as below:

$$\begin{array}{c} \chi\_{G} = \frac{x\_{1}(z\_{2} - z\_{G}) - x\_{2}(z\_{1} - z\_{G})}{z\_{2} - z\_{1}}\\ \mathcal{Y}\_{G} = \frac{y\_{1}(z\_{2} - z\_{G}) - y\_{2}(z\_{1} - z\_{G})}{z\_{2} - z\_{1}}\\ \mathcal{\theta}\_{\mathcal{Y}} = \frac{x\_{2} - x\_{1}}{z\_{2} - z\_{1}}\\ \mathcal{\theta}\_{\mathcal{X}} = \frac{y\_{1} - y\_{2}}{z\_{2} - z\_{1}} \end{array} \tag{1}$$

The shaft displacement at the spindle nose (*z* = *znose* = −40 mm) can be extrapolated using Equation (2):

$$\begin{array}{l} \mathcal{X}\_{\text{M08st}} = \frac{x\_1(z\_2 - z\_{n\text{nor}}) - x\_2(z\_1 - z\_{n\text{nor}})}{z\_2 - z\_1} \\ y\_{\text{nost}} = \frac{y\_1(z\_2 - z\_{n\text{nor}}) - y\_2(z\_1 - z\_{n\text{nor}})}{z\_2 - z\_1} \end{array} \tag{2}$$

#### *2.3. Static Stiffness on the Nose*

The static stiffness on the nose was evaluated by extrapolating the displacement produced through the application of a 30 ± 1 N load at *z* = *znose*. A dynamometer was employed to impose the load, and the spindle displacement at *z* = *znose* was computed by using Equation (2). Table 2 lists the resulting stiffness, measured at different supply pressures.

**Table 2.** Resulting static stiffness on the nose.


#### *2.4. Experimental Characterization at Null Rotational Speed*

The damped natural frequencies, *ωd*, and the damping factors, *ζ*, were experimentally investigated by means of impulse response tests. A vertical impulse was applied to the nose in correspondence of the tool artifact, and the spindle trajectory was detected by the capacitive sensors (capaNCDT CS05 by Microepsilon) (see Figure 1a). The tests were performed at null rotational speed and absolute supply pressures *ps* = 0.3, 0.5, 0.7 and 0.9 MPa. Each test was repeated four times, and signals were sampled with sampling frequency *fs* = 50 kHz.

The spectra of signals *y*<sup>1</sup> and *y*2, measured by sensors on the front and rear measuring planes, were evaluated through fast Fourier transform (FFT) in order to visualize the number of modes captured and their frequencies. Figure 3 shows one of the computed spectra (at 0.5 MPa supply pressure), while the related time signals are depicted in Figure 4. In all the investigated cases, only one rigid mode shape is visible. Meanwhile, the higher frequency can be attributed to the flexural mode of the spindle since it does not depend on the supply pressure of the air bearings. As can be seen from Figure 4, the first mode shape is conical since *y*<sup>1</sup> and *y*<sup>2</sup> present a 180◦ phase shift.

**Figure 3.** Example of spectra measured at 0.5 MPa supply pressure.

**Figure 4.** Example time signals at 0.5 MPa supply pressure.

#### 2.4.1. Single-DOF Modal-Identification Analysis through LDM and HPM

In view of these considerations, as a first approach, two single-DOF methods were employed to identify the modal parameters of the system: the logarithmic decrement method (LDM) [29] and the half-power method (HPM) [30]. Table 3 lists the mean values of the estimations of the damped frequencies, *ωd*, and the damping ratios, *ζ*, computed through HPM and LDM. The values of the damped natural frequencies increase with the supply pressure, whereas the damping ratios exhibit an opposite trend. The accuracy of these estimations was verified by comparing the experimental signals with those reconstructed through the identified modal parameters (see Figure 4) by making use of the following 1 DOF analytical formula:

$$\begin{array}{c} y\_{\rm rec}(t) \ = \ [a \cos(\omega\_d t) \ + b \sin(\omega\_d t)] \cdot e^{-\frac{\tau}{b}\omega\_n t} \\\ a = y\_{\rm rec}(0) \\\ b = \frac{\dot{y}\_{\rm rec}(0) + \xi \omega\_n y\_{\rm rec}(0)}{\omega\_d} \end{array} \tag{3}$$

where constants a and b depend on the initial conditions, and the undamped natural frequency is *ω<sup>n</sup>* = *ωd*/ <sup>+</sup><sup>1</sup> <sup>−</sup> *<sup>ζ</sup>*2.

**Table 3.** Experimental natural frequency and damping ratio of signal *y*<sup>2</sup> at different supply pressures; the results were obtained with HPM and LDM.


Unfortunately, as can be seen from Figure 4, the reconstructed signals do not exhibit a satisfactory degree of accuracy if the estimated parameters are not manually corrected with a trial-and-error procedure. This mismatch was mainly due to the fact that in some cases, the analyzed modes presented a relatively high damping factor (small number of oscillations in the time signal) and a flat peak in the spectrum. In fact, despite the large amount of energy provided by hammering the spindle nose, in most cases, the rotor stopped after a few oscillations. To overcome these problems, the estimations were repeated by means of the least-squares complex exponential method (LSCEM).

#### 2.4.2. Multi-DOF Modal Analysis with LSCEM

This is a multi-DOF method that works in the time domain [31]. It is based on the impulse-response function (IRF) of the MDOF system with N couples of complex and conjugate poles:

$$h\_{\vec{ij}}(t) \;= \sum\_{r=1}^{2N} A\_{r,\vec{ij}} \cdot e^{s\_r t} \tag{4}$$

where *sr* = *ζrωn*,*<sup>r</sup>* ± *i ωn*,*<sup>r</sup>* <sup>+</sup><sup>1</sup> <sup>−</sup> *<sup>ζ</sup>*<sup>2</sup> *<sup>r</sup>* and *Ar*,*ij* are the *rth* poles and the modal constants of the system and *i* and *j* refer to the nodes where the impulse is applied and where the system response is measured, respectively. The advantage of this method is that the nonlinear solution for the eigenvalues is computationally very straightforward since it performs an exponential fitting that requires no starting values. To obtain a solution, it is only necessary to supply the algorithm with the impulse-response function, *hij*(*t*), and the order of the model (*N*). A particular solution strategy was adopted in these cases. In order to provide the algorithm with the possibility to fit-noisy signals, the order of the model was chosen equal to 20, and frequency stabilization diagrams were used to distinguish the "computational modes" from the real ones. This kind of diagram was obtained by iterating the LSCEM from a minimum (3) to a maximum order (20) and superimposing the obtained frequencies on the obtained spectrum.

In order to obtain a cleaner stabilization diagram, before plotting the obtained damped frequencies, *ωd*, they were filtered by considering only the first 4 modes presenting the higher modal constant (*ω*1,*N*, *ω*2,*N*, *ω*3,*<sup>N</sup>* and *ω*4,*N*). Figures 5 and 6 show the spectrum and time signal reconstructed by using the modal parameters obtained from this procedure. In Figure 5, the results correspond to the "+" symbols, and they are coloured in blue, red, cyan and black from the higher to the lower modal constant. A further filter was imposed by considering a threshold of 5% on the values of the normalized modal constants in an attempt to exclude the computational modes. The results that satisfied this threshold are circled in green (see Figure 5). Figure 5 also shows the comparison between the experimental spectrum and that reconstructed through the LSCEM. Figure 6 shows the same comparison in the time domain.

**Figure 5.** The experimental and reconstructed spectra at 0.5 MPa supply pressure.

**Figure 6.** The experimental and reconstructed time signals at 0.5 MPa supply pressure.

As can be seen, in this case, the accuracy of the method is satisfactory. The results obtained with the LSCEM method are summed up in Table 4, where the damped natural frequencies and damping ratios are reported, together with their phase, *φ*, and modal constant, A. Only the first rigid mode was considered, as the second mode exhibited a natural frequency of about 220 krpm, which was almost constant with the supply pressure and had modal constants much smaller than the first mode. The evaluation of the phase shifts, Δ*φ* = |*φ*(*y*1) − *φ*(*y*2)|, confirms that the mode identified is of the tilting type.

**Table 4.** Experimental natural frequency and damping ratio of signals *y*<sup>1</sup> and *y*<sup>2</sup> at different supply pressures (only the first rigid mode is reported); the results were obtained with LSCEM.


#### 2.4.3. LSCEM vs. LDM and HPM (Experimental)

The damped frequencies obtained with the 1 DOF methods are in good accordance with those obtained with LSCEM; the error is less than 4%. The damping factors, apart from cases with *ps* = 0.3 MPa, differ by less than 36% for HPM and 61% for LDM.

#### **3. Numerical Model**

In this section, the numerical model of the spindle described in Section 2.1 is presented. The effect of the double-face thrust bearing is not taken into account, as it is negligible in the determination of radial and tilting stiffness with respect to the contribution of journal bearings.

The time-dependent Reynolds equation is considered in the fluid domain of journal bearings:

$$\frac{\partial}{\partial z} \left( \frac{ph^3}{12\mu RT} \frac{\partial p}{\partial z} \right) + \frac{\partial}{r\partial\theta} \left( \frac{ph^3}{12\mu RT} \frac{\partial p}{r\partial\theta} \right) + \mathcal{g}\_{\text{in}} = \frac{\omega r}{2RT} \frac{\partial (ph)}{r\partial\theta} + \frac{1}{RT} \frac{d(ph)}{dt} \tag{5}$$

where *gin* is the input flow per unit surface that crosses the supply orifices. The input flow is estimated using the analytical expression of the supply hole's discharge coefficient [32] without considering the Reynolds number dependence:

$$c\_d = 0.85 \left( 1 - e^{-8.2 \frac{h}{d\_s}} \right) \tag{6}$$

where *ds* is the supply hole's diameter and *h* is the local air gap under the orifice. The RE is discretized with finite difference technique, and the pressure distribution at time *t* + 1 is obtained from the previous pressure distribution with the Euler explicit method:

$$p\_{ij}^{t+1} = p\_{ij}^t + \Delta t \cdot f\left(p^t, \ p^{t+1}\right) \tag{7}$$

where *f* is a non-linear function that involves the pressure values in node *i*,*j* and in the adjacent nodes. The pocket downstream of each supply orifice is also considered, as it influences the dynamic characteristics of bearings. Once pressure *pij* is calculated at time *t* + 1, this value is immediately used to calculate *pi*+1,*<sup>j</sup>* at the same time. This improves the stability of the numerical scheme, as it facilitates the transmission of the information in the adjacent nodes. The film thickness depends on the rotor center-of-mass position (*xG*, *yG*) and on the rotor axis orientation (*ϑx*, *ϑy*):

$$h(z, \theta) := h\_0 - \left[x\_G + \vartheta\_y(z - z\_G)\right] \cos \theta - \left[y\_G - \vartheta\_x(z - z\_G)\right] \sin \theta \tag{8}$$

Once the pressure distribution is calculated, the shear stress on the rotor surface is given by:

$$
\pi = \frac{h}{2} \frac{\partial p}{r \partial \theta}\_{\text{el}} + \mu \frac{\omega r}{h} \tag{9}
$$

The 4 DOF equations of motion for a rigid rotor are considered:

$$\begin{array}{c} m\_r \ddot{\mathbf{x}}\_G - F\_{\mathbf{cx}} - F\_{\mathbf{ext}, \mathbf{x}} = m\_r e \omega^2 \cos \omega t\\ m\_r \ddot{\mathbf{y}}\_G - F\_{\mathbf{cy}} - F\_{\mathbf{ext}, \mathbf{y}} = m\_r e \omega^2 \sin \omega t\\ I \ddot{\theta}\_x + I\_p \omega \dot{\theta}\_y - M\_{\mathbf{cx}} - F\_{\mathbf{ext}, \mathbf{y}} (z\_G - z\_{n \text{osc}}) = (I\_p - I) \gamma \omega^2 \sin(\omega t + q\_1)\\ I \ddot{\theta}\_y - I\_p \omega \dot{\theta}\_x - M\_{\mathbf{cy}} + F\_{\mathbf{ext}, x} (z\_G - z\_{n \text{osc}}) = (I - I\_p) \omega^2 \gamma \cos(\omega t + q\_1) \end{array} \tag{10}$$

where *e* and *γ* are the static and the dynamic rotor unbalances, respectivelye; *ϕ*<sup>1</sup> is the angle that locates the dynamic unbalance plane with respect to the static unbalance plane; *Fext*,*<sup>x</sup>* and *Fext*,*<sup>y</sup>* are the external-force components acting on the rotor at coordinates *znose*; *Fcx*, *Fcy*, *Mcx* and *Mcy* are the reaction forces and moments of the bearings, calculated by integrating the pressure distribution and the shear stress in the fluid domain.

The moments are expressed with respect to the rotor center of mass:

$$
\begin{bmatrix} F\_{c\chi} \\ F\_{c\chi} \end{bmatrix} = \int\_0^L \int\_0^{2\pi} \left( -p \begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix} + \pi \begin{bmatrix} \sin\theta \\ -\cos\theta \end{bmatrix} \right) r d\theta dz \tag{11}
$$

The moments are expressed with respect to the rotor center of mass:

$$
\begin{bmatrix} M\_{\rm cx} \\ M\_{\rm cy} \end{bmatrix} = \int\_0^L \int\_0^{2\pi} \begin{bmatrix} (z - z\_G)(p\sin\theta + \tau\cos\theta) \\ (z - z\_G)(-p\cos\theta + \tau\sin\theta) \end{bmatrix} r d\theta dz \tag{12}
$$

The time integration of these equations of motion are carried out with the first order Euler method:

$$q\_{t+1} = q\_t + \dot{q}\_t \Delta t \tag{13}$$

)

where *q* is the state-space vector defined by

$$q = \begin{bmatrix} & \mathbf{x}\_G & \mathbf{y}\_G & \boldsymbol{\theta}\_x & \boldsymbol{\theta}\_y & \cdots & \dot{\mathbf{x}}\_G & \dot{\mathbf{y}}\_G & \boldsymbol{\theta}\_x & \dot{\boldsymbol{\theta}}\_y \end{bmatrix}$$

The algorithm performed for the coupled integration of the ODE and the PDE system is described below. After assuming an initial pressure distribution, for each time iteration, the following steps are followed:


*I*

This method is known as the "orbit method" and has the advantage of also taking into account the non-linearities of bearings [22]. Conversely, it is more time-consuming compared to linearized methods.

For static problems, the rotor position is fixed, and only the pressure distribution is considered time-dependent. In this case, only steps 1 and 5 are considered in order to find the steady-state solution. The iterative process is stopped when the load-capacity relative variation decreases below a given tolerance:

$$\left|\frac{F\_{\mathcal{c}}^{t+1} - F\_{\mathcal{c}}^{t}}{F\_{\mathcal{c}}^{t}}\right| < 10^{-6} \tag{14}$$

#### **4. Numerical Results**

The numerical model developed was validated using both literature benchmarks and numerical results. In particular, the numerical load capacity, *W*, and attitude angle, Φ, were compared with the analytical values for plain dynamic journal bearings. In the comparison, the static stiffness of the spindle measured in correspondence of the nose was also considered, together with the damped natural frequencies and the damping ratio at null rotational speed.

#### *4.1. Static Validation with Analytical Solution of Plain Dynamic Journal Bearings*

The mathematical model was validated in static conditions, comparing it with the analytical solution that exists for the plain dynamic journal bearing [7]. The dimensionless load capacity is expressed with a real and imaginary parts, which represent the in-phase and the out-of-phase components with respect to the line of centers:

$$\overline{W} = \frac{W}{\varepsilon p\_a LD} = \frac{i\Lambda}{1 + i\Lambda} \frac{\pi}{2} \left[ 1 - \frac{\tanh\left(\sqrt{1 + i\Lambda} \frac{L}{D}\right)}{\sqrt{1 + i\Lambda} \frac{L}{D}} \right] \tag{15}$$

where

$$
\Lambda = \frac{6\mu\omega}{p\_a} \left(\frac{r}{\mathbb{C}}\right)^2 \tag{16}
$$

is the bearing number, and

$$
\epsilon = \frac{\epsilon}{\mathbb{C}} \tag{17}
$$

the eccentricity ratio. The front and rear journal bearings were simulated with no supply orifices for comparison with the dynamic-bearing analytical solution. Figure 7 shows a perfect match, both regarding the load capacity and the attitude angle, Φ.

#### *4.2. Choice of Grid Resolution*

The influence of the grid-point number on the numerical solution is evaluated both in static and dynamic conditions. A total of 31 nodes were chosen along the axial direction for each bushing. Along the circumferential direction, 48 or 96 nodes were considered. Figures 8–10 compare the pressure distribution, load capacity and air consumption obtained with these computational grids in the presence of an eccentricity of *e* = 1 μm. As can be seen, the difference is minimal. As a result of these considerations, a 31 × 48 grid was used.

**Figure 8.** Axial (**a**) and circumferential (**b**) pressure distributions obtained with different numbers of grid points along the circumferential direction (48 and 96); *ps* = 0.5 MPa.

**Figure 9.** Load capacity vs. number of iterations obtained with different grid numbers along the circumferential direction.

**Figure 10.** Air flow vs. number of iterations obtained with different grid numbers along the circumferential direction.

#### *4.3. Static Stiffness of Journal Bearings*

The stiffness of journal bearings at zero speed was obtained from static simulations. A given displacement was imposed on the shaft (1 μm on radius), and the reaction forces in the bearings were computed in order to estimate the radial stiffness, *k <sup>f</sup>* and *kr*, of the front and rear bearings, respectively. The stiffness values are listed in Table 5, considering minimal and maximal radial clearances and taking into account the manufacturing tolerances, as shown in Section 2.1. Of course, stiffness values increase with supply pressure and decrease with increasing air gaps. The front bearing shows an increment of stiffness with the air gap only with *ps* = 0.3 MPa. The following explanation can be given after investigating the pressure distribution: the reason is that the front bearing is near the saturation condition; that is to say that the pressure drop through the input orifices is quite small (less than 0.03 MPa).


**Table 5.** Radial stiffness of journal bearings at different supply pressures and radial clearances.

The tilting stiffness of each journal bearing, calculated with respect to its center, is also evaluated in order to investigate whether its contribution to the overall tilting stiffness of the rotor bearings system is negligible or not. Table 6 compares the contribution of each bearing to the tilting stiffness obtained from the radial stiffness (left column) with the tilting stiffness of each bearing with respect to its center (right column). From this comparison, it is evident that in this case, this last bearing can be neglected so the overall tilting stiffness can be estimated from the radial stiffness of each journal bearing.


**Table 6.** Comparison between different contributions to tilting stiffness; results refer to geometrical configuration with *hf* = 19 μm, *hr* = 23 μm and supply pressure *ps* = 0.7 MPa.

#### *4.4. Stiffness on the Nose*

The overall stiffness, *knose*, measured on the nose at *z* = *zF* = −40 mm, is a function of the bearing stiffness according to:

$$k\_{\text{mse}} = k\_f \frac{z\_r - z\_f}{z\_r - z\_F} \tag{18}$$

Table 7 lists the radial stiffness on the nose extrapolated from the journal bearings stiffness calculated in the centered position, with 1 μm rotor eccentricity. In case an external radial load of 30 N is applied on the nose, the eccentricities on bearings are greater and the resulting nose stiffness is lower. Table 8 lists these values.

**Table 7.** Radial stiffness on the nose at different supply pressures and radial clearances, calculated at small eccentricity (1 μm).


**Table 8.** Radial stiffness of the nose at different supply pressures and radial clearances; these values were obtained with a 30 N load applied on the nose.


#### *4.5. Simplified Natural Frequencies*

At null rotational speed, in case the damping effects are neglected, the cylindrical natural frequency is expressed by:

$$
\omega\_{cyl} = \sqrt{\frac{k\_f + k\_r}{m\_r}}\tag{19}
$$

The conical natural frequency is given by:

$$
\omega\_{\rm corr} = \sqrt{\frac{\left(k\_f l\_f^2 + k\_r l\_r^2\right)}{I}} \tag{20}
$$

where *lf* and *lr* are the distances between the rotor center of mass and the centers of the front and the rear journal bearings, respectively. In this formulation, the eventual tilting stiffness of each journal bearing is neglected. Table 9 reports the undamped natural frequencies of the rigid rotor-bearings system at null rotational speed.


**Table 9.** Undamped natural frequencies at different supply pressures and radial clearances.

#### *4.6. Choice of the Time Step*

The optimal time step, *dt*, must be chosen for simulations in order to guarantee a stable solution and independence of the dynamic solution on the time step. The entity of the damping factor was found to be dependent on the time step, *dt*, and on the grid refinement. Figure 11 compares the values obtained for the damping factor at different rotational speeds and time steps, *dt*, with a 31 × 96 grid. As shown, the convergence is reached with a time step of *dt* = 10−<sup>8</sup> s, and a further decrease does not significantly change the damping factor. The effect of the grid on the damping factor is visible in Figure 12. When doubling the grid nodes along the circumferential direction, it is sufficient to halve the time step, *dt*, in order to obtain a very similar solution. As a result of these considerations, a 31 × 48 grid with *dt* = 2 × <sup>10</sup>−<sup>8</sup> s was chosen for dynamic simulations, as it was found to be a good compromise between accuracy and computational time.

**Figure 11.** Damping ratio vs. rotational speed obtained with different time steps; the grid is 31 × 96.

**Figure 12.** *Cont*.

**Figure 12.** Comparison of vibrations obtained on front and rear bearings with 31 × 96 and 31 × 48 grids; the two curves match fairly well in case the time step is halved when doubling the grid points.

#### *4.7. Damped Natural Frequencies and Damping-Factor Identification*

The damped natural frequencies and the damping factors were evaluated for the displacements along the x and y directions by the rotor in correspondence of planes at *z* = 0 (at the right edge of the front journal bearing) and at *z* = *zL* = 126 mm (at the left edge of the rear journal bearing). Moreover, to better identify both the cylindrical and the conical modes, the center-of-mass displacement, *yG*, and the tilting angles, *ϑ<sup>x</sup>* and *ϑy*, were considered. Two different initial conditions were simulated since the frequency content of the spectra of the system response may depend on how the system is excited. Condition 1 is defined by:

$$\begin{array}{l} \chi(z=0, t=0) = 0 \ \text{ $y(z=0, t=0)$ } = 0\\ \chi(z=z\_L, t=0) = 0, \ \text{ $y(z=z\_L, t=0)$ } = 0\\ \dot{\chi}(z=0, t=0) = 0 \ \text{ $\dot{y}(z=0, t=0)$ } = -0.01 \ \text{m/s} \\ \dot{\chi}(z=z\_L, t=0) = 0 \ \text{ $\dot{y}(z=z\_L, t=0)$ } = 0 \end{array}$$

Condition 1 was adopted in an attempt to reproduce the performed experimental tests. Vertical impulses were applied to the spindle nose, and the rotor oscillations were measured on the rear and front plane. Since the results obtained with the initial conditions (condition 1) were mainly governed by the conical mode, a second condition (condition 2) was simulated in an attempt to better figure out the cylindrical mode-shape of the rotor.

$$\begin{array}{ll} \mathbf{x}(z=z\_{G},t=0) & = 1 \text{ } \upmu\text{m}, \quad \mathbf{y}(z=z\_{G},t=0) & = 0\\ \boldsymbol{\theta}\_{\mathbf{x}}(t=0) & = 0, \quad \boldsymbol{\theta}\_{\mathbf{y}}(t=0) & = 0\\ \dot{\mathbf{x}}(z=z\_{G},t=0) & = 0, \dot{\mathbf{y}}(z=z\_{G},t=0) & = -0.01 \text{ m/s} \\ \dot{\boldsymbol{\theta}}\_{\mathbf{x}}(t=0) & = 0, \dot{\boldsymbol{\theta}}\_{\mathbf{y}}(t=0) & = 0 \end{array}$$

Condition 2 was chosen to simultaneously produce a whirl and a translation of the rotor. Due to the non-symmetry of the rotor-bearings system, both modes were visible both with initial condition 1 and with condition 2. The frequency content of rotor displacements puts in evidence the general presence of two natural frequencies, as seen in Figure 13. The prominence of one with respect to the other depends on the initial condition.

**Figure 13.** Examples of frequency spectra of signals *yG*, *y*0, *θ<sup>x</sup>* and *yL* obtained from the orbit method.

4.7.1. Logarithmic Decrement Method (LDM)

Although the LSCEM proved to be the most reliable method, for the sake of completeness, the LDM method was also used to identify the main damping factors (with the largest modal constant) and natural frequencies of the numerical signals *y*0, *yL*, *yG* and *θx*. The LDM was used, considering time intervals Δ*T* of 5 and 10 ms. The resulting values of the damping factors were compared to evaluate the sensitivity of the method with respect to the signal length. The main frequency of each signal was also evaluated by measuring the average period of the oscillations present in the considered time intervals. This frequency can be compared with the peaks visible in the related spectrum. Tables A1–A4 in Appendix A list both the natural frequencies resulting from the FFT analysis and the frequencies and the damping factor of the signal resulting from the logarithmic decrement method. When the oscillation was highly damped (*ps* = 0.3 and 0.5 MPa), only a 5 ms time window was considered. The following considerations can be made based on the data analysis:


#### 4.7.2. LSCE Method

The frequency content of the numerical signals *yL* and *y*<sup>0</sup> was investigated by adopting a procedure similar to that adopted in Section 2.4. To obtain more reliable frequencystabilization diagrams, the adopted procedure was slightly different than that used for the experimental data. In this case, to favor the exponential fitting, the model order used in the LSCEM was 80, and a random noisy signal was added to the numerical signals [31]. Figures 14 and 15 show the comparison obtained between the original and reconstructed data. The damped frequencies and damping factors evaluated on signals *y*<sup>0</sup> and *yL* are listed in Table 10, together with their phase shift and modal constants.

**Figure 14.** Stabilization diagram.

**Figure 15.** Comparison in time domain of the reconstructed signal with the numerical signal.


**Table 10.** Modal parameters resulting from LSCEM.

#### 4.7.3. Comparison of LSCEM and LDM for Numerical Simulations

When LDM and LSCEM are compared for the numerical simulations, the damped frequencies obtained with LSCEM differ by less than 1.6% with respect the frequencies obtained with LDM, apart from the case with *ps* = 0.3 MPa. The estimations of damping factors differ by much more—up to 80%. This parameter is more sensitive to variations in the air gap and pressure.

#### **5. Comparison between Experimental and Numerical Results**

The first critical speed of the spindle was experimentally measured by analyzing the amplitude of the synchronous whirl at different rotational speeds. It was found to be 72 krpm and 94 krpm at gauge supply pressure *ps* = 0.5 and 0.7 MPa, respectively. If compared with the theoretical frequencies resulting from the simplified model of Section 4.5, the difference is less than 6%.

The experimental stiffness on the nose at different pressures is lower than the theoretical value calculated near the rotor-centered condition (between 27% and 39% with respect to the theoretical value). Considering that the 30 N external load imposed on the nose involves non-negligible eccentricity ratios on bearings (especially on the front bearing and with low pressure), a non-linear estimation of the stiffness is preferred (see Section 4.4). If compared with the stiffness obtained in this case, the experimental stiffness is more in accordance with the prediction (between 52% and 86% with respect to the theoretical value).

The experimental damped frequencies and the damping factors evaluated on the shaft vertical vibration at *z* = *zL* are compared with the numerical frequencies in Figure 16. The experimental damped natural frequencies are in good accordance with the numerical frequencies, as they fall in the range defined on the basis of the tolerance considered on the front (*hf* = 17.5 ± 1.5 μm) and rear (*hr* = 21.5 ± 1.5 μm) clearances. Conversely, the computed experimental damping factors are higher than the numerical simulations. The single-DOF methods were found to be less accurate than LSCEM to capture the modal parameters of experimental signals, especially the damping factors.

**Figure 16.** Comparison of experimental and numerical damped frequencies (**a**) and damping factors (**b**) on rear plane.

#### **6. Conclusions**

This paper describes the experimental and numerical investigation of an electrospindle supported by aerostatic bearings. The numerical model was validated against literature results in the case of plain aerodynamic bearings and against static and dynamic experimental tests carried out on the electro-spindle at null rotational speed. Both single-DOF methods and a multi-DOF method were proposed to estimate the damped natural frequencies and the damping factors. The accuracy of the model results with respect to the experimental data was discussed. In particular, the damping factor evaluated with LSCEM was in good agreement with the experimental data, while the single-DOF methods overestimated it. Future investigations will regard the dynamic behavior of the spindle at different rotational speeds in order to study its unbalance response and stability.

**Author Contributions:** Conceptualization, F.C., L.L., A.T., T.R. and V.V.; methodology, F.C. and L.L.; software, F.C. and L.L.; writing—review and editing, F.C., L.L. and A.T.; supervision, V.V.; project administration, T.R.; funding acquisition, T.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Carbomech Company.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

#### **Nomenclature**


#### **Appendix A**

Tables A1–A4 report both the natural frequencies resulting from the FFT analysis, both the frequencies and the damping factor of the signal resulting from the logarithmic decrement method.

**Table A1.** Damped frequencies and damping factors resulting from FFT analysis and from LDM; *ps* = 0.3 MPa.



**Table A2.** Damped frequencies and damping factors resulting from FFT analysis and from LDM; *ps* = 0.5 MPa.

**Table A3.** Damped frequencies and damping factors resulting from FFT analysis and from LDM; *ps* = 0.7 MPa.



**Table A3.** *Cont.*

**Table A4.** Damped frequencies and damping factors resulting from FFT analysis and from LDM; *ps* = 0.9 MPa.



**Table A4.** *Cont.*

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34 www.mdpi.com

ISBN 978-3-0365-5608-6