Mathematical Modeling on Statics and Dynamics of Aerostatic Thrust Bearing with External Combined Throttling and Elastic Orifice Fluid Flow Regulation

As aerostatic bearings are used in high-speed metal-cutting machines to increase machining accuracy, there is the need to improve their characteristics, including compliance, which is usually high. In practical applications, a significant reduction of bearing compliance is often necessary, sometimes down to zero and even negative values, to ensure automatic compensation of the elastic deformation in the machine technological system. A decrease in compliance leads to deterioration in the dynamic performance of the bearing, so it is necessary to develop new designs that meet the above requirements. This article considers an aerostatic bearing, in which decrease in compliance is ensured by the use of air throttling with elastic orifices. To ensure its stability, the principle of combined external throttling was applied, which can substantially improve the dynamics of conventional aerostatic bearings. A mathematical model of the elastic orifice deformation was developed, together with the flow rate performance calculation method. The method ensured full qualitative and satisfactory quantitative agreement with the experimental data. The model was used in the mathematical modeling of the aerostatic bearing movement. The article also proposes a method to calculate the static load capacity and compliance of a bearing, as well as a numerical method for fast computation of its dynamic performance, which allows for real-time multi-parameter optimization by the bearing dynamic performance criteria. The study showed that there is an optimal set of design parameters for which low, zero, and negative static compliance of the bearing is ensured, with the necessary stability margin, high speed, and the non-oscillatory nature of the transient processes.


Introduction
The exceptional advantages of aerostatic bearings are the high rotation speed and extremely low friction [1][2][3][4][5]. The bearings are lubricated with compressed gas from the source into the flow path, which usually uses passive elements to limit the gas flow to create a certain pressure which is necessary to maintain the current load and to control the displacement of the moving part of the structure, thereby ensuring its stiffness [4][5][6][7][8][9][10][11][12][13][14][15][16]. And although the development of technology has contributed to the improvement of bearings, they still have a drawback-low stiffness [8,10,17,18].
In 1960, John H. Laub and Frank Batsch of the California Institute of Technology's Jet Propulsion Laboratory announced experiments with elastic orifices as a simple means to control the flow of gas lubricant [19], that helped to improve the performance of structures by increasing stiffness, experimental studies at the Stanford Research Institute demonstrated the dependence of an elastic orifice performance on the geometry and material properties [20,21].
However, the characteristics of such bearings and the optimal operating conditions have not been studied and as well the methods for their calculation have not been developed. Attempts at mathematical modeling and calculation of elastic orifices did not lead Newgard P. M. and Kiang R. L. to results that would give satisfactory agreement with experimental data [22].
The thrust bearing design described by these researchers, apparently, is not very suitable for practical use, because with a small chamber depth the elastic orifice cannot be deformed, while with a large depth the dynamic performance of the bearing will inevitably lead to instability.
Bearings with external combined throttling of gas have certain prospects, which make it possible to ensure stable operation of structures even in the presence of a relatively large volume of gas enclosed in an intermediate chamber [14,23,24].
The literature indicates that replacing a passive throttling diaphragm in the prototype with an elastic orifice produces a bearing which may have a potential to reduce stiffness. The mathematical modeling and theoretical study of the static and dynamic characteristics of such a bearing is actually the goal of this study.
The article proposes a method for calculating the static bearing capacity and bearing compliance, as well as a numerical method for quickly calculating its dynamic characteristics, which allows performing multi-parameter optimization of the criteria for the quality of bearing dynamics in real time. The study showed that there is an optimal set of design parameters for which a low, zero, and negative static bearing is ensured, with the necessary stability margin, high speed and non-oscillatory nature of the transients.

Design Scheme of the Bearing
The layout of a circular type bearing with a throttling elastic orifice and damping annular diaphragms is shown in Figure 1. The design consists of a base 1 and a movable part 2, separated by a lubricant gap of thickness h. Gas under constant pressure ps through a throttling diaphragm of radius rp of the elastic orifice 3 is supplied to the chamber 4 of volume vp, where the pressure becomes equal to pp < ps. Then, through damping annular diaphragms of diameter dk evenly spaced along a circle of radius r1, gas under pressure pk < pp enters the bearing layer of the bearing, having overcome it and then flows into the environment at pressure pa.
The bearing operates as follows. A change in the external force f on the bearing first causes a change in pressure pk, then a pressure pp, as a result of which the elastic orifice is deformed and the radius rp of the throttling diaphragm changes, as shown in Figure 2. As a result, an additional change in gas flow through the elastic orifice occurs. By a targeted choice of the elasticity of the elastic orifice, The design consists of a base 1 and a movable part 2, separated by a lubricant gap of thickness h. Gas under constant pressure p s through a throttling diaphragm of radius r p of the elastic orifice 3 is supplied to the chamber 4 of volume v p , where the pressure becomes equal to p p < p s . Then, through damping annular diaphragms of diameter d k evenly spaced along a circle of radius r 1 , gas under pressure p k < p p enters the bearing layer of the bearing, having overcome it and then flows into the environment at pressure p a .
The bearing operates as follows. A change in the external force f on the bearing first causes a change in pressure p k , then a pressure p p , as a result of which the elastic orifice is deformed and the radius r p of the throttling diaphragm changes, as shown in Figure 2. As a result, an additional change in gas flow through the elastic orifice occurs. By a targeted choice of the elasticity of the elastic orifice, it is possible to obtain the desired shape of the load characteristic and bearing stiffness in the calculation area. it is possible to obtain the desired shape of the load characteristic and bearing stiffness in the calculation area. In a theoretical analysis of the bearing under study, it is convenient to replace the stiffness function with its inverse compliance function, which is the transfer function of the elastic system and is determined by the limit of the ratio of the output quantity (change in the thickness h of the carrier gap) to the input value (load change f).

Mathematical Modeling of Elastic Orifice Deformation
The design scheme of the elastic orifice deformation is shown in Figure 2. The modeling of elastic orifice deformation is based on the theory of large deflections of thick plates. Below, is a developed mathematical model.

Elastic Orifice Deformation Model
As can be seen from Figure 2, the radius of the orifice can be calculated by the formula 0 0 ψ δ sin arctg where rp0 is the radius of the undeformed orifice, r is the current radius, ψ = ψ(r) is the deflection function, u = u (r) is the tensile function of the middle surface of the elastic orifice. Considering the elastic orifice as a round plate with large deflections, we obtain the resolved system for the function in the form [25] ( where φ = φ(r) is the Airy stress function, rs is the outer radius of the elastic orifice, E is the elastic modulus of the orifice material, ν is Poisson's ratio [26]. During operation, the outer edge of the elastic orifice does not experience deflection, therefore In a theoretical analysis of the bearing under study, it is convenient to replace the stiffness function with its inverse compliance function, which is the transfer function of the elastic system and is determined by the limit of the ratio of the output quantity (change in the thickness h of the carrier gap) to the input value (load change f ).

Mathematical Modeling of Elastic Orifice Deformation
The design scheme of the elastic orifice deformation is shown in Figure 2. The modeling of elastic orifice deformation is based on the theory of large deflections of thick plates. Below, is a developed mathematical model.

Elastic Orifice Deformation Model
As can be seen from Figure 2, the radius of the orifice can be calculated by the formula where r p0 is the radius of the undeformed orifice, r is the current radius, ψ = ψ(r) is the deflection function, u = u (r) is the tensile function of the middle surface of the elastic orifice. Considering the elastic orifice as a round plate with large deflections, we obtain the resolved system for the function in the form [25] where ϕ = ϕ(r) is the Airy stress function, r s is the outer radius of the elastic orifice, E is the elastic modulus of the orifice material, ν is Poisson's ratio [26].
During operation, the outer edge of the elastic orifice does not experience deflection, therefore ψ(r s ) = 0 (4) On the inner and outer radii of the orifice there are no bending moments and normal stresses in the radial direction, therefore Within the framework of Hooke's law, the relation between the tensile function and the Airy function has the form [25] u r The system of Equations (2), (3), (9), taking into account the boundary conditions (4)- (8), is closed and allows one to determine the functions ψ and u. Using these functions, one can find the radius r p of the deformed elastic orifice by the formula (1). The proposed calculation model (1)-(9) of elastic orifice deformation is a first approximation. In a more rigorous formulation, it is necessary to take into account the change in the tensile function not only in the radius but also in the height of the orifice.
The calculation of the characteristics of the elastic orifice are carried out in a dimensionless form using dimensionless quantities Differential Equations (2), (3) are nonlinear, therefore, to solve the problem (1)-(9), an iterative process was constructed, at each step of which a system of linearized differential equations was considered d dR with boundary conditions where n = 0,1, 2, . . . is the iteration number, The initial approximation for the function T(R) is its value for an undeformed elastic orifice. The linearized boundary value problems (10), (14), (15) and (11)- (13) were solved at each iteration using the finite-difference sweep method [27]. The iterative process (10)- (15) was stopped when the condition max T (n+1) − T (n) < eps (16) where eps is some small positive number (in the calculations eps = 10 −4 is accepted). The dimensionless radius of the elastic orifice diaphragm, taking into account (1), (9), (15), was calculated by the formula Based on the presented formulas, the radius of the orifice was calculated under the action of the differential pressure in the bearing on the elastic orifice.

Analysis of Elastic Orifice Calculation Results
The results of the described method for calculating the elastic orifice were compared with the experimental data of Newgard P. M. and Kiang R. L. [22] for two elastic orifices: 1-for r p0 = 0.125 mm, δ = 2.64 mm; 2-for r p0 = 0.19 mm, δ = 1.14 mm. The comparison results are shown in Figure 3.  where cp = 0.85 the empirical coefficient of an orifice diaphragm [17], R0 = 287 m 2 /(c 2 K) is the gas constant, T0 = 293 K is the absolute temperature, Γ = 4.313 is the adiabatic constant of air [17]. It can be seen from Figure 3 that these curves have a complete qualitative and satisfactory The graph shows the theoretical (solid lines) and experimental (dashed lines) sensitivity curves of various-sized elastic orifices with Shore hardness A − 45 (E = 2.16 MN/m 2 , ν = 0.5) at p s = 238 kPa. The pressure difference p s − p p is plotted along the x-axis, and the mass flow q of gas through the orifice is plotted along the y-axis. The calculations used the formula where c p = 0.85 the empirical coefficient of an orifice diaphragm [17], R 0 = 287 m 2 /(c 2 K) is the gas constant, T 0 = 293 K is the absolute temperature, Γ = 4.313 is the adiabatic constant of air [17]. It can be seen from Figure 3 that these curves have a complete qualitative and satisfactory quantitative agreement. Therefore, the described method can be used to a first approximation in the calculation of supports with an elastic orifice.
The discrepancy between the theoretical and experimental curves at large pressure drops (at large deformations of the elastic orifice) can be explained by the fact that with an increase in the strain, the error increases in the approximate formula (1), in Equations (2), (3), which do not take into account the changes in the functions ψ and u in the height of the elastic orifice, as well as in relation (9) based on Hooke's law, which for rubber at large deformations gives a progressive error. At pressure drop values from low to moderate the discrepancy is no more than 15%-20%, which is quite acceptable. High values of the discrepancy appear only at excessively large differences, however, these modes are not of interest, because in this case the acceptable load of the bearing is way too low.
When calculating the compensating effect of an elastic orifice, the dimensionless input parameters varied in the region P s ∈ [1.5, 6], Figures 4 and 5 shows the calculation results for the elastic orifice, showing that R p substantially depends on ∆ and ε and only slightly depends on ν.
The working radius R p is an almost linear function of the differential pressure, therefore, it can be approximated with the expression where K e is the elastic ratio of the orifice. (18) approximately determine the values of K e . For the curves corresponding to ε = 15, we have approximate Equation (18), which implies that K e ≈ 0.95. For the curves corresponding to ε = 30, a similar equation has the form 0.0119 ≈ 0.02(1 − 1.2K e ), whence it follows that K e ≈ 0.34.

Bearing Mathematical Model
To model the bearings presented on Figure 1, we used a mathematical model of motion of a gas-static bearing with combined external throttling [23], which based on Reynolds differential equation [3,17], describes the pressure distribution in a thin lubricating gap, and uses the position of "continuous pressurization line" to replace discrete annular diaphragms with a continuous pressure function on the line of their locations [28].
The model takes into account the difference regarding the input throttling resistance, which is an elastic orifice. The dimensionless flow rate through this resistance is expressed by the formula where the change in the radius of the passage of the elastic orifice is determined by formula (18).

Bearing Mathematical Model
To model the bearings presented on Figure 1, we used a mathematical model of motion of a gas-static bearing with combined external throttling [23], which based on Reynolds differential equation [3,17], describes the pressure distribution in a thin lubricating gap, and uses the position of "continuous pressurization line" to replace discrete annular diaphragms with a continuous pressure function on the line of their locations [28].
The model takes into account the difference regarding the input throttling resistance, which is an elastic orifice. The dimensionless flow rate through this resistance is expressed by the formula where the change in the radius of the passage of the elastic orifice is determined by formula (18).

Bearing Static Model
The static model will include two equations of equality of the flow rate Q p through the elastic orifice and Q k through the damping annular diaphragms, the flow rates Q k and Q h at the entrance to the bearing layer of the dimensionless thickness H, as well as the equation of force equilibrium of the movable element 2, connecting the bearing capacity W and the load F to the movable element where A h = −1/ ln R c , Π is the Prandtl expiration function [23]. Obviously, the method of calculating the load characteristic H(F) used for conventional bearings, in which the input value is convenient to use H and the output F = W, is not suitable in this case, since, as will be shown below, this function can be multi-valued, i.e., several values of H may correspond to one value of F. Moreover, for certain values of this variable, the value of F may not be found.
Therefore, a parametric method was used to determine the load F and the gap H, as a function of pressure P p F = F(P p ), H = H(P p ).
When calculating the static characteristics, we used the boost pressure P s , radius R c , the tuning factors of the combined external throttling system χ, and the elastic ratio of the elastic orifice K e .
The method for calculating static characteristics includes determination of the "calculated point" parameters corresponding to the dimensionless thickness of the bearing layer H = 1, -calculation of flow rate Q(F) and load H(F) characteristics of the bearing.
In the "calculated point" mode, the coefficients A k , A p of expenses (21), (22) are determined. First, we calculate the coefficient A h and determine the pressure P k , P p , corresponding to the value of the calculated clearance H = 1 then we calculate the coefficients When calculating the flow and load characteristics, we sequentially set n p pressure values P p ∈ [1, P s ] according to the formula P p = 1 + i p (P s − 1)/n p , i p = 0, 1, . . . , n p . For each of its values, we can find the flow Q p according to formula (19).
The first two Equations (20) allow us to write a nonlinear equation for an unknown pressure P k ∈ 1, P p in the form where The solution of Equation (25) was found by the bisection method [27]. The integral in formula (23) of the bearing capacity W was calculated by the Simpson rule [18].
Bearing compliance is determined by the formula Equation (20) make it possible to establish the values of the elastic ratio K e , at which the bearing reaches zero compliance (K = 0).
To determine the compliance of the bearing, we give the load F a small perturbation ∆F. The response will be disturbances ∆H, ∆P k , ∆P p , quantities H, P k , P p , respectively. We write the system of Equation (20) with respect to disturbances in the matrix form where Solving (28) by the Cramer rule [27], we obtain the formula for bearing compliance Using (29), it is possible to set the value of the elastic ratio K e , at which the bearing has zero compliance (K = 0). Most simply, this can be done for the "calculated point" mode H = 1. Using (19), (20), (22), (29), we find In calculations it is sometimes convenient to use an inverse relationship It has been found that for relatively large K e , the bearing can lose stiffness. This is the case when the denominator (29) vanishes. It is easy to verify that this occurs when Note that the critical value of the elastic ratio K e calculated by formula (32) occurs when R p > 0, that is, when the orifice is still able to pass gas during deformation. It closes only when the pressure drop is too large and K e ≥ 1 P s −P p . Figure 6 shows the load characteristics H(F) of an elastic orifice bearing for different values of elastic ratio K e . All curves intersect at one point corresponding to the "calculated point" mode H = 1.

Static Characteristics of the Bearing and Their Discussion
It has been found that for relatively large Ke, the bearing can lose stiffness. This is the case when the denominator (29) vanishes. It is easy to verify that this occurs when Note that the critical value of the elastic ratio Ke calculated by formula (32) occurs when Rp > 0, that is, when the orifice is still able to pass gas during deformation. It closes only when the pressure drop is too large and 1 .  It can be seen that with increasing Ke, the curvature of the lines changes significantly. When Ke = 0, that is for the rigid orifice 3, the compliance of the bearing K is always positive. With increasing Ke in some parts of the load characteristic, the compliance decreases to zero (K = 0), which corresponds to infinite stiffness, and even a negative value (K < 0). So, at Ke = 0.22 and F ≈ 1.5, the bearing has zero compliance. At Ke = 0.22, one can find the region 0.35 < F < 4.4, in which the bearing has negative static It can be seen that with increasing K e , the curvature of the lines changes significantly. When K e = 0, that is for the rigid orifice 3, the compliance of the bearing K is always positive. With increasing K e in some parts of the load characteristic, the compliance decreases to zero (K = 0), which corresponds to infinite stiffness, and even a negative value (K < 0). So, at K e = 0.22 and F ≈ 1.5, the bearing has zero compliance. At K e = 0.22, one can find the region 0.35 < F < 4.4, in which the bearing has negative static compliance (K < 0). Figure 7 shows graphs of expenditure flow rate dependencies Q(F). It is seen that with an increase in the coefficient K e , the nature of the curves changes sharply. Comparison of graphs Figures 6 and 7 shows that the less the compliance of the bearing, the lower the mass flow rate of air.   This is explained by the fact that at low loads the pressure difference Ps − Pp is large, and as follows from formulas (18), (19) for Ke > 0, the diaphragm radius Rp will be small, and, therefore, the flow rate Qp will be small. As the load F increases, the pressure difference Ps − Pp decreases and, consequently, the radius Rp and the flow rate Qp increase, thereby providing an additional supply of lubricant to the bearing, which contributes to a more intensive decrease in bearing compliance K.
An analysis of the curves shows that the dependences K(ς) for the modes of positive compliance are monotonous in nature, while for zero and negative compliance they are unimodal in nature, where the minimum of the dependencies falls on a certain value of ς. Moreover, the smaller this minimum, the greater the corresponding ς.     This is explained by the fact that at low loads the pressure difference Ps − Pp is large, and as follows from formulas (18), (19) for Ke > 0, the diaphragm radius Rp will be small, and, therefore, the flow rate Qp will be small. As the load F increases, the pressure difference Ps − Pp decreases and, consequently, the radius Rp and the flow rate Qp increase, thereby providing an additional supply of lubricant to the bearing, which contributes to a more intensive decrease in bearing compliance K.
An analysis of the curves shows that the dependences K(ς) for the modes of positive compliance are monotonous in nature, while for zero and negative compliance they are unimodal in nature, where the minimum of the dependencies falls on a certain value of ς. Moreover, the smaller this minimum, the greater the corresponding ς. This is explained by the fact that at low loads the pressure difference P s − P p is large, and as follows from formulas (18), (19) for K e > 0, the diaphragm radius R p will be small, and, therefore, the flow rate Q p will be small. As the load F increases, the pressure difference P s − P p decreases and, consequently, the radius R p and the flow rate Q p increase, thereby providing an additional supply of lubricant to the bearing, which contributes to a more intensive decrease in bearing compliance K.
An analysis of the curves shows that the dependences K(ς) for the modes of positive compliance are monotonous in nature, while for zero and negative compliance they are unimodal in nature, where the minimum of the dependencies falls on a certain value of ς. Moreover, the smaller this minimum, the greater the corresponding ς. Figure 9 shows the curves of the dependence of the elastic parameter K e0 on the coefficient ς of the tuning of the damping annular diaphragms for the "calculated point" mode H = 1, at which the bearing has zero compliance. For smaller χ and ς, smaller K e = K e0 are required. At the same time, when setting to too small χ, there may be no modes in which the bearing reaches zero compliance. Graphs of the dependence of the critical values of the elastic ratio Ked(ς), at which the bearing loses its static stiffness, and, consequently, stability, are shown in Figure 10. The Ked values determined by formula (31) give the upper boundary on the admissible values of the elastic ratio Ke = Ked. For Ke > Ked, the bearing is also statically unstable.

Bearing Dynamic Model
The mathematical model of the bearing dynamics is based on a bearing motion model with combined external throttling [23]. The only difference is that the formula for the deviation of the Laplace transform of the gas flow through the throttling simple diaphragm is replaced with the corresponding flow formula through the elastic orifice, which has the form Graphs of the dependence of the critical values of the elastic ratio K ed (ς), at which the bearing loses its static stiffness, and, consequently, stability, are shown in Figure 10. Graphs of the dependence of the critical values of the elastic ratio Ked(ς), at which the bearing loses its static stiffness, and, consequently, stability, are shown in Figure 10. The Ked values determined by formula (31) give the upper boundary on the admissible values of the elastic ratio Ke = Ked. For Ke > Ked, the bearing is also statically unstable.

Bearing Dynamic Model
The mathematical model of the bearing dynamics is based on a bearing motion model with combined external throttling [23]. The only difference is that the formula for the deviation of the Laplace transform of the gas flow through the throttling simple diaphragm is replaced with the corresponding flow formula through the elastic orifice, which has the form The K ed values determined by formula (31) give the upper boundary on the admissible values of the elastic ratio K e = K ed . For K e > K ed , the bearing is also statically unstable.

Bearing Dynamic Model
The mathematical model of the bearing dynamics is based on a bearing motion model with combined external throttling [23]. The only difference is that the formula for the deviation of the Laplace transform of the gas flow through the throttling simple diaphragm is replaced with the corresponding flow formula through the elastic orifice, which has the form where R p , P p are static parameters.

Method for Determining the Transfer Function Coefficients of Dynamic Bearing Compliance
In the calculations, a new method was used to determine the criteria for the quality of bearing dynamics, the essence of which is described below.
The transfer function (TF) of the dynamic compliance of the bearing, as a system of automatic control and regulation can be represented in the form . . + a n s n (33) where n > 0, m > 0, n > m, s is the Laplace transform variable [29]. The denominator (33) is the characteristic polynomial (CP) of a linear dynamic system with distributed parameters, the use of which in the study of its dynamics quality allows us to determine the stability margin and evaluate the system's speed from the roots of the characteristic equation [18,29].
For a system with distributed parameters, it is necessary to use a method that would ensure its presentation in the form (33) based on the calculation of the criteria of stability margin and speed of devices with a predetermined accuracy.
The representation of the TF in the form of (33) falls under the classical problem of rational interpolation [30], the solution of which, however, does not provide an exhaustive answer to the question of the accuracy of the system stability criteria obtained by root methods using the characteristic equation, because the value of the degree of CP is unknown in advance.
Existing methods of linear interpolation are based on solving a linear system of equations with respect to coefficients (33), which contains n + m equations [30]. Such systems can be solved using general methods, for example, the Gauss-Jordan method, which has a cubic order of complexity (n + m) 3 (hereinafter, the order of complexity of the computational method means the time complexity of the algorithm that implements it [30,31]). For large n and m, this can entail a significant expenditure of computer time in the process of multi-parameter optimization of a dynamic system. This article presents a quick method to find the coefficients (33). It is based on solving systems of linear equations of a special form of a substantially smaller order, which allows one to find their solution using fast methods with a quadratic order of complexity m(n + m), which contributes to a significant acceleration of the optimization procedure for dynamic systems. If it is required to find only CP coefficients, then the order of complexity of the method is n 2 .
When conducting rational interpolation, the degrees n and m of polynomials (33) should be known. However, their acceptable values can be obtained only on the basis of satisfactory accuracy in determining the quality criteria of the dynamics of the system.
If it is impossible to determine the values of these parameters without calculating the mentioned criteria, then finding their difference is not difficult.
Indeed, if a n 0 and b m 0, then the infinite limit where p = n − m. Usually, the difference p is calculated in several units, more often p = 1 or 2, therefore, you can find this limit and p can be fast enough.
We set the value of m, define n = m + p and find the coefficient We rewrite (1) in the form a 1 + a 2 s + . . . + a n s n−1 where Equation (35) contains k = n + m unknown coefficients.
Substituting s = s j , (j = 1, 2, . . . , k) in sequence (35), we obtain a system of linear equations for unknown coefficients (33) Mx = y where We represent the matrix of the system in the form where U is the matrix, Φ is the matrix of the discrete Fourier transform [32] Its inverse matrix is determined by the formula Multiplying Φ −1 by (37), we bring the system (35) to the form For m > 0, the matrix U has a cellular structure of the form where E and 0 are the identity and zero matrices of size n × n and m × n, C and D are Toeplitz matrices of size n × m and m × m, respectively [31]. Indeed, the blocks of matrices Φ −1 and M cells are mutually inverse matrices of the discrete Fourier transform, therefore, their product will give the identity matrix E. Elements of the block of cells are obtained by multiplying the rows of the matrix Φ −1 and the columns of the matrix M, which are also elements of the direct and inverse matrices Fourier transform. The sums of their products, giving off-diagonal elements of the identity matrix, will be zeros by analogy with the way this holds for the zero elements of the block E located above them.
The nature of the matrices A and B is explained by the fact that the elements of the columns of the matrix M for j > n are formed by sums of products of displaced elements of the matrix Φ and the elements of vector Γ that are different from unity. In such cases, their scalar products yield Toeplitz matrices [31].
From Equations (37) and (38) it follows that for m > 0 the vector b of the coefficients of the TF (33) satisfies the system of equations where d is a vector. Compared to the original Equation (35) of order k = n + m, Equation (39) has a significantly lower order of m < k/2. Therefore, the solution for it can be obtained much faster. Equation (39) is a standard problem with an asymmetric Toeplitz matrix of a special form [33], which can be solved both using general methods, for example, the Gauss-Jordan method [34], the complexity of which is proportional to m 3 , and using special fast methods taking into account the features of Equation (39) and having complexity proportional to m 2 . The latter include the methods described in [33,34].
With this in mind, we can dispense with the matrix C and, using the solution of system (10), quickly find the CP coefficients using the complexity nm formula With this in mind, the total complexity order of this method for finding the coefficients (33) is m(n + m).

Quality Criteria for Bearing Dynamics
To assess the quality of the dynamics of linear systems, root criteria are often used [29]: The degree of stability η characterizes the speed of the system, that is, the rate of attenuation of its free oscillations.
The criterion ξ of damping of oscillations over a period can be applied to the estimation of the system stability margin. The smaller ξ, the greater the transient response will have an oscillation, and the system will have a smaller margin of stability. It is believed that a dynamic system is well damped if ξ ≥ 90% [29].
At the beginning, using the algorithm for calculating the values of the transfer function, the degree difference is determined by the TF polynomial p = n − m based on (34).
Further calculations are carried out using the following iterative process.
Step 1. Put i= 1 and m = 1, η 0 = inf, ξ 0 = inf, where inf is a large number (for example, inf = 10 10 ), set the accuracy of determining the degree of stability ε η and the damping of oscillations for the period ε ξ . Step 2. Calculate n = p + m and, after performing rational interpolation find the vector a of the CP coefficients.
Step 3. Determine the roots of the characteristic equation, find among them the root with the largest real part and calculate the criteria η i and ξ i .
Step 4. Verify that the iterative process converges to a solution Step 5. If conditions (12) are fulfilled, then the quality criteria of the system dynamics are determined with the required accuracy, otherwise the process should be continued. To do this, increase the values of the iteration counter i and degree m by one and go to step 2.
The method described above was used to calculate the quality criteria of the dynamics of the bearing under consideration.

Bearing Dynamic Characteristics and Discussion
Among the parameters that do not affect the static characteristics, but affect the dynamic characteristics of the bearing, are the "compression number" σ and the volume of the intermediate cavity V p [23]. The influence of these parameters is of particular interest, because they are a resource to optimize the quality of a dynamic system. Table 1 shows the results of bearing optimization according to the performance criterion, which consisted in finding such pairs (σ, V p ) whose values deliver the maximum degree of stability η. Optimization was carried out for the values of other input parameters: P s = 5, R c = 0.5, χ = 0.48, ς = 0.3, H = 1. The parameter iK was varied, with the help of which the static compliance K was calculated by the formula where K 0 is the static compliance of the bearing corresponding to hard orifice 3 (K e = 0). Corresponding K values of the coefficient K e were calculated by the formula (30). An analysis of the calculation results showed that for each set of source data there is exactly one pair of values of the parameters σ, V p that deliver the maximum performance criterion η. From Table 1 it follows that with a decrease in compliance, the optimal σ tends to increase, while the optimal V p tends to decrease. With an increase in K e , the performance of the bearing decreases somewhat, but it remains stable even with negative values of K that are in absolute value superior to those of a conventional bearing (K e = 0).
A detailed idea of the influence of the parameters σ, V p on the degree of stability is given by the graph curves of Figure 11, which are constructed for the same initial data and K e = 0.32, at which the bearing has zero static compliance (K = 0). For σ < σmin, the system is unstable, and for σ > σmin the stability region is divided into two parts σmin < σ < σopt and σ > σopt, where σopt is the value of the parameter σ at which the function η(σ) reaches its maximum.
It was found that in the first part of this region, the transient characteristics are always oscillatory in nature, and in the second, they can be both oscillatory and aperiodic.
This follows from the analysis of graphs in Figure 12, where the dependences of the decay criterion for the period ξ on the compression number σ are presented. For Vp < Vpopt, the criterion ξ indicates the oscillatory nature of the transients (ξ < 100%), where Vp opt are the values of the parameter Vp at which the function η(Vp) reaches its maximum. The graphs show that for Vp < Vp opt all the curves, regardless of the value of σ, indicate the vibration nature of the transients (ξ < 100%), for Vp > Vpopt the transients become aperiodic (ξ = 100%). It is also seen that with an increase in Vp, the range of aperiodicity of the transient characteristics expands, that is, an increase in the volume of Vp contributes to an increase in the dynamic stability margin of the bearing as a dynamic system.  The graphs confirm the conclusion that the dependences η(σ) and η(V p ) are extreme. It can be seen that for the regime of zero compliance, it is necessary to choose the values σ > σ min , where σ min is the value of the parameter σ at which the system reaches the stability boundary (η = 0).
For σ < σ min , the system is unstable, and for σ > σ min the stability region is divided into two parts σ min < σ < σ opt and σ > σ opt , where σ opt is the value of the parameter σ at which the function η(σ) reaches its maximum.
It was found that in the first part of this region, the transient characteristics are always oscillatory in nature, and in the second, they can be both oscillatory and aperiodic.
This follows from the analysis of graphs in Figure 12, where the dependences of the decay criterion for the period ξ on the compression number σ are presented. For V p < V popt , the criterion ξ indicates the oscillatory nature of the transients (ξ < 100%), where V p opt are the values of the parameter Vp at which the function η(V p ) reaches its maximum. The graphs show that for V p < V p opt all the curves, regardless of the value of σ, indicate the vibration nature of the transients (ξ < 100%), for V p > V popt the transients become aperiodic (ξ = 100%). It is also seen that with an increase in V p , the range of aperiodicity of the transient characteristics expands, that is, an increase in the volume of V p contributes to an increase in the dynamic stability margin of the bearing as a dynamic system. Figure 11. Dependences of the degree of stability η on the "compression number" σ for various values of the volume Vp.  The data obtained allow us to give recommendations on the choice of values of the parameters σ and V p . The values of σ and V p should be considered the best, which are slightly larger than the optimal σ opt and V popt . For example, for the mode of zero static compliance K = 0, σ = 14 and V p = 7 will be the best. In this case, the bearing will have close to maximum speed and an aperiodic nature of the transition characteristics.

Conclusions
A mathematical model of statics and dynamics of a thrust gas-static bearing with an elastic orifice was formulated, a model for the deformation of an elastic orifice and a method to calculate the radius of its throttling diaphragm are proposed, a fast numerical method to calculate the quality indicators of the dynamics of a linear system with distributed parameters was developed, methods for determining the static and dynamic characteristics of the bearing were developed.
Calculations using the proposed model of elastic orifice deformation showed that its characteristics have full qualitative and satisfactory quantitative agreement with the results of experiment [22]. The flaws of the model were revealed and it was concluded that the proposed model and the method described can be used to a first approximation to calculate the radius of the supply diaphragm of the elastic orifice. The results of the study of the deformation of the elastic orifice are the basis for the model of the static and dynamic state of the bearing.
A formula was derived that allows one to calculate the value of the elastic ratio at which the bearing reaches zero compliance, and a formula to calculate the value of this parameter, at which the bearing loses its static stability.
The use of the fast numerical method to calculate the bearing dynamics quality criteria made it possible to repeatedly increase the speed of the multi-parameter design optimization algorithm by the performance criterion. A numerical experiment showed that, compared with the previously used method, the speed of the algorithm increased by about 12 times.
It was shown that the use of an elastic orifice allows a multiple reduction of positive compliance of the bearing, and also allows reduction of the compliance of the bearing to zero and even negative values.
An analysis of the quality criteria of structural dynamics confirms that with an optimal choice of the compression number σ and the volume of the inter-throttle chamber V p of the throttle chamber, the bearing can provide high speed and a guaranteed stability margin, while keeping the transition characteristics non-oscillatory.
Author Contributions: Conceptualization, validation, and supervision: S.S.; methodology, V.K. and S.S.; software, writing-review and editing, and project administration: V.K.; formal analysis, S.S. and A.K.; investigation and resources: Y.P.; data curation, V.K., Y.P. and A.K.; writing-original draft preparation, V.K. and A.K.; visualization, A.K. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
The work was carried out within the framework of the scientific research budget themes «Methods of modeling and optimizing of quality control information systems on the basis of intellectual technologies» at the Department of Standardization, Metrology and Quality Management and Department of Design and Technology Support for Engineering of the Polytechnic Institute of the Siberian Federal University.

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

CP
characteristic polynomial E elastic modulus of the orifice material f, F dimensional and dimensionless bearing loads h, H dimensional and dimensionless gaps K dimensionless bearing compliance K e dimensionless elastic ratio n, m transfer function polynomial orders p k , P k dimensional and dimensionless air pressures at the exit of annular orifice plates p p , P p dimensional and dimensionless air pressures at the chamber 4 p s , P s dimensional and dimensionless air supply pressures Q h , Q p dimensionless mass flow rates r, R dimensional and dimensionless radii r 1 , R 1 dimensional and dimensionless radii of arrangement of annular diaphragms r p , R p dimensional and dimensionless radii of elastic orifice diaphragm r p0 , R p0 dimensional and dimensionless radii of elastic orifice s Laplace transform variable TF transfer function u,U dimensional and dimensionless tensile functions of the middle surface of the elastic orifice v p , V p dimensional and dimensionless values of chamber 4 ∆H, ∆F Laplace transformants of the deviation of dynamic functions from their static values η degree of stability ν Poisson's ratio ξ damping of oscillations over a period Π Prandtl expiration function σ "compression number" ς coefficient of resistance adjustment of damping annular diaphragms ϕ, Φ dimensional and dimensionless Airy stress functions χ elastic orifice resistance adjustment coefficient ψ, Ψ dimensional and dimensionless deflection functions