$$\varepsilon^{2}\operatorname{Re}\tilde{\rho}\left(\frac{\partial\mathfrak{u}\_{\boldsymbol{x}}}{\partial\boldsymbol{\tilde{t}}}+\mathfrak{u}\_{\boldsymbol{x}}\frac{\partial\mathfrak{u}\_{\boldsymbol{x}}}{\partial\boldsymbol{\tilde{x}}}+\mathfrak{u}\_{\boldsymbol{y}}\frac{\partial\mathfrak{u}\_{\boldsymbol{x}}}{\partial\boldsymbol{\tilde{y}}}+\mathfrak{u}\_{\boldsymbol{z}}\frac{\partial\mathfrak{u}\_{\boldsymbol{x}}}{\partial\boldsymbol{\tilde{z}}}\right) = \frac{\varepsilon^{2}\operatorname{Re}}{\operatorname{Fr}^{2}}\tilde{\rho}\tilde{f}\_{\boldsymbol{x}}-\frac{\partial\tilde{\rho}}{\partial\boldsymbol{\tilde{x}}}+\frac{\partial}{\partial\boldsymbol{\tilde{y}}}\left(\tilde{\mu}\frac{\partial\mathfrak{u}\_{\boldsymbol{x}}}{\partial\boldsymbol{\tilde{y}}}\right), \tag{3.18a}$$

$$\varepsilon^{4}\operatorname{Re}\tilde{\rho}\left(\frac{\partial\tilde{u}\_{y}}{\partial\tilde{t}}+\tilde{u}\_{x}\frac{\partial\tilde{u}\_{y}}{\partial\tilde{x}}+\tilde{u}\_{y}\frac{\partial\tilde{u}\_{y}}{\partial\tilde{y}}+\tilde{u}\_{z}\frac{\partial\tilde{u}\_{y}}{\partial\tilde{z}}\right)=\frac{\varepsilon^{3}\operatorname{Re}}{\operatorname{Fr}^{2}}\tilde{\rho}\tilde{f}\_{y}-\frac{\partial\tilde{p}}{\partial\tilde{y}}.\tag{3.18b}$$

$$\varepsilon^{2}\operatorname{Re}\tilde{\rho}\left(\frac{\partial\mathfrak{u}\_{z}}{\partial\tilde{t}}+\mathfrak{u}\_{x}\frac{\partial\mathfrak{u}\_{z}}{\partial\tilde{x}}+\mathfrak{u}\_{y}\frac{\partial\mathfrak{u}\_{z}}{\partial\tilde{y}}+\mathfrak{u}\_{z}\frac{\partial\mathfrak{u}\_{z}}{\partial\tilde{z}}\right)=\frac{\varepsilon^{2}\operatorname{Re}}{\operatorname{Fr}^{2}}\tilde{\rho}f\_{z}-\frac{\partial\tilde{\rho}}{\partial\tilde{z}}+\frac{\partial}{\partial\tilde{y}}\left(\mathfrak{f}\frac{\partial\mathfrak{u}\_{z}}{\partial\tilde{y}}\right).\tag{3.18c}$$

Among the remaining terms, the relative importance of inertial forces and of body forces is then determined by a characteristic Reynolds number and a characteristic Froude number that are, respectively, defined as

$$\text{Re} = \frac{\rho\_0 \ell\_0^2 \omega\_0}{\mu\_0} \,\,\,\,\text{er} \qquad\qquad\qquad\text{Fr} = \omega\_0 \sqrt{\frac{\ell\_0}{\mathcal{g}}} \,\,\,\,\,\tag{3.19}$$

Looking at typical values for the involved physical quantities such as introduced in Chapter 6, the scaled Reynolds number *ε* 2 Re in both Eqs. (3.18a) and (3.18c) ranges between 10−2 and 10−1 , whereas *ε* 2 Re Fr−2 is located in the region of 10−5 to 10−4 . With the exception of excessively high rotational speeds, this confirms the essential hypothesis by Reynolds [Rey86] that inertial forces and body forces acting upon a thin fluid film appear insignificant compared to pressure forces and viscous forces.

Finally, the coefficients *ε* 4 Re and *ε* 3 Re Fr−2 in Eq. (3.18b) being smaller than even *ε* 2 , the back-transformed expression that remains from the original PDE simply writes

$$\frac{\partial p}{\partial y} = 0 \tag{3.20}$$

and suggests that the pressure field *p* is independent of the cross-film coordinate *y*. In classical lubrication theory, it is additionally hypothesized that not only pressure

but also density, viscosity, temperature, and other properties are constant across a sufficiently thin fluid film (see, e.g., [Cop49]). This simplification is overcome, at least in theory, with the generalized Reynolds equation proposed by Dowson [Dow62], which accounts for possible cross-film variations of density and viscosity. In practice, however, it can be extremely cumbersome or even impossible to rigorously describe the transient behavior of an entire bearing – rotor system starting from such a partial integro-differential equation. In the literature, cross-film variations of properties are considered almost exclusively in very detailed steady-state analyses [Gar12, MLS16], whereas comparable rotordynamic studies [Bou17, LMS18] understand especially density as cross-film averaged quantity. As recently outlined by Lehn et al. [LMS18], the generalized Reynolds equation has of course its justification, but rather when applied to oil-lubricated bearings because the magnitude of non-uniformities across gas films is often insignificant. With this in mind, the present thesis supposes that the density field is a priori independent of the cross-film coordinate in analogy to what Eq. (3.20) indicates for pressure. Having said this, the thermodynamic state principle mentioned in Section 3.1.1 requires that, since the two properties pressure and density do not depend on *y*, any other properties such as viscosity do neither. Setting accordingly

$$\frac{\partial \rho}{\partial y} = 0,\tag{3.21}$$

$$\frac{\partial \mu}{\partial y} = 0,\tag{3.21}$$

the back-transformed terms remaining from Eqs. (3.18a) and (3.18c) simplify to give

$$\frac{\partial p}{\partial x} = \frac{\partial}{\partial y} \left( \mu \frac{\partial u\_x}{\partial y} \right) = \mu \frac{\partial^2 u\_x}{\partial y^2},\tag{3.22a}$$

$$\frac{\partial p}{\partial z} = \frac{\partial}{\partial y} \left( \mu \frac{\partial u\_z}{\partial y} \right) = \mu \frac{\partial^2 u\_z}{\partial y^2}. \tag{3.22b}$$

To actually obtain flow velocity profiles depending on prevailing pressure gradients, both Eqs. (3.22a) and (3.22b) must be integrated twice with respect to *y* and can then be adapted to boundary conditions that are related to the motion of adjacent walls.

On the inner fluid film boundary, which is found at *y* = −*h*∗ according to Eq. (2.14), the tangential component of the wall velocity is dominated by the rotation rather than by the translation of the shaft journal and points approximately in the **e***x*-direction for a sufficiently thin film (see, e.g., [Sze10]). Together with the simplifications that any tangential motion of the foil structure on the outer boundary at *y* = −*q* appears negligible [Bau18] and that neither of the bounding walls moves axially, this yields

$$\mu\_{\mathbf{x}}\Big|\_{\mathbf{y}=-h\_{\bullet}} = \mathrm{R}\omega\_{\mathbf{0}}\qquad\qquad\qquad\qquad\mu\_{\mathbf{x}}\Big|\_{\mathbf{y}=-q} = \mathbf{0},\tag{3.23a}$$

$$\left.u\_{z}\right|\_{y=-h\_{\*}} = 0, \qquad\qquad\qquad\qquad\left.u\_{z}\right|\_{y=-q} = 0.\tag{3.2\text{\textquotedbl{}}}$$

Possible refinements of such no-slip boundary conditions, which could be extended to consider slip flow as proposed by Burgdorfer [Bur59], are beyond the scope of the present thesis. Nevertheless, it should be kept in mind that there are possible applications with elevated Knudsen numbers, in which no-slip boundary conditions might become inaccurate because of beginning molecular rarefaction effects [HD83]. Here, keeping simply with Eq. (3.23a), the circumferential flow velocity field

$$
\mu\_{\rm x} = \frac{1}{2\mu} \frac{\partial p}{\partial \mathbf{x}} (h\_\* + y) \left(q + y\right) - \mathcal{R}\omega\_0 \left(\frac{q + y}{h\_\* - q}\right) \tag{3.24}
$$

matches the prescribed velocity *Rω*0 on the inner fluid film boundary at *y* = −*h*∗ and vanishes on the outer boundary at *y* = −*q*, whereas the axial flow velocity field

$$
\mu\_z = \frac{1}{2\mu} \frac{\partial p}{\partial z} (h\_\* + y) \left(q + y\right) \tag{3.25}
$$

satisfies Eq. (3.23b) and is zero at both boundaries. Generally speaking, the effective circumferential flow that is described by Eq. (3.24) results from the superposition of a pressure-driven Poiseuille flow and a shear-driven Couette flow, whereas the axial flow described by Eq. (3.25) corresponds to a purely pressure-driven Poiseuille flow.

#### **3.2.2 Integration of Continuity Equation**

To proceed with the derivation of an appropriate lubricant model, the Navier–Stokes equations from Eq. (3.14) are completed by the continuity equation for compressible fluids in the general form

$$\frac{\mathbf{D}\rho}{\mathbf{D}t} + \rho \operatorname{div} \mathbf{u} = \mathbf{0},\tag{3.26}$$

which again is found in all common textbooks on fluid mechanics (see, e.g., [Hue98]).

With the definition of the material derivative from Eq. (3.15) and with a generalized product rule, the PDE from Eq. (3.26) can be stated in Cartesian coordinates as

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} \left( \rho u\_x \right) + \frac{\partial}{\partial y} \left( \rho u\_y \right) + \frac{\partial}{\partial z} \left( \rho u\_z \right) = 0,\tag{3.27}$$

where all terms are of the same order of magnitude if the characteristic length and velocity scales from Eqs. (3.16a) and (3.16b) hold. Since no further information about the radial flow velocity *uy* can here be obtained from the Navier–Stokes equations, an essential step consists in the integration of Eq. (3.27) across the fluid film, giving

$$\int\_{-h\_\*}^{-q} \frac{\partial \rho}{\partial t} \, \mathrm{d}y + \int\_{-h\_\*}^{-q} \frac{\partial}{\partial \mathbf{x}} \left(\rho u\_x\right) \, \mathrm{d}y + \left(\rho u\_y\right) \Big|\_{-h\_\*}^{-q} + \int\_{-h\_\*}^{-q} \frac{\partial}{\partial \mathbf{z}} \left(\rho u\_z\right) \, \mathrm{d}y = 0. \tag{3.28}$$

31

In Eq. (3.28), the unknown *uy* appears only at the boundaries of the fluid film, where it can be directly related to the motion of the shaft journal and of the foil structure.

The normal components of the wall velocities at the respective boundaries, which point both approximately in the **e***y*-direction for a sufficiently thin film, are given by

$$\left.u\_y\right|\_{y=-h\_\*} = -\frac{\partial h\_\*}{\partial t} - \mathcal{R}\omega\_0 \frac{\partial h\_\*}{\partial \mathbf{x}}\,\prime \qquad\qquad \qquad \left.u\_y\right|\_{y=-q} = -\frac{\partial q}{\partial t} \tag{3.29}$$

and consider the rotation and the translation of the shaft journal (see, e.g., [Sze10]) as well as the deformation rate of the foil structure [Bau18]. Since the density field *ρ* is supposed to be constant across the film according to Eq. (3.21), the third term on the left-hand side of Eq. (3.28) can directly be expressed with the help of Eq. (3.29) as

$$\rho \left( \rho u\_{\mathcal{Y}} \right) \Big|\_{ -h\_\*}^{ -q } = \rho \left( u\_{\mathcal{Y}} \Big|\_{y=-q} - u\_{\mathcal{Y}} \Big|\_{y=-h\_\*} \right) = -\rho \frac{\partial q}{\partial t} + \rho \frac{\partial h\_\*}{\partial t} + \rho \text{R}\omega\_0 \frac{\partial h\_\*}{\partial \mathbf{x}} \,. \tag{3.20}$$

With Leibniz' rule for differentiation under the integral sign, which is obviously applicable to each of the three integrals on the left-hand side of Eq. (3.28), it comes

$$\int\_{-h\_\*}^{-q} \frac{\partial \rho}{\partial t} \, \mathbf{d}y = \frac{\partial}{\partial t} \left( \rho \int\_{-h\_\*}^{-q} \mathbf{d}y \right) + \rho \frac{\partial q}{\partial t} - \rho \frac{\partial h\_\*}{\partial t} \, \mathbf{,}\tag{3.31a}$$

$$\int\_{-h\_{\ast}}^{-q} \frac{\partial}{\partial \mathfrak{x}} \left(\rho u\_{\mathfrak{x}}\right) \mathrm{d}y = \frac{\partial}{\partial \mathfrak{x}} \left(\rho \int\_{-h\_{\ast}}^{-q} u\_{\mathfrak{x}} \mathrm{d}y\right) - \rho \mathrm{R}\omega\_{0} \frac{\partial h\_{\ast}}{\partial \mathfrak{x}} \; , \tag{3.31b}$$

$$\int\_{-\mathbf{h}\_{\ast}}^{-q} \frac{\partial}{\partial z} \left( \rho u\_{z} \right) \, \mathrm{d}y = \frac{\partial}{\partial z} \left( \rho \int\_{-\mathbf{h}\_{\ast}}^{-q} u\_{z} \, \mathrm{d}y \right), \tag{3.31c}$$

where density is constant across the film in accordance with Eq. (3.21) and velocities at the integration limits are obtained from Eqs. (3.23a) and (3.23b), respectively.

When finally rewriting Eq. (3.28) using these conversions, all terms from Eq. (3.30) cancel out with the corresponding terms in Eqs. (3.31a) and (3.31b), so it only remains

$$\frac{\partial}{\partial t} \left( \rho \int\_{-h\_\*}^{-q} \mathbf{d}y \right) + \frac{\partial}{\partial \mathbf{x}} \left( \rho \int\_{-h\_\*}^{-q} u\_{\mathbf{x}} \, \mathbf{d}y \right) + \frac{\partial}{\partial \mathbf{z}} \left( \rho \int\_{-h\_\*}^{-q} u\_{\mathbf{z}} \, \mathbf{d}y \right) = 0. \tag{3.32}$$

Remembering the relation *h* = *h*∗ − *q* from Eq. (2.13c), the trivial integral in the first term of Eq. (3.32) simply yields the fluid film thickness field *h* itself. Additionally, the two remaining integrals can be compactly expressed as *hu*¯*x* and *hu*¯*z* after calculating from Eqs. (3.24) and (3.25) the cross-film averaged bulk flow velocity fields

$$\mathfrak{u}\_{\mathbf{x}} = \frac{1}{h} \int\_{-h\_{\mathbf{+}}}^{-q} u\_{\mathbf{x}} \mathbf{d}y = -\frac{h^2}{12\mu} \frac{\partial p}{\partial \mathbf{x}} + \frac{\mathbf{R}\omega\_0}{2} = \mathfrak{u}\_{\mathbf{x},\mathbf{P}} + \mathfrak{u}\_{\mathbf{x},\mathbf{C}} \tag{3.33a}$$

$$
\bar{u}\_z = \frac{1}{h} \int\_{-h\_\*}^{-q} u\_z \, \mathrm{d}y = -\frac{h^2}{12\mu} \frac{\partial p}{\partial z} \,. \tag{3.32b}
$$

According to Eq. (3.33a), the expression for the circumferential velocity field can be split into a first term *u*¯*x*,P representing the Poiseuille flow and a second term *u*¯*x*,C representing the Couette flow. In the case of a pure Poiseuille flow as in Eq. (3.33b) for the axial direction, such a decomposition of the velocity field is of course needless.

To obtain a more descriptive and more common form of the resulting shortform PDE

$$\frac{\partial}{\partial t} \left( \rho h \right) + \frac{\partial}{\partial \mathbf{x}} \left( \rho h \overline{u}\_x \right) + \frac{\partial}{\partial z} \left( \rho h \overline{u}\_z \right) = 0,\tag{3.34}$$

the bulk flow velocities from Eqs. (3.33a) and (3.33b) are explicitly inserted to yield

$$\frac{\partial}{\partial t} \left( \rho h \right) + \frac{\text{R} \omega\_0}{2} \frac{\partial}{\partial \mathbf{x}} \left( \rho h \right) = \frac{\partial}{\partial \mathbf{x}} \left( \frac{\rho h^3}{12\mu} \frac{\partial p}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial z} \left( \frac{\rho h^3}{12\mu} \frac{\partial p}{\partial z} \right) \tag{3.35}$$

or, making use of the arc-length transformation *x* = *Rϕ* from Eq. (2.6), to finally give

$$\frac{\partial}{\partial t} \left( \rho h \right) + \frac{\omega\_0}{2} \frac{\partial}{\partial \varphi} \left( \rho h \right) = \frac{1}{R^2} \frac{\partial}{\partial \varphi} \left( \frac{\rho h^3}{12\mu} \frac{\partial p}{\partial \varphi} \right) + \frac{\partial}{\partial z} \left( \frac{\rho h^3}{12\mu} \frac{\partial p}{\partial z} \right). \tag{3.36}$$

Based on the fundamental work of Reynolds [Rey86] about incompressible lubricants, this PDE is subsequently referred to as the Reynolds equation for compressible fluids.

#### **3.2.3 Closure of Coupled Equation System**

In the most general form of Eq. (3.36) to be considered, the involved field quantities

$$p = p(\boldsymbol{\varrho}, \boldsymbol{z}, \mathbf{t}), \qquad \boldsymbol{\varrho} = \boldsymbol{\varrho}(\boldsymbol{\varrho}, \mathbf{z}, \mathbf{t}), \qquad \boldsymbol{\mu} = \boldsymbol{\mu}(\boldsymbol{\varrho}, \mathbf{z}, \mathbf{t}), \qquad \boldsymbol{h} = \boldsymbol{h}(\boldsymbol{\varrho}, \mathbf{t})\tag{3.37}$$

may depend on the time *t* as well as on the circumferential coordinate *ϕ* and, with the exception of the fluid film thickness in the last expression, on the axial coordinate *z*. With regard to the number of unknown fields in Eq. (3.37), the Reynolds equation needs to be closed by some constitutive laws for the lubricant, which are a thermal equation of state to eliminate the pressure field *p*(*ϕ*, *z*, *t*) according to Eq. (3.2) and a viscous equation of state as in Eq. (3.7) to express the viscosity field *µ*(*ϕ*, *z*, *t*). Concerning the refrigerant R-245fa, for which the perfect gas assumption according to Eqs. (3.9) and (3.13) is clearly insufficient, the consideration of equations of state as shown in Figs. 3.2a and 3.2e adds an important nonlinearity to the resulting PDE. Moreover, as explained in Section 2.3, it is finally the insertion of *h*(*ϕ*, *t*) from either Eq. (2.13a) or Eq. (2.13b) into Eq. (3.36) that introduces the coupling quantities from the foil structure submodel and from the rotor submodel into the fluid submodel. Provided the coupling quantities are adequately described by the equations of the other two submodels, the Reynolds equation is ultimately formulated only in terms

of the selected primary variables, which are the density field *ρ*(*ϕ*, *z*, *t*) and the temperature field *ϑ*(*ϕ*, *z*, *t*). Of course, a PDE that physically describes the evolution of the latter is missing up to this point, which is the motivation for additionally deriving a temperature equation in Section 3.3 based on the internal energy.

In connection with the Reynolds number2 from Eq. (3.19) and notwithstanding the assumption of negligible inertial forces, it remains to clarify the role of fluid instability mechanisms that lead to the formation of Taylor vortices or turbulent flow. Already in the pioneering theoretical investigation of turbulent gas lubrication by Constantinescu [Con64], it actually turns out that the influence of such phenomena is only of a quantitative nature and less significant in gas bearings than in oil bearings. For this reason and for the sake of simplicity, the widely used class of zero-equation turbulence models leaves governing PDEs such as the Reynolds equation basically unchanged3 and introduces an additional turbulent eddy viscosity instead, which is algebraically related to quantities such as flow velocity gradients (see, e.g., [FPS20]). With such an approach, more precisely with the Ng –Pan –Elrod model [NP65, EN67], the numerical investigation by Garcia [Gar12] reveals that turbulence does occur in refrigerant-lubricated GFBs, but typically at locations that are largely insignificant for their overall performance. Although effects on bearing–rotor systems are thus limited, it is noted that the Reynolds equation in Eq. (3.36) principally proves suitable for a qualitative, but also for a quantitative description of turbulent flow provided the viscosity field is altered appropriately. Unsurprisingly, such a model extension would be relatively straightforward [Gar12, Bou17], but is not pursued further here to keep the number of empirical parameters small when analyzing qualitative trends.

## **3.3 Derivation of Thin-Film Temperature Equation**

While it might be sufficient to generally assume isothermal conditions for the analysis of GFBs lubricated with ambient air, this is certainly not a reasonable supposition when dealing with refrigerant-lubricated GFBs [Gar12]. Consequently, the continuum hypothesis according to Eq. (3.1) and the thin-film assumption from Eq. (2.1) are used hereafter to derive a temperature equation similar to the Reynolds equation.

2 To be precise, the characteristic Reynolds number must be completed by a characteristic Taylor number to correctly distinguish between parallel flow instabilities and centrifugal instabilities (see, e.g., [Sze10]).

3 The form of the Reynolds equation for turbulent flow is unchanged since the underlying Reynoldsaveraged Navier–Stokes equations also have the same form as the original equations (see, e.g., [FPS20]).

#### **3.3.1 Approximated Internal Energy Equation**

According to usual textbooks on fluid mechanics (see, e.g., [Hue98]), the combination of the Navier–Stokes equations from Eq. (3.14) with the continuity equation from Eq. (3.26) is typically accompanied by the internal energy equation for compressible Newtonian fluids in the general form

$$\rho \frac{\text{DE}}{\text{D}t} = -p \operatorname{div} \mathbf{u} - \operatorname{div} \boldsymbol{\Phi} + \lambda \left( \operatorname{div} \mathbf{u} \right)^2 + 2\mu \left\| \frac{\operatorname{grad} \mathbf{u} + \operatorname{grad}^{\mathsf{T}} \mathbf{u}}{2} \right\|^2. \tag{3.38}$$

Using the material derivative from Eq. (3.15), this PDE expresses the principle of conservation of energy by relating changes within the specific internal energy field *E* to different energy sources and energy sinks. While the first term of the right-hand side of Eq. (3.38) is associated with compression work, the second term contains the heat flux vector field φ and is therefore supposed to represent thermal conduction. Most importantly, the last two terms4 correspond to viscous dissipation work, which depends on both the dynamic viscosity field *µ* and the volume viscosity field *λ*.

Transferring the reformulated exact differential d*E* that is given by Eq. (3.4c) into the corresponding material derivative D*E*/ D*t* and then substituting D*ρ*/ D*t* = −*ρ* div **u** by virtue of the continuity equation from Eq. (3.26), it comes immediately

$$
\rho \frac{\text{DE}}{\text{D}t} = \rho c\_{\text{V}} \frac{\text{D}\vartheta}{\text{D}t} + \frac{p}{\rho} (1 - \alpha\_{\text{V}}\vartheta) \frac{\text{D}\rho}{\text{D}t} \tag{3.39a}
$$

$$\mathbf{u} = \rho \mathbf{c}\_{\mathrm{V}} \frac{\mathrm{D}\theta}{\mathrm{D}t} - p \,\mathrm{div}\,\mathbf{u} + \alpha\_{\mathrm{V}} \theta p \,\mathrm{div}\,\mathbf{u}.\tag{3.39b}$$

With this relationship, the internal energy equation from Eq. (3.38) takes a convenient temperature-explicit form to describe changes within the temperature field *ϑ* by

$$\rho c\_{\rm v} \frac{\mathbf{D}\theta}{\mathbf{D}t} = -a\_{\rm v} \theta p \,\mathrm{div}\,\mathbf{u} - \mathrm{div}\,\boldsymbol{\Phi} + \lambda \left(\mathrm{div}\,\mathbf{u}\right)^{2} + 2\mu \left\| \frac{\mathrm{grad}\,\mathbf{u} + \mathrm{grad}^{\mathsf{T}}\mathbf{u}}{2} \right\|^{2}.\tag{3.40}$$

This reformulated PDE now involves the specific isochoric heat capacity field *c*v on the left-hand side and the isochoric compressibility coefficient field *α*v in the first term on the right-hand side, whereas contributions due to thermal conduction and dissipation work remain unchanged. An alternative temperature-explicit form of this PDE is actually shown in [LBS19], where the specific isobaric heat capacity and the isobaric compressibility coefficient occur instead of isochoric properties.

4 The mathematical expression for the viscous dissipation work makes use of the Frobenius norm k. . .k, which is calculated as the square root of the sum of the squared tensor components (see, e.g., [GV13]).

However, such a derivation via the specific enthalpy leads to an appearance of the time derivative of pressure, which cannot be substituted by virtue of the continuity equation and thus makes the solution of transient problems much less intuitive.

In the following step, in which the PDE from Eq. (3.40) is finally represented in the Cartesian coordinates (*x*, *y*, *z*) that describe the lubrication gap, the substitutions from Eqs. (3.16a) and (3.16b) allow again for a normalization of all coordinates and flow velocity components. Moreover, the various scalar fields of fluid properties as well as the heat flux components *φx*, *φy*, and *φz* can be transformed to normalized variables (denoted again with a tilde) by means of Eq. (3.17a) and by setting

$$\mathfrak{c}\_{\rm V} = \mathfrak{c}\_{0}\mathfrak{x}\_{\rm V} \qquad \qquad \mathfrak{d} = \frac{\mu\_{0}\omega\_{0}}{\varepsilon^{2}\rho\_{0}c\_{0}}\mathfrak{d}\_{\rm V} \qquad \qquad \mathfrak{a}\_{\rm V} = \frac{\varepsilon^{2}\rho\_{0}c\_{0}}{\mu\_{0}\omega\_{0}}\mathfrak{x}\_{\rm V} \tag{3.41a}$$

$$\phi\_{\mathbf{x}} = \mathfrak{c}\_{0}\mu\_{0}\frac{\Delta\theta\_{0}}{\ell\_{0}}\tilde{\phi}\_{\mathbf{x}} \qquad \qquad \phi\_{\mathbf{y}} = \mathfrak{c}\_{0}\mu\_{0}\frac{\Delta\theta\_{0}}{\varepsilon\ell\_{0}}\tilde{\phi}\_{\mathbf{y}} \qquad \qquad \phi\_{\mathbf{z}} = \mathfrak{c}\_{0}\mu\_{0}\frac{\Delta\theta\_{0}}{\ell\_{0}}\tilde{\phi}\_{\mathbf{z}}.\tag{3.41b}$$

After introducing in Eq. (3.41a) a characteristic specific heat capacity *c*0, it is natural to construct the temperature scale *µ*0*ω*0*ε* −2*ρ* −1 0 *c* −1 0 and, reciprocally, the compressibility coefficient scale, in analogy to the adopted viscous pressure scale (see, e.g., [Cop49]). In Eq. (3.41b), assuming that heat fluxes generally obey Fourier's law φ ∝ − grad *ϑ*, the constant of proportionality is of the order of *c*0*µ*0 according to the kinetic theory of gases (see, e.g., [Ken38]) and the temperature gradient may be assessed by relating a characteristic temperature difference *∆ϑ*0 to the respective length scales `<sup>0</sup> and *ε*`0.

From a subsequent order-of-magnitude estimation for the normalized internal energy equation in Cartesian coordinates, it turns out that most of the terms related to viscous dissipation work are of orders greater or equal to *ε* 2 and can thus be neglected. In doing so, it remains only the reduced PDE

$$
\begin{split}
\tilde{\rho}\tilde{c}\_{\text{v}}\left(\frac{\partial\tilde{\theta}}{\partial\tilde{t}} + \tilde{u}\_{\text{x}}\frac{\partial\tilde{\theta}}{\partial\tilde{\mathbf{x}}} + \tilde{u}\_{\text{y}}\frac{\partial\tilde{\theta}}{\partial\tilde{y}} + \tilde{u}\_{\text{z}}\frac{\partial\tilde{\theta}}{\partial\tilde{z}}\right) &= -\tilde{u}\_{\text{v}}\tilde{\theta}\tilde{\rho}\left(\frac{\partial\tilde{u}\_{\text{x}}}{\partial\tilde{\mathbf{x}}} + \frac{\partial\tilde{u}\_{\text{y}}}{\partial\tilde{y}} + \frac{\partial\tilde{u}\_{\text{z}}}{\partial\tilde{z}}\right) \\ & \quad - \frac{1}{\text{Ec}}\left(\varepsilon^{2}\frac{\partial\tilde{\phi}\_{\text{x}}}{\partial\tilde{x}} + \frac{\partial\tilde{\phi}\_{\text{y}}}{\partial\tilde{y}} + \varepsilon^{2}\frac{\partial\tilde{\phi}\_{\text{z}}}{\partial\tilde{z}}\right) + \tilde{\mu}\left[\left(\frac{\partial\tilde{u}\_{\text{x}}}{\partial\tilde{y}}\right)^{2} + \left(\frac{\partial\tilde{u}\_{\text{z}}}{\partial\tilde{y}}\right)^{2}\right],\tag{3.42}
\end{split}
$$

in which the relative importance of thermal conduction can directly be judged by a characteristic Eckert number defined as

$$\mathbf{Ec} = \frac{\ell\_0^2 \omega\_0^2}{c\_0 \Delta \theta\_0}.\tag{3.43}$$

Since Ec−1 is expected of the order of unity when inserting typical parameter values according to Chapter 6, the cross-film heat flux must clearly be retained in Eq. (3.42), while heat fluxes along the film are smaller by the factor of *ε* 2 and appear negligible. On account of the state principle and with an intent to describe internal energy at the same level of model complexity as pressure and density in Eqs. (3.20) and (3.21), it is assumed hereafter that temperature, specific heat capacity, and the compressibility coefficient are all independent of the cross-film coordinate *y*, in other words

$$
\frac{\partial \vartheta}{\partial y} = 0,\qquad\qquad\qquad\frac{\partial c\_{\rm v}}{\partial y} = 0,\qquad\qquad\qquad\frac{\partial a\_{\rm v}}{\partial y} = 0.\tag{3.44}
$$

Even though this simplification seems entirely justified for investigating qualitative effects that concern the overall system, it should of course not be ignored that the accuracy of quantitative predictions can still be improved by allowing temperature and other properties to vary across the fluid film [Gar12, Bou17]. In order not to overcomplicate the fluid model without necessity, however, the present thesis sticks to Eq. (3.44), especially since there is both analytical [Sny65] and numerical [Leh17] evidence that non-uniformities across gas films are much less significant than they are across liquid films. Ultimately, this makes the cross-film temperature gradient disappear on the left-hand side of Eq. (3.42) and the back-transformed PDE becomes

$$\begin{split} \rho c\_{\mathrm{v}} \left( \frac{\partial \theta}{\partial t} + u\_{x} \frac{\partial \theta}{\partial x} + u\_{z} \frac{\partial \theta}{\partial z} \right) &= -u\_{\mathrm{v}} \theta p \left( \frac{\partial u\_{x}}{\partial x} + \frac{\partial u\_{y}}{\partial y} + \frac{\partial u\_{z}}{\partial z} \right) \\ &\quad - \frac{\partial \phi\_{y}}{\partial y} + \mu \left[ \left( \frac{\partial u\_{x}}{\partial y} \right)^{2} + \left( \frac{\partial u\_{z}}{\partial y} \right)^{2} \right]. \end{split} \tag{3.45}$$

The flow velocity components and the retained heat flux component may still depend on *y* as they are not intensive thermodynamic properties obeying the state principle.

#### **3.3.2 Integration of Internal Energy Equation**

As no further information is retrieved from any equation about the fields *uy* and *φy*, the most convenient remedy is to integrate Eq. (3.45) across the film, which leads to

$$\begin{split} \int\_{-h\_{\ast}}^{-q} \rho c\_{\mathrm{V}} \left( \frac{\partial \theta}{\partial t} + u\_{\mathrm{x}} \frac{\partial \theta}{\partial \mathrm{x}} + u\_{\mathrm{z}} \frac{\partial \theta}{\partial \mathrm{z}} \right) \mathrm{d}y \\ &= - \int\_{-h\_{\ast}}^{-q} a\_{\mathrm{V}} \theta p \frac{\partial u\_{\mathrm{x}}}{\partial \mathrm{x}} \, \mathrm{d}y - \left( a\_{\mathrm{v}} \theta p u\_{\mathrm{y}} \right) \Big|\_{-h\_{\ast}}^{-q} - \int\_{-h\_{\ast}}^{-q} a\_{\mathrm{v}} \theta p \frac{\partial u\_{\mathrm{z}}}{\partial \mathrm{z}} \, \mathrm{d}y \\ &- \phi\_{\mathrm{y}} \Big|\_{-h\_{\ast}}^{-q} + \int\_{-h\_{\ast}}^{-q} \mu \left( \frac{\partial u\_{\mathrm{x}}}{\partial \mathrm{y}} \right)^{2} \mathrm{d}y + \int\_{-h\_{\ast}}^{-q} \mu \left( \frac{\partial u\_{\mathrm{z}}}{\partial \mathrm{y}} \right)^{2} \mathrm{d}y \end{split} \tag{3.46}$$

in the very same manner as previously seen in Eq. (3.28) for the continuity equation.

After integration, the unknown cross-film heat flux *φy* appears only at the boundaries of the fluid film, where it is prescribed by simple heat transfer boundary conditions

$$\left.\phi\_{y}\right|\_{y=-h\_{\ast}}=-\beta\_{0,\mathrm{R}}\left(\theta-\theta\_{\infty}\right),\qquad\left.\phi\_{y}\right|\_{y=-q}=\beta\_{0,\mathrm{S}}\left(\theta-\theta\_{\infty}\right).\tag{3.47}$$

Here, imagining the entire bearing – rotor system as a network of thermal resistances, heat fluxes from or into the environment are assumed proportional to the difference between the local fluid temperature *ϑ* and the constant ambient temperature *ϑ*∞. With this empirical approach (see, e.g., [ÇG15]), all contributions due to convection, conduction, and radiation are combined into two equivalent heat transfer coefficients, of which *β*0,R refers to heat fluxes through the rotor shaft and *β*0,S to heat fluxes through the foil structure and bearing sleeve. Even though this concept obviously represents a strong idealization of actual heat transfer mechanisms, the in-depth analysis of GFB temperatures by San Andrés et al. [SK10] shows that the complexity of such a model is largely sufficient to reproduce experimental results. Consequently, the actual difficulty arises rather from the correct determination of empirical heat transfer coefficients, which is of course not the primary interest of the present thesis. Finally, the first term in the third line of Eq. (3.46) can be compactly expressed as

$$\left.\phi\_{y}\right|\_{-h\_{\ast}}^{-q} = \left.\phi\_{y}\right|\_{y=-q} - \left.\phi\_{y}\right|\_{y=-h\_{\ast}} = \beta\_{0}\left(\theta - \theta\_{\infty}\right) \tag{3.48}$$

when summing up the heat transfer coefficients in Eq. (3.47) to give *β*0 = *β*0,S + *β*0,R.

The further evaluation of the integrated energy equation turns out to be relatively simple because, owing to Eqs. (3.20), (3.21), and (3.44), all remaining field quantities apart from the components of the flow velocity are independent of the *y*-coordinate. Having this in mind and starting with the most obvious part, the integral describing temperature changes in the first line of Eq. (3.46) immediately becomes

$$\begin{split} \int\_{-h\_{\ast}}^{-q} \rho c\_{\mathrm{V}} \left( \frac{\partial \theta}{\partial t} + u\_{\mathrm{x}} \frac{\partial \theta}{\partial \mathrm{x}} + u\_{\mathrm{z}} \frac{\partial \theta}{\partial \mathrm{z}} \right) \mathrm{d}y \\ = \rho c\_{\mathrm{V}} \left( \frac{\partial \theta}{\partial t} \int\_{-h\_{\ast}}^{-q} \mathrm{d}y + \frac{\partial \theta}{\partial \mathrm{x}} \int\_{-h\_{\ast}}^{-q} u\_{\mathrm{x}} \, \mathrm{d}y + \frac{\partial \theta}{\partial \mathrm{z}} \int\_{-h\_{\ast}}^{-q} u\_{\mathrm{z}} \, \mathrm{d}y \right). \end{split} \tag{3.49}$$

Furthermore, the terms related to compression work in the second line of Eq. (3.46) possess a structure that is somewhat similar to the integrated continuity equation, from which one benefits by slightly modifying Eqs. (3.30), (3.31b), and (3.31c) to get

$$\int\_{-h\_{\ast}}^{-q} \mathfrak{a}\_{\mathrm{V}} \mathfrak{d}p \frac{\partial \mathfrak{u}\_{\mathrm{x}}}{\partial \mathfrak{x}} \mathrm{d}y = \mathfrak{a}\_{\mathrm{V}} \mathfrak{d}p \frac{\partial}{\partial \mathfrak{x}} \int\_{-h\_{\ast}}^{-q} \mathfrak{a}\_{\mathrm{x}} \mathrm{d}y - \mathfrak{a}\_{\mathrm{V}} \mathfrak{d}p \mathrm{R}\omega\_{0} \frac{\partial \mathfrak{h}\_{\ast}}{\partial \mathfrak{x}} \,. \tag{3.50a}$$

$$\left(a\_{\rm V}\theta p u\_{\rm y}\right)\Big|\_{-\hbar\_{\rm \ast}}^{-q} = -a\_{\rm V}\theta p \frac{\partial q}{\partial t} + a\_{\rm V}\theta p \frac{\partial h\_{\rm \ast}}{\partial t} + a\_{\rm V}\theta p R\omega\_{0} \frac{\partial h\_{\rm \ast}}{\partial \mathbf{x}} \,. \tag{3.50b}$$

$$\int\_{-\hbar\_{\ast}}^{-q} \mathfrak{a}\_{\mathrm{V}} \mathfrak{d}p \frac{\partial \mathfrak{u}\_{\mathrm{z}}}{\partial \boldsymbol{z}} \, \mathrm{d}y = \mathfrak{a}\_{\mathrm{V}} \mathfrak{d}p \frac{\partial}{\partial \boldsymbol{z}} \int\_{-\hbar\_{\ast}}^{-q} \mathfrak{u}\_{\mathrm{z}} \, \mathrm{d}y. \tag{3.50c}$$

On closer inspection of the found expressions, all of Eqs. (3.49) and (3.50a) to (3.50c) can be stated in terms of the effective fluid film thickness *h* = *h*∗ − *q* from Eq. (2.13c) and of the bulk flow velocities *u*¯*x* and *u*¯*z* from Eqs. (3.33a) and (3.33b), respectively.

Coming to the integrals of viscous dissipation work in the third line of Eq. (3.46), however, no such simplification seems possible without prior knowledge about the involved flow velocity fields. Consequently, it is only after calculating the squared derivatives of Eqs. (3.24) and (3.25) that carrying out the actual integration yields

$$\int\_{-h\_{\ast}}^{-q} \mu \left(\frac{\partial u\_{\mathrm{x}}}{\partial y}\right)^{2} \mathrm{d}y = \overbrace{\frac{(h\_{\ast} - q)^{3}}{12\mu} \left(\frac{\partial p}{\partial \mathrm{x}}\right)^{2}}\_{\mathrm{-}} + \underbrace{\mathrm{R}^{2} \omega\_{0}^{2} \frac{\mu}{h\_{\ast} - q}}\_{\mathrm{-}}.\tag{3.51a}$$

$$\int\_{-h\_\*}^{-q} \mu \left(\frac{\partial u\_z}{\partial y}\right)^2 dy = \underbrace{\frac{\left(h\_\* - q\right)^3}{12\mu} \left(\frac{\partial p}{\partial z}\right)^2}\_{=\frac{12\mu}{h\_\* - q} a\_z^2}.\tag{3.51b}$$

Surprisingly, each summand in Eqs. (3.51a) and (3.51b) resembles the square of a particular bulk flow velocity contribution when setting *h* = *h*∗ − *q* and comparing the resulting expressions to *u*¯ 2 *x*,P and *u*¯ 2 *x*,C from Eq. (3.33a) and to *u*¯ 2 *z* from Eq. (3.33b). Since there is no general relationship between the integrated square of the derivative of an arbitrary function and the squared integrals of its individual summands, this is only a neat side effect for flow velocity profiles given by second degree polynomials.

Altogether, having substituted Eqs. (3.48), (3.49), (3.50a) to (3.50c), (3.51a), and (3.51b) into Eq. (3.46), the final temperature equation is given by the handy shortform PDE

$$\begin{split} \rho c\_{\rm v} h \Big( \frac{\partial \theta}{\partial t} + \bar{u}\_{\rm x} \frac{\partial \theta}{\partial \mathbf{x}} + \bar{u}\_{z} \frac{\partial \theta}{\partial z} \Big) &= -a\_{\rm v} \theta p \Big[ \frac{\partial h}{\partial t} + \frac{\partial}{\partial x} \Big( h \bar{u}\_{\rm x} \Big) + \frac{\partial}{\partial z} \Big( h \bar{u}\_{z} \Big) \Big] \\ &- \beta\_{0} \Big( \theta - \theta\_{\infty} \Big) + \frac{4\mu}{h} \Big( 3\bar{u}\_{\rm x,P}^{2} + \bar{u}\_{\rm x,C}^{2} + 3\bar{u}\_{z}^{2} \Big) . \end{split} \tag{3.52}$$

Alternatively, when inserting the bulk flow velocities from Eqs. (3.33a) and (3.33b) and additionally exploiting the arc-length transformation *x* = *Rϕ* stated by Eq. (2.6), the temperature equation in Eq. (3.52) can be brought to the more descriptive form

$$\begin{split} \rho \, c\_{\mathrm{V}} \hbar \left[ \frac{\partial \theta}{\partial t} + \frac{\omega\_{0}}{2} \frac{\partial \theta}{\partial \varphi} - \frac{\hbar^{2}}{12\mu} \left( \frac{1}{R^{2}} \frac{\partial p}{\partial \varphi} \frac{\partial \theta}{\partial \varphi} + \frac{\partial p}{\partial z} \frac{\partial \theta}{\partial z} \right) \right] \\ = -a\_{\mathrm{V}} \theta p \left[ \frac{\partial h}{\partial t} + \frac{\omega\_{0}}{2} \frac{\partial h}{\partial \varphi} - \frac{1}{R^{2}} \frac{\partial}{\partial \varphi} \left( \frac{h^{3}}{12\mu} \frac{\partial p}{\partial \varphi} \right) - \frac{\partial}{\partial z} \left( \frac{h^{3}}{12\mu} \frac{\partial p}{\partial z} \right) \right] \\ - \beta\_{0} \left( \theta - \theta\_{\mathrm{\infty}} \right) + R^{2} \omega\_{0}^{2} \frac{\mu}{h} + \frac{h^{3}}{12\mu} \left[ \frac{1}{R^{2}} \left( \frac{\partial p}{\partial \varphi} \right)^{2} + \left( \frac{\partial p}{\partial z} \right)^{2} \right]. \tag{3.53} \end{split}$$

Together with the particular form of the Reynolds equation in Eq. (3.36), the derived temperature equation governs the transient behavior of any compressible lubricant.

#### **3.3.3 Closure of Coupled Equation System**

In addition to the field quantities listed in Eq. (3.37), the PDE in Eq. (3.53) involves

$$\theta = \theta(\varphi, z, t), \qquad \mathsf{c}\_{\mathrm{V}} = \mathsf{c}\_{\mathrm{V}}(\varphi, z, t), \qquad \mathsf{a}\_{\mathrm{V}} = \mathsf{a}\_{\mathrm{V}}(\varphi, z, t), \tag{3.54}$$

all of which depend once more on the time *t* as well as on the coordinates *ϕ* and *z*. Here, to close the temperature equation, exactly the same constitutive laws as for the closure of the Reynolds equation are required, which are a thermal equation of state and a viscous equation of state according to Eqs. (3.2) and (3.7), respectively. While the former already allows for an elimination of the compressibility coefficient field *α*v(*ϕ*, *z*, *t*) owing to Eq. (3.6), an additional caloric equation of state as in Eq. (3.3) must be known to express the heat capacity field *c*v(*ϕ*, *z*, *t*) as Eq. (3.5) stipulates. Moreover, the temperature equation is not only coupled to the Reynolds equation, but also to the other submodels via the insertion of *h*(*ϕ*, *t*), which makes the included coupling quantities the only unknowns besides the fields *ρ*(*ϕ*, *z*, *t*) and *ϑ*(*ϕ*, *z*, *t*).

### **3.4 Statement of Boundary Value Problems**

Remembering that the Reynolds equation in Eq. (3.36) and the temperature equation in Eq. (3.53) are based on the thin-film assumption from Eq. (2.1), neither of them remains valid beyond the axial bearing edges or in the vicinity of foil fixation gaps. The PDEs are thus solved only within certain open and bounded domains *Ωn* ⊂ **R**2 , whereas additional constraints are provided on the corresponding boundaries *∂Ωn*. Altogether, this creates one or more boundary value problems (BVPs) with solutions that belong to spaces of sufficiently many times continuously differentiable functions. At this point, it can already be anticipated that the considered number of pads *N*P determines not only the geometric shapes of *Ωn* and *∂Ωn*, but also the total number of individual BVPs into which the overall fluid problem is subsequently divided.

#### **3.4.1 Classification of Partial Differential Equations**

In the theory on linear PDEs of at most second order, the unambiguous classification into equations of either elliptic, parabolic, or hyperbolic type decides on appropriate boundary and initial conditions that lead to a well-posed problem (see, e.g., [Are+18]). For the nonlinear PDEs of fluid mechanics, on the other hand, such a distinction is much more delicate because the expected characteristics for these equations are often the result of combined elliptic, parabolic, and hyperbolic behavior (see, e.g., [FPS20]).

With regard to the Reynolds equation in Eq. (3.36), it can in particular be observed that the convective terms on the left-hand side are hyperbolic in themselves, whereas the equally important viscous terms on the right-hand side have a parabolic character. The temperature equation in Eq. (3.53) is obviously dominated by convective terms of hyperbolic character, which means that only the few terms that are similar to the right-hand side of the Reynolds equation can be associated with parabolic behavior. In possibly existing regions of a flow where pressure gradients physically disappear, for example due to an isothermal condensation at prescribed vapor pressure [GS19], the parabolic contributions in both of the PDEs may locally even be completely lost.

Besides differences in the choice of suitable solution methods, the crucial point is that well-posed parabolic problems require boundary conditions at all edges, whereas for hyperbolic problems it reveals unphysical to impose restrictions at outflow edges. To enable such a distinction, each boundary is formally partitioned according to

$$
\partial \Omega\_{\mathfrak{n}} = \partial \Omega\_{\mathfrak{n},\mathfrak{i}} \cup \partial \Omega\_{\mathfrak{n},\mathfrak{o}} \cup \partial \Omega\_{\mathfrak{n},\mathfrak{s}} \tag{3.55}
$$

into three pairwise disjoint and not necessarily connected regions, namely an inflow boundary *∂Ωn*,i, an outflow boundary *∂Ωn*,o, and finally a symmetry boundary *∂Ωn*,s. In the following, the assignment of edge points to either the inflow or the outflow boundary is made on the basis of the respective bulk flow velocity components, which implies that the partitioning of the boundaries is a priori unknown.

#### **3.4.2 Identification of Domains and Boundaries**

In the case of a full-arc bearing with *N*P = 0, there is one single domain denoted by

$$
\Omega\_0 = \mathbb{T}^1 \times \left(0, \frac{L}{2}\right),
\tag{3.56}
$$

which corresponds topologically to the lateral surface of a cylinder with height *L*/2. With the definition of the circumferential coordinate *ϕ* on the circle **T**1 = **R**/(2*π***Z**), only 2*π*-periodic solutions are permitted as the schematic sketch in Fig. 3.4a indicates. In the axial direction, it is sufficient to consider only half the bearing length owing to the symmetry that results from identical conditions on both sides of the bearing and from supposing the film thickness in Eq. (2.13) to be independent of the *z*-coordinate. Accordingly, the symmetry boundary of a zero-pad GFB writes

$$
\partial \Omega\_{0,s} = \mathbb{T}^1 \times \{0\} \tag{3.57}
$$

and is directly identified with the circumferential centerline of the bearing at *z* = 0. At the open bearing edge at *z* = *L*/2, the fluid has the possibility to either leave or enter the lubrication gap, which can be judged by whether the axial pressure gradient *∂p*/*∂z* makes the velocity component *u*¯*z* from Eq. (3.33b) positive or negative. Relying on this distinctive criterion, the mathematical expressions

$$\partial \Omega\_{0, \mathbf{i}} = \left\{ (\boldsymbol{\varrho}, \boldsymbol{z}) \in \mathbb{T}^1 \times \left\{ \frac{L}{2} \right\} \; \middle| \; \frac{\partial p}{\partial \boldsymbol{z}} \ge \mathbf{0} \right\}, \tag{3.58a}$$

$$\partial \Omega\_{0,o} = \left\{ (\rho, z) \in \mathbb{T}^1 \times \left\{ \frac{L}{2} \right\} \, \middle| \, \frac{\partial p}{\partial z} < 0 \right\} \tag{3.58b}$$

describe, respectively, the inflow and outflow boundary of a full-arc bearing without of course any possibility of fluid flowing through the non-existent foil fixation gaps.

Coming to the far more important case of a multi-pad bearing with *N*P ≥ 1 pads, the situation is complicated by the fact that there are now *N*P unconnected domains

$$
\Omega\_{\rm n} = (\chi\_{\rm n}, \chi\_{\rm n} + \Delta \chi) \times \left(0, \frac{L}{2}\right), \quad n = 1, \dots, N\_{\rm P} \tag{3.59}
$$

which are separated by the *N*P foil fixation gaps as the sketch in Fig. 3.4b shows. Based on Eqs. (2.2) and (2.3), the circumferential coordinate *ϕ* for the *n*-th domain must range within the leading gap position *χn* and the trailing gap position5 *χn*+*∆χ*. Given that the abstracted fixation gaps themselves have virtually no spatial extension,

5 The notation *χn*+1 is avoided because the *N*P-th gap is followed by the first and not by an (*N*P+1)-th gap.

**(a)** Full-arc bearing with *N*P = 0 pads. **(b)** Multi-pad bearing with *N*P = 3 pads.

**Figure 3.4:** Sketches of exemplary domains on which the relevant fluid PDEs can be solved by referring to boundary conditions for inflow (solid lines), outflow (dotted lines), and symmetry (dot-dashed lines).

the trailing edge of each domain coincides with the leading edge of the next domain, where the *N*P-th pad is followed again by the first pad. When splitting the circle **T**1 in Eqs. (3.57), (3.58a), and (3.58b) accordingly, this leads to *N*P symmetry boundaries

$$
\partial \Omega\_{\rm n,s} = (\chi\_{\rm n} \chi\_{\rm n} + \Delta \chi) \times \{0\}, \quad \textit{n} = \mathbf{1}, \ldots, \textit{N}\_{\rm P} \tag{3.60}
$$

for any multi-pad GFB and, analogously, to as many inflow and outflow boundaries

$$\partial \Omega\_{\mathsf{n},\mathsf{i}} = \left\{ (\mathsf{q}, z) \in \left( \chi\_{\mathsf{n}}, \chi\_{\mathsf{n}} + \Delta \chi \right) \times \left\{ \frac{L}{2} \right\} \, \middle| \, \frac{\partial p}{\partial z} \ge 0 \right\} \cup \left\{ \chi\_{\mathsf{n}} \right\} \times \left( 0, \frac{L}{2} \right), \tag{3.61a}$$

$$\partial \Omega\_{\mathsf{H},0} = \left\{ (\mathsf{q}, z) \in \left( \chi\_{\mathsf{H}}, \chi\_{\mathsf{H}} + \Delta \chi \right) \times \left\{ \frac{L}{2} \right\} \; \middle| \; \frac{\partial p}{\partial z} < 0 \right\} \cup \left\{ \chi\_{\mathsf{H}} + \Delta \chi \right\} \times \left( 0, \frac{L}{2} \right). \tag{3.61b}$$

In Eqs. (3.61a) and (3.61b), the additional inflow regions at *ϕ* = *χn* and the additional outflow regions at *ϕ* = *χn*+*∆χ* can be included without any case distinction since the velocity component *u*¯*x* in Eq. (3.33a) is always positive6 for sufficiently large *ω*0.

#### **3.4.3 Selection of Boundary Conditions**

On the inflow boundaries, the thermodynamic state of the fluid is supposed to be entirely known, which dictates the ambient density *ρ*∞, the ambient temperature *ϑ*∞, and all other properties owing to thermal, caloric, and viscous equations of state. On the outflow boundaries, which are only relevant for parabolic terms, the ambient

6 This supposes a positive rotational speed, but it is trivial to simply swap the edges in the contrary case.

pressure *p*∞ alone acts back on the flow in the upstream direction, while density, temperature, and all other properties are extrapolated from the domain as required. Besides these Dirichlet boundary conditions, homogeneous Neumann boundary conditions are introduced on the symmetry edges and make disappear the gradients of all properties across the bearing centerline. Unsurprisingly, this entire approach shows many parallels to the usual procedure for dealing with general compressible flow problems (see, e.g., [FPS20]), which is transferred here to fluid film lubrication.

While the specification of constant Dirichlet boundary conditions is common practice at the axial bearing edges (see, e.g., [Sze10]), the validity of such an assumption may well be questioned considering the foil fixation gaps. This uncertainty is due to the possible mixing of fluid from the environment with the partially recirculating fluid in the lubrication gap, which would require an empirical mixing parameter that is difficult to determine and strongly related to the exact fixation gap geometry [Gar12]. However, knowing that a good thermal design should ensure the most complete fluid exchange possible according to San Andrés et al. [SK10], such ambient conditions at the fixation gaps are easily justifiable when thinking of idealized bearings anyway.

Combining all of the previous assumptions and referring to domains and boundaries that are defined in Eqs. (3.56), (3.57), (3.58a), (3.58b), (3.59), (3.60), (3.61a), and (3.61b), the fluid BVPs to be solved can ultimately be stated for any number of pads as

$$\begin{cases} \text{Eq. (3.59),} \\ \text{Eq. (3.53),} & \text{in } \Omega\_{\text{n}} \\ \text{Eqs. of state} \\ \rho = \rho\_{\text{co}} \\ \theta = \theta\_{\text{co}} & \text{on } \partial\Omega\_{\text{n},\text{i}} \\ \text{Eq. of state} \\ p = p\_{\text{co}} & \text{on } \partial\Omega\_{\text{n},\text{o}} \\ \partial/\partial z = 0 & \text{on } \partial\Omega\_{\text{n},\text{o}} \end{cases} \quad \text{=BVP}\_{\text{n}\prime} \quad n = \begin{cases} 0 & \text{if } N\_{\text{P}} = 0, \\ 1, \dots, N\_{\text{P}} & \text{if } N\_{\text{P}} \ge 1. \end{cases} \tag{3.62}$$

Here, it should of course not be forgotten that the radial boundary conditions from Eqs. (3.23a), (3.23b), (3.29), and (3.47) are already incorporated into the respective PDEs as an outcome of them being integrated from one bounding wall to the other. From Eq. (3.62), the pressure distribution along the entire circumference as required for the calculation of bearing forces in Eq. (2.25) is obtained by sequencing all the domain-specific solutions together with the values of Dirichlet boundary conditions. To finally complete the problem formulation, suitable initial conditions are specified for the density and temperature fields, so that all other properties can be deduced via the respective equations of state. In subsequent investigations, this implies that the ambient state from *∂Ωn*,i in Eq. (3.62) initially prevails throughout the domains, which should of course not be mistaken for the description of an actual GFB start-up.

# **4 Efficient Foil Models with Dry Friction Characteristics**

As the most obvious feature, the compliant metal foil structure gives such flexibility to the outer fluid film boundary that GFBs can compensate manufacturing tolerances and may adapt to different operating conditions much better than rigid gas bearings. Under fluctuating pressure load, the foil strips are subject to minute relative motion accompanied by dry friction, thus providing an adaptive damping mechanism that reduces or even prevents undesirable vibrations of the entire bearing – rotor system. To reproduce the relevant nonlinear effects while staying computationally efficient, this chapter first reduces the number of parameters to a minimum before proposing a lumped-element bump foil model of reasonable complexity. Based on this approach, it is feasible to describe friction in various ways reaching from a regularization of Coulomb's law of friction to dynamic bristle models capturing stick –slip transitions.

## **4.1 Lumped-Element Modeling Approach**

Besides the lubricating fluid, the foil structure represents the main difficulty but also the main interest in modeling GFBs, which is mostly due to its particular design with a multitude of resulting effects. Since there is consequently no lack of modeling attempts in the literature, many of which seem overly complex and rather inefficient, it is all the more important to first clarify the relevant requirements. In this thesis, one of the objectives consists in a systematic understanding of how friction affects and possibly improves the overall system dynamics. Therefore, the primary purpose of the established minimal model is to qualitatively reproduce this strongly nonlinear influence while introducing the smallest possible number of empirical parameters. In doing so, given that friction is directly related to the occurrence of relative motion, it becomes of course essential to take into account not only local displacements of the compliant structure but also the far-reaching interaction between deforming bumps.

The schematic sketch in Fig. 4.1 shows the considered lumped-element model, which approximates the basic geometry and thus also the kinematics of a bump foil strip by

**Figure 4.1:** Sketch of the lumped-element foil structure model for an exemplary first pad with five bumps, in which the bump foil (in light green) is reduced to an arrangement of rods, springs, and mass particles, whereas the top foil (in dark green) is imitated by simply interpolating linearly between adjacent bumps.

representing each bump as a triangular arrangement of two rigid rods and a spring. As already explained in Section 2.2 in a more general context, this structure model is not only imagined "unwrapped" in conformity with the fluid model, but it is also a plane model that is sufficient for an axially almost constant deformation field. Referring to dimensions that are frequently found in the literature, the undeformed bump length1 2*b*B ≤ *R∆ψ* and the undeformed bump height *h*B define the geometry, which yields the rod length √ (*b* 2 B + *h* 2 B ) by consideration of an undeformed bump2 . Idealized as roller supports in the sketch, the leading edge and the connecting rods between the triangles can move horizontally but cannot detach from the surface, whereas the trailing edge rests on a pinned support in imitation of the foil fixation3 . In this greatest possible simplicity, which nevertheless allows to include friction and bump –bump interaction, the described model is similar to that of Feng et al. [FK10].

The vertical forces acting on the bump apices are supposed to represent the resulting pressure forces that, in reality, are transmitted via areas in contact with the top foil. Depending on the local differences between the pressure in the lubrication gap and the assumed ambient pressure between and below the foils, the top foil adapts to the deformed shape of the bumps and thereby sags into the clear space between them. In a preliminary work by the author [Lei15], the aforementioned foil sagging effect

1 It may surprise in view of Eq. (2.4) that the undeformed bump length is not simply equated with *R∆ψ*. However, as can be seen in the course of the nondimensionalization in Section 5.1, the new quantity 2*b*B is actually introduced on the radial length scale *C* rather than on the circumferential length scale *R*. Hence, both lengths are independent but should fulfill the given compatibility condition for consistency.

2 Note that the fourth bump from the left in Fig. 4.1 is actually shown undeformed because of *x*1,4 = *x*1,5.

3 It goes without saying that the model discussed here is not suitable for describing full bearings without any fixation gap, which makes the limiting zero-pad case *N*P = 0 only applicable with a rigid structure.

is thoroughly investigated by modeling the top foil as a segmented beam structure. The found results show that the wavy shape of the deformed top foil is reflected in a waviness of the fluid pressure field even though this has little or no influence on the qualitative rotordynamic behavior. When discretizing the lubrication gap, however, the local character of this effect requires an exaggeratedly fine computational grid, which seems hardly justifiable compared to the lumped-element bump foil model. In the present thesis, the top foil segments are imagined as being extensible rods that straightly connect bump apices and thereby represent the outer fluid film boundary. For the sake of completeness, Appendix A shows the effect of top foil deformation on the results from Section 6.3.1 and thus builds a bridge between this work and [Lei15].

From an engineering perspective, lumped-element modeling of continuous structures is only permissible if no natural frequencies near the expected excitation frequencies are lost (see, e.g., [DF20]). According to numerical and experimental results reported by Bin Hassan et al. [BB17], the first natural frequency of a typical bump foil is found in the region of 2 kHz, that of an isolated single bump is even about ten times higher. Given that the frequency range of interest correlates with the rotor speed and thus is comparatively low, the foil mass seems to be completely negligible at first glance. However, practical experience suggests that uncertainties involved in the treatment of DAEs are far more delicate than challenges arising from solving systems of stiff but ordinary differential equations [OLS17]. Hence, the roughly approximated mass of individual bumps may be concentrated in discrete mass particles with

$$
\sigma\_{\mathsf{B}} = \rho\_{\mathsf{B}} \pi b\_{\mathsf{B}} \sigma\_{\mathsf{B}} L. \tag{4.1}
$$

The exact value of Eq. (4.1) being insignificant, the foil density *ρ*B is multiplied by the bump-like volume of a half-cylinder shell of radius *b*B, thickness *σ*B, and height *L*, where *σ*B is interpreted as the foil thickness and *L* denotes the bearing length.

The search for an appropriate definition of the spring stiffness is also limited to a rough estimation of an order of magnitude, thus providing a range for parameter variations without restricting the analysis to a particular foil design and geometry. Based on the pioneering approach by Walowit et al. [WA75], which is still widely used despite its surprising simplicity and the many suggestions for improvement, each bump is modeled as an elastic circular arc under axial plane-strain conditions. Depending on Young's modulus *E*B and on Poisson's ratio *ν*B, the stiffness in relation to a force applied at the apex is then given in first approximation by the expression

$$k\_{\rm B} = \frac{E\_{\rm B}L}{2 - 2\nu\_{\rm B}^2} \left(\frac{\sigma\_{\rm B}}{b\_{\rm B}}\right)^3. \tag{4.2}$$

In fact, the simplification behind Eq. (4.2) is essentially based on the finding that the curvature of bumps can be neglected without introducing substantial inaccuracies.

**Table 4.1:** Properties of a typical foil structure.

The inclusion of nonlinear bump stiffness formulations as proposed by Baum [Bau18] is straightforward and may improve quantitative results for practical applications.

Concerning realistic dimensions of the bearing and of the foil structure, the report by Ruscitto et al. [RMG78] is among the first to give an overview of such properties. Ever since, this data is of regular use in the literature4 and is understood as reference for the comparison of results, so that no deviation from this practice is intended here. The report mentions the nickel –chromium INCONEL alloy X-750 as the foil material, which finally yields Table 4.1 when referring to the manufacturer data sheet [Spe04]. Using these parameter values, the bump mass *m*B equals a few tenths of a gram and the bump stiffness *k*B is around 1 N/µm, which is both in agreement with [BB17]. Lastly, the rod length is only slightly larger than *b*B for the small bump height *h*B, meaning that this geometry parameter is expected in the region of a few millimeters.

### **4.2 Foil Strip Kinematics and Interpolation**

Seeing that the lumped-element model of the *n*-th pad5 has *N*B degrees of freedom, the configuration of a given foil strip can be described by *N*B generalized coordinates. Here, the circumferential displacements6 *xn*,*m* of the mass particles shown in Fig. 4.1 seem the most natural choice for a minimum number of independent coordinates. Compared to other selections, each generalized coordinate *xn*,*m* with the resulting generalized velocity *x*˙*n*,*m* later allows for an immediate characterization of friction between the bump foil region that precedes the *m*-th bump and the bearing sleeve.

4 The original dimensions being reported in inches, values in SI units differ slightly between publications. In the present thesis, the converted properties listed in Table 4.1 are also rounded as deemed reasonable.

5 Since all pads are modeled identically, it is sufficient to consider an arbitrary *n*-th pad in isolation. Consequently, the index *n* is always mentioned hereafter but it remains unchanged until further notice.

6 Note that *xn*,*m* is introduced in the negative **e***x*-direction, which later reverses the signs of friction forces.

Moreover, even though not explicitly needed for the described modeling approach, the circumferential displacements and velocities of bump apices may be obtained by just averaging the two concerned generalized coordinates or velocities, respectively.

As indicated by the sketch in Fig. 4.1, the radial displacements *yn*,*m* of the bump apices are expressed relative to their respective positions when all springs are relaxed. Using the Pythagorean theorem to determine the deformed height of the *m*-th bump, the apex displacement depends on generalized coordinates by a kinematic relation

$$y\_{n,m} = \frac{1}{2}\sqrt{4\left(b\_{\mathsf{B}}^2 + h\_{\mathsf{B}}^2\right) - \left(2b\_{\mathsf{B}} + \mathbf{x}\_{n,m} - \mathbf{x}\_{n,m+1}\right)^2} - h\_{\mathsf{B}}.\tag{4.3}$$

To avoid case distinctions in Eq. (4.3) and other equations to come, the symbol *xn*,*N*B+1 is an alias for the vanishing displacement at the foil fixation and must be set to zero. Requiring a mapping at the differential level, it reveals useful to define the Jacobian

$$\mathbf{J}\_{n} = \begin{bmatrix} \frac{\partial y\_{n,1}}{\partial \mathbf{x}\_{n,1}} & \cdots & \frac{\partial y\_{n,1}}{\partial \mathbf{x}\_{n,N\_{\mathrm{B}}}} \\ \vdots & \ddots & \vdots \\ \frac{\partial y\_{n,N\_{\mathrm{B}}}}{\partial \mathbf{x}\_{n,1}} & \cdots & \frac{\partial y\_{n,N\_{\mathrm{B}}}}{\partial \mathbf{x}\_{n,N\_{\mathrm{B}}}} \end{bmatrix} = \begin{bmatrix} I\_{n,1} & -I\_{n,1} & & & 0 \\ & I\_{n,2} & -I\_{n,2} & & \\ & & I\_{n,3} & -I\_{n,3} \\ & & & \ddots & \ddots \\ 0 & & & & I\_{n,N\_{\mathrm{B}}} \end{bmatrix},\tag{4.4}$$

which corresponds to an upper bidiagonal matrix because Eq. (4.3) involves only the generalized coordinates *xn*,*m* and *xn*,*m*+1 of the two directly adjacent mass particles. Moreover, given that Eq. (4.3) is anti-symmetric in the two generalized coordinates, the matrix in Eq. (4.4) is entirely determined in terms of only *N*B diagonal entries

$$J\_{\rm n,m} = \frac{\partial y\_{\rm n,m}}{\partial \mathbf{x}\_{\rm n,m}} = -\frac{1}{2} \frac{2b\_{\rm B} + \mathbf{x}\_{\rm n,m} - \mathbf{x}\_{\rm n,m+1}}{\sqrt{4\left(b\_{\rm B}^2 + h\_{\rm B}^2\right) - \left(2b\_{\rm B} + \mathbf{x}\_{\rm n,m} - \mathbf{x}\_{\rm n,m+1}\right)^2}}.\tag{4.5}$$

Taking finally advantage of Eq. (4.4) when analyzing foil motion, the transformation

$$\begin{bmatrix} \dot{\mathcal{Y}}\_{n,1} \\ \dot{\mathcal{Y}}\_{n,2} \\ \dot{\mathcal{Y}}\_{n,3} \\ \vdots \\ \dot{\mathcal{Y}}\_{n,N\_{\rm B}} \end{bmatrix} = \mathbf{J}\_{n} \begin{bmatrix} \dot{\mathcal{X}}\_{n,1} \\ \dot{\mathcal{X}}\_{n,2} \\ \dot{\mathcal{X}}\_{n,3} \\ \vdots \\ \dot{\mathcal{X}}\_{n,N\_{\rm B}} \end{bmatrix} \tag{4.6}$$

yields all the radial velocities *y*˙*n*,*m* of the bump apices, which depend not only on the generalized velocities, but through Eq. (4.5) also on the coordinates themselves.

Geometrically speaking, the top foil corresponds to straight lines that permanently connect the bump apices, so that position and velocity of any point can be obtained

through linear interpolation based on *yn*,*m* and *y*˙*n*,*m* of the two concerned bumps. At fixation gaps between neighboring pad edges, it seems fair to extrapolate values of the closer bump apex, thus creating a regularizable discontinuity in the middle. With Eqs. (4.3) and (4.6), this makes the deformation and deformation velocity fields

$$q(\boldsymbol{\varrho},t) = q(\dots,\boldsymbol{x}\_{\boldsymbol{n},\boldsymbol{m}\_{\boldsymbol{\nu}}}\dots),\qquad\qquad\frac{\partial}{\partial t}q(\boldsymbol{\varrho},t) = \frac{\mathbf{d}}{\mathbf{d}t}q(\dots,\boldsymbol{x}\_{\boldsymbol{n},\boldsymbol{m}\_{\boldsymbol{\nu}}}\dots)\tag{4.7}$$

functions of the generalized coordinates and velocities, the piecewise definitions of which are not explicitly stated above because they are of purely technical nature7 .

### **4.3 Derivation of Equations of Motion**

The established foil structure model corresponds to *N*P systems of *N*B mass particles, for which *N*P systems of *N*B coupled equations of motion need to be formulated. With the numerous constraints and the repetitive geometry, it seems advantageous to proceed for the *n*-th pad via the Euler–Lagrange equations (see, e.g., [See18])

$$\frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}\_n}{\partial \dot{\boldsymbol{\chi}}\_{n,m}} \right) - \left( \frac{\partial \mathcal{L}\_n}{\partial \boldsymbol{\chi}\_{n,m}} \right) = \mathcal{Q}\_{n,m} \quad m = 1, \ldots, N\_{\mathrm{B}} \tag{4.8}$$

which involve the Lagrangian L*n* and also a number *N*B of generalized forces Q*n*,*m*. The kinetic and the potential energies for the *n*-th pad being particularly simple, it comes

$$\mathcal{L}\_{\rm n} = \frac{m\_{\rm B}}{2} \sum\_{m=1}^{N\_{\rm B}} \dot{\boldsymbol{x}}\_{n,m}^{2} - \frac{k\_{\rm B}}{2} \sum\_{m=1}^{N\_{\rm B}} \left(\boldsymbol{x}\_{n,m} - \boldsymbol{x}\_{n,m+1}\right)^{2} \tag{4.9}$$

an expression in which gravity does not appear for the "unwrapped" structure model. Moreover, formulating the virtual work of the non-conservative fluid forces *Fn*,*m* and friction forces *fn*,*m* in generalized coordinates, a comparison of coefficients yields the generalized forces

$$\begin{bmatrix} \mathcal{Q}\_{n,1} \\ \mathcal{Q}\_{n,2} \\ \mathcal{Q}\_{n,3} \\ \vdots \\ \mathcal{Q}\_{n,N\_{\rm B}} \end{bmatrix} = \mathbf{J}\_{n}^{\rm T} \begin{bmatrix} F\_{n,1} \\ F\_{n,2} \\ F\_{n,3} \\ \vdots \\ F\_{n,N\_{\rm B}} \end{bmatrix} + \begin{bmatrix} f\_{n,1} \\ f\_{n,2} \\ f\_{n,3} \\ \vdots \\ f\_{n,N\_{\rm B}} \end{bmatrix}. \tag{4.10}$$

7 The interested reader is provided with detailed formulae of such a linear top foil interpolation in [Bar16], which represents a preliminary work inspired and co-supervised by the author of the present thesis.

Most notably, Eq. (4.10) takes advantage of the transposed Jacobian **J** T *n* from Eq. (4.4) to relate pressure forces that act on the bump apices to the generalized coordinates.

Inserting Eqs. (4.9) and (4.10) into Eq. (4.8) and calculating the respective derivatives, the obtained outcome for the *n*-th pad is the system of coupled equations of motion

$$
\begin{aligned}
\begin{bmatrix} \mathbb{X}\_{n,1} \\ \mathbb{X}\_{n,2} \\ \mathbb{X}\_{n,3} \\ \vdots \\ \mathbb{X}\_{n,N\_{\mathrm{B}}} \end{bmatrix} + \omega\_{\mathrm{B}}^{2} \begin{bmatrix} 1 & -1 & & & 0 \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & \ddots & \ddots & \ddots \\ 0 & & & -1 & 2 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_{n,1} \\ \mathbf{x}\_{n,2} \\ \mathbf{x}\_{n,3} \\ \vdots \\ \mathbf{x}\_{n,N\_{\mathrm{B}}} \end{bmatrix} \\ &= -\frac{1}{m\_{\mathrm{B}}} \begin{bmatrix} \begin{bmatrix} \boldsymbol{\mu}\_{n,1} & -\boldsymbol{\boldsymbol{\mu}}\_{n,1} & & & \boldsymbol{0} \\ & \boldsymbol{\mu}\_{n,2} & -\boldsymbol{\boldsymbol{\mu}}\_{n,2} & & \\ & & \boldsymbol{\mu}\_{n,3} & -\boldsymbol{\mu}\_{n,3} \\ & & \ddots & \ddots & \ddots \\ & & & \ddots & \ddots \\ 0 & & & & \boldsymbol{\mu}\_{n,N\_{\mathrm{B}}} \end{bmatrix} \end{bmatrix} \begin{bmatrix} \boldsymbol{\mu}\_{n,1} \\ \boldsymbol{\mu}\_{n,2} \\ \boldsymbol{\mu}\_{n,3} \\ \vdots \\ \boldsymbol{\mu}\_{n,N\_{\mathrm{B}}} \end{bmatrix} - \frac{1}{m\_{\mathrm{B}}} \begin{bmatrix} \boldsymbol{\mu}\_{n,1} \\ \boldsymbol{\mu}\_{n,2} \\ \boldsymbol{\mu}\_{n,3} \\ \vdots \\ \boldsymbol{\mu}\_{n,N\_{\mathrm{B}}} \end{bmatrix}, \end{aligned} \tag{4.110}
$$

in which the inherent dynamics are characterized by the natural angular frequency

$$
\omega\_{\mathbf{B}} = \sqrt{\frac{k\_{\mathbf{B}}}{m\_{\mathbf{B}}}}.\tag{4.12}
$$

While the linear left-hand side of Eq. (4.11) represents a harmonic oscillator chain, the right-hand side corresponds to nonlinear fluid forces and frictional damping.

To assess the resulting fluid forces that are exerted on the respective bump apices, the locally load-distributing function of a top foil is roughly imitated by integrating elementary pressure forces over the associated area of negligible curvature. With the ambient pressure *p*∞ underneath the foil structure counteracting the pressure field *p* in the lubrication gap, this becomes

$$F\_{n,m} = 2\int\_0^{\frac{L}{2}} \int\_{\psi\_{n,m}}^{\psi\_{n,m} + \Delta\psi} (p - p\_{\infty}) \operatorname{Rd}\varphi \,\mathrm{d}z. \tag{4.13}$$

Since any circumferential displacements of the foil structure are judged negligible compared to the length scale `0 of overall bearing dimensions, the integration limits in Eq. (4.13) rely on Eqs. (2.4) and (2.5) that refer to the undeformed foil geometry. Eventually, considering the fluid forces from Eq. (4.13) for the solution of Eq. (4.11) and inserting the kinematic relationships from Eq. (4.7) into Eqs. (2.13) and (2.15), the occurring fluid –structure interaction is entirely captured in both directions.

## **4.4 Discussion of Dry Friction Formulations**

As a fundamental design feature of GFBs, which are prone to undesirable vibrations due to the low viscosity of gases, the compliant and movable foil structure has the beneficial property of introducing dry friction as additional damping to the system. For the description of friction-induced effects, which is of course not only a matter in GFBs but also in many other mechanical systems, the available literature offers a very large selection of models in varying levels of complexity for different purposes. A general overview of relevant approaches is given in the thesis of Kapelke [Kap19], while a comprehensive review of many friction models can be found in [Mar+16].

Relative motion and therefore dry friction occur in bump-type GFBs principally in two different types of contact regions, which are located between bearing sleeve and bump foil on the one hand and between bump foil and top foil on the other hand. In connection with the present work, a preliminary finite-element analysis [Alv18] distinguishes exactly between the respective influences of friction at this or that spot, which leads to the simple conclusion that all observed effects are practically the same. Consequently, and also in the very spirit of a lumped-element modeling approach, the only concentrated friction forces *fn*,*m* in Fig. 4.1 act directly on the mass particles. Setting the normal forces equal to the nearby pressure forces *Fn*,*m* from Eq. (4.13), a combined expression for Coulomb's law of dynamic friction (with coefficient *µ*B) and for a corresponding regime of static friction can be given in the set-valued form

$$f\_{n,m} \in \begin{cases} \left\{-\mu\_{\rm B} F\_{n,m} \right\} & \text{if } \pounds\_{n,m} < 0, \\\left[ -\mu\_{\rm B} F\_{n,m} + \mu\_{\rm B} F\_{n,m} \right] & \text{if } \pounds\_{n,m} = 0, \\\left\{ +\mu\_{\rm B} F\_{n,m} \right\} & \text{if } \pounds\_{n,m} > 0. \end{cases} \tag{4.14}$$

In Eq. (4.14) and also in other upcoming friction formulations, the normal force *Fn*,*m* is assumed to be positive, so that the concerned friction force is set to zero otherwise. As a matter of fact, this can be considered a basic representation of foil detachment, knowing that more sophisticated contact mechanics approaches [AB19] are difficult to be integrated into a model that ultimately results in a single dynamical system.

To further complicate the problem, the insertion of Eq. (4.14) into Eq. (4.11) leads to a system of differential inclusions, which can roughly be seen as multi-valued ODEs. For such Filippov-type equations, especially the study of bifurcations is an own field of research that is thoroughly discussed, for instance, in the thesis of Leine [Lei00]. Since the GFB–rotor model in the present thesis reveals far too complicated to be treated using such a method, it is discarded in favor of some regularized formulation. In this context, it is common practice to neglect static friction and to express Eq. (4.14)

by a sign function, which can then be approximated by a hyperbolic tangent function. Introducing a small regularization velocity *u*B *x*˙*n*,*m*, this transforms Eq. (4.14) into

$$f\_{n,m} = \mu\_{\rm B} F\_{n,m} \tanh\left(\frac{\dot{\mathfrak{X}}\_{n,m}}{\mu\_{\rm B}}\right) \tag{4.15}$$

and can thus be included in Eq. (4.11) without affecting the solution requirements. On the downside, it must of course not be forgotten that this convenient formulation has some deficits in view of possible sticking conditions, however, it still represents an important model improvement when compared to equivalent viscous damping.

For the sake of completeness, it must be mentioned that more advanced formulations from the class of bristle models can of course also be considered, including notably the Dahl model [Dah76], LuGre model [Can+95], and elasto-plastic model [Dup+02]. The idea behind such approaches basically consists in the introduction of a tangential contact compliance to account for pre-sliding deformations of the surface asperities, which thereby provides some kind of physically motivated regularization [Kap19]. Applied to the foil structure model, each bump needs an additional coordinate8 *zn*,*m* to describe the linearly elastic deflection (with stiffness *kz*) of the associated bristle. Assuming for brevity a most simple evolution rule for this internal bristle variable following the Dahl model [Dah76], the friction equations coupled to Eq. (4.11) write

$$f\_{\rm n,m} = k\_z z\_{\rm n,m}, \qquad \qquad \qquad \dot{z}\_{\rm n,m} = \dot{x}\_{\rm n,m} - \left| \dot{x}\_{\rm n,m} \right| \frac{f\_{\rm n,m}}{\mu\_{\rm B} F\_{\rm n,m}}. \tag{4.16}$$

Especially when the right-hand side of the ODE in Eq. (4.16) is modified by adding a suitable switch function as in the elasto-plastic model by Dupont et al. [Dup+02], such a friction model reveals capable of actually reproducing stick –slip transitions.

Nevertheless, no model that truly accounts for sticking can circumvent the inherent and physically originated difficulties that inevitably arise when facing dry friction. Even though Eq. (4.16) does not lead to a differential inclusion as does Eq. (4.14), another directly related issue is observed upon consideration of the evolution rule. More precisely, the previously mentioned switch function in the elasto-plastic model must be defined such that Eq. (4.16) becomes *z*˙*n*,*m* = *x*˙*n*,*m* during phases of sticking. In this case, the concerned bristle effectively loses its additional degree of freedom, so that suitable solution strategies have to deal with varying problem dimensions. Since it is not intended to perform stability and bifurcation analyses for such systems within the scope of this thesis, most of the presented findings are based on Eq. (4.15), while some time integration results with a bristle model are also published in [LSB18].

8 The symbol *zn*,*m* is chosen in reference to the literature and has nothing to do with the axial direction.

# **5 Computational Analysis of GFB –Rotor Systems**

The comprehensive fluid–structure–rotor modeling approach that is described in previous chapters leads to a multi-dimensional system of coupled PDEs and ODEs, so that an analytical treatment does not seem to be a very promising undertaking. Consequently, this chapter outlines the way toward a numerical solution strategy, according to which the equations are nondimensionalized and spatially discretized to give a mathematical dynamical system that can be investigated computationally.

### **5.1 Nondimensional Problem Formulation**

Owing to the coupled nature of the problem, it appears indispensable to establish a nondimensional formulation that remains consistent throughout all three submodels. First of all, the introduction of a characteristic time scale *t*0 yet to be determined defines the nondimensional time

$$
\pi = \frac{t}{t\_0},
\tag{5.1}
$$

with the corresponding nondimensional time derivative denoted by (. . .) 0 for short. Coming to the independent Cartesian coordinates that describe the lubrication gap, a suitable nondimensional circumferential coordinate *ϕ* is already given in Eq. (2.6) and is completed, with the bearing length *L*, by a nondimensional axial coordinate

$$Z = \frac{z}{L} \,. \tag{5.2}$$

In the direction of the film thickness, the radial clearance *C* is the natural reference to define nondimensional forms of all concerned kinematic quantities in Chapter 2, most notably foil structure deformations and rotor displacements, according to

$$H = \frac{h}{\mathbb{C}}\,'\qquad Q = \frac{q}{\mathbb{C}}\,'\qquad X = \frac{\tilde{\xi}}{\mathbb{C}}\,'\qquad Y = \frac{\eta}{\mathbb{C}}\,'\qquad \varepsilon = \frac{e}{\mathbb{C}}\,'\qquad \gamma = \Gamma.\tag{5.3}$$

With Eq. (5.3), any dependencies on the dimensional coordinates are also transformed to their nondimensional counterparts, which implies *γ*(*τ*) = *Γ*(*τt*0) for instance.

For the nondimensionalization of fluid BVPs, the relevant thermodynamic properties in Chapter 3 are related to the respective values under normal conditions by setting1

$$P = \frac{p}{p^{\circ}}, \qquad D = \frac{\rho}{\rho^{\circ}}, \qquad T = \frac{\theta}{\theta^{\circ}}, \qquad V = \frac{\mu}{\mu^{\circ}}, \qquad W = \frac{w}{w^{\circ}}.\tag{5.4}$$

Then, it seems a natural choice to express the nondimensional derivative properties

$$A\_{\rm V} = \theta^{\diamond} a\_{\rm V} \tag{5.5} \tag{5.6}$$

using the same reference values as in Eq. (5.4) and combining them appropriately. Even though not stated explicitly for the sake of brevity, such nondimensionalizations apply analogously to fluid boundary conditions such as *T*∞ = *ϑ*∞/*ϑ* ◦ and others. Taking advantage of Eqs. (5.1) to (5.5), the PDE from Eq. (3.36) results immediately in the nondimensional Reynolds equation for compressible fluids

$$\frac{\partial}{\partial \tau} \left( DH \right) + \frac{\Lambda\_0}{2} \frac{\partial}{\partial \varphi} \left( DH \right) = \frac{\partial}{\partial \varphi} \left( \frac{DH^3}{2V} \frac{\partial P}{\partial \varphi} \right) + \kappa\_0^2 \frac{\partial}{\partial Z} \left( \frac{DH^3}{2V} \frac{\partial P}{\partial Z} \right), \tag{5.6}$$

just like the PDE from Eq. (3.53) yields the nondimensional temperature equation

$$\begin{split} D\mathbb{C}\_{\mathrm{V}}H\left[\frac{\partial T}{\partial \boldsymbol{\tau}} + \frac{\Lambda\_{0}}{2}\frac{\partial T}{\partial \boldsymbol{\rho}} - \frac{H^{2}}{2V}\left(\frac{\partial P}{\partial \boldsymbol{\rho}}\frac{\partial T}{\partial \boldsymbol{\rho}} + \kappa\_{0}^{2}\frac{\partial P}{\partial \boldsymbol{Z}}\frac{\partial T}{\partial \boldsymbol{Z}}\right)\right] \\ = -A\_{\mathrm{V}}T\mathbb{P}\left[\frac{\partial H}{\partial \boldsymbol{\tau}} + \frac{\Lambda\_{0}}{2}\frac{\partial H}{\partial \boldsymbol{\rho}} - \frac{\partial}{\partial \boldsymbol{\rho}}\left(\frac{H^{3}}{2V}\frac{\partial P}{\partial \boldsymbol{\rho}}\right) - \kappa\_{0}^{2}\frac{\partial}{\partial \boldsymbol{Z}}\left(\frac{H^{3}}{2V}\frac{\partial P}{\partial \boldsymbol{Z}}\right)\right] \\ - \Theta\_{0}\Big(T - T\_{\infty}\Big) + \frac{\Lambda\_{0}^{2}}{6}\frac{V}{H} + \frac{H^{3}}{2V}\left[\left(\frac{\partial P}{\partial \boldsymbol{\rho}}\right)^{2} + \kappa\_{0}^{2}\left(\frac{\partial P}{\partial \boldsymbol{Z}}\right)^{2}\right]. \end{split} \tag{5.7}$$

Having a closer look at Eqs. (5.6) and (5.7), it becomes clear that their particular form fixes not only the time scale

$$t\_0 = \frac{6\mu^\circ}{p^\circ} \left(\frac{R}{\mathcal{C}}\right)^2,\tag{5.8}$$

but also results in an appearance of the three nondimensional similarity parameters

$$
\epsilon\_0 = \frac{R}{L}\,'\qquad \Lambda\_0 = t\_0\omega\_0 = \frac{6\mu^\circ}{p^\circ} \left(\frac{R}{\mathcal{C}}\right)^2 \omega\_0\qquad \Theta\_0 = \frac{\theta^\circ}{(p^\circ)^2 \mathcal{C}^3} \left(\frac{R}{\mathcal{C}}\right)^2 \beta\_0. \tag{5.9}
$$

1 This supposes the vapor quality to be larger than 0 % under normal conditions, which holds for R-245fa.

In Eq. (5.9), the first parameter *κ*0 is a purely geometric number referring to the ratio between radius and length of a GFB, which distinguishes short from long bearings. The second parameter *Λ*0 resembles the Sommerfeld number (see, e.g., [Sze10]) and is subsequently denoted as bearing number or nondimensional rotational speed. As an important matter of fact, the bearing number is the only place in the overall nondimensional problem in which the angular frequency *ω*0 of the rotor is involved, while in particular the adopted time scale in Eq. (5.8) is independent of the speed. For run-up and coast-down simulations or when considering *Λ*0 the bifurcation parameter in solution continuation, this ensures to not misleadingly bias the results. Finally, the third parameter *Θ*0 in Eq. (5.9) is called thermal management number and describes the general ability of a GFB to transfer heat to the environment.

Concerning the foil structure model, the nondimensional form of displacements should be compatible with the nondimensional deformation field *Q* from Eq. (5.3), so the clearance *C* is used in the circumferential and in the radial direction to give

$$X\_{\rm n,m} = \frac{x\_{\rm n,m}}{\rm C}, \quad Y\_{\rm n,m} = \frac{y\_{\rm n,m}}{\rm C}, \quad Z\_{\rm n,m} = \frac{z\_{\rm n,m}}{\rm C}, \quad B\_{\rm B} = \frac{b\_{\rm B}}{\rm C}, \quad H\_{\rm B} = \frac{h\_{\rm B}}{\rm C}.\tag{5.10}$$

Any loading of the foil model being related to surface forces due to pressure, nondimensional forms of Eqs. (4.13) to (4.16) are given by the respective expressions

$$
\Pi\_{\rm n,m} = \frac{F\_{\rm n,m}}{p^\diamond \mathcal{R}L}, \qquad \qquad \qquad \qquad \pi\_{\rm n,m} = \frac{f\_{\rm n,m}}{p^\diamond \mathcal{R}L}. \tag{5.11}
$$

Altogether, this transforms the system of equations of motion from Eq. (4.11) into2

$$
\begin{bmatrix} \mathbf{X}\_{n,1}^{\prime\prime} \\ \mathbf{X}\_{n,2}^{\prime\prime} \\ \mathbf{X}\_{n,3}^{\prime\prime} \\ \vdots \\ \mathbf{X}\_{n,N\_{\mathrm{B}}}^{\prime\prime} \end{bmatrix} + \Omega\_{\mathrm{B}}^{2} \left[ \dots \right] \begin{bmatrix} \mathbf{X}\_{n,1} \\ \mathbf{X}\_{n,2} \\ \mathbf{X}\_{n,3} \\ \vdots \\ \mathbf{X}\_{n,N\_{\mathrm{B}}} \end{bmatrix} = -\frac{1}{M\_{\mathrm{B}}} \left[ \dots \right] \begin{bmatrix} II\_{n,1} \\ II\_{n,2} \\ II\_{n,3} \\ \vdots \\ II\_{n,N\_{\mathrm{B}}} \end{bmatrix} - \frac{1}{M\_{\mathrm{B}}} \begin{bmatrix} \pi\_{n,1} \\ \pi\_{n,2} \\ \pi\_{n,3} \\ \vdots \\ \pi\_{n,N\_{\mathrm{B}}} \end{bmatrix}, \tag{5.12}
$$

so that the foil behavior is characterized in terms of two nondimensional parameters

$$\Omega\_{\mathsf{B}} = t\_0 \omega\_{\mathsf{B}} = \frac{6\mu^{\diamond}}{p^{\diamond}} \left(\frac{\mathsf{R}}{\mathsf{C}}\right)^2 \omega\_{\mathsf{B}\prime} \qquad \qquad M\_{\mathsf{B}} = \frac{p^{\diamond}}{36(\mu^{\diamond})^2 L} \left(\frac{\mathsf{C}}{\mathsf{R}}\right)^5 m\_{\mathsf{B}\prime} \tag{5.13}$$

referred to as nondimensional bump frequency number and bump mass number. Besides Eq. (5.13), the nondimensional friction coefficient *µ*B in Eqs. (4.14) to (4.16)

2 Since the two large matrices in Eq. (4.11) are formally identical to their nondimensional equivalents, they are not repeated here. In the same way, the large matrix in Eq. (2.23) is later omitted for the rotor.

represents a native nondimensional number that is unaffected by any reformulation. When considering a bristle friction model, however, the ODEs from Eq. (4.16) become

$$
\pi\_{\rm n,m} = \mathcal{K}\_{\rm Z} \mathcal{Z}\_{\rm n,m}, \qquad \qquad \mathcal{Z}'\_{\rm n,m} = \mathcal{X}'\_{\rm n,m} - \left| \mathcal{X}'\_{\rm n,m} \right| \frac{\pi\_{\rm n,m}}{\mu\_{\rm B} \Pi\_{\rm n,m}}, \tag{5.14}
$$

in which the definitions of nondimensional displacements and forces in Eqs. (5.10) and (5.11) make the bristle stiffness number *KZ* = *C*/(*p* ◦*RL*)*kz* another parameter.

Regarding the rotor model, the same nondimensionalization as in Eq. (5.3) leads to

$$X\_{\mathbb{R}} = \frac{\tilde{\xi}\_{\mathbb{R}}}{\mathbb{C}}\,' \, \qquad \qquad \qquad \qquad \qquad \boldsymbol{\Upsilon}\_{\mathbb{R}} = \frac{\eta\_{\mathbb{R}}}{\mathbb{C}}\,' \, \qquad \qquad \qquad \boldsymbol{\epsilon}\_{\mathbb{R}} = \frac{\boldsymbol{e}\_{\mathbb{R}}}{\mathbb{C}} \, \tag{5.15}$$

and the same consideration as in Eq. (5.11) gives the nondimensional bearing forces

$$
\Pi\_X = \frac{F\_\xi}{p^\circ RL} \, , \qquad \qquad \qquad \Pi\_Y = \frac{F\_\eta}{p^\circ RL} \, . \tag{5.16}
$$

Applying this to Eq. (2.23) results in a system of nondimensional equations of motion

$$\begin{aligned} \begin{bmatrix} \mathbf{X}^{\prime\prime} \\ \mathbf{Y}^{\prime\prime} \\ \mathbf{X}^{\prime}\_{\mathbf{R}} \\ \mathbf{Y}^{\prime\prime}\_{\mathbf{R}} \end{bmatrix} + \frac{2D\_{\mathbf{R}}}{1 - \Xi\_{\mathbf{R}}} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{Y}^{\prime}\_{\mathbf{R}} \\ \mathbf{Y}^{\prime}\_{\mathbf{R}} \end{bmatrix} + \frac{\Omega^{2}\_{\mathbf{R}}}{1 - \Xi\_{\mathbf{R}}} \begin{bmatrix} \mathbf{X} \\ \dots \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \\ \mathbf{X}\_{\mathbf{R}} \\ \mathbf{Y}\_{\mathbf{R}} \end{bmatrix} \\ &= \frac{2}{\Xi\_{\mathbf{R}}M\_{\mathbf{R}}} \begin{bmatrix} II\_{X} \\ II\_{Y} \\ 0 \\ 0 \\ \end{bmatrix} - \epsilon\_{\mathbf{R}}\Lambda\_{0}^{2} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \sin(\Lambda\_{0}\boldsymbol{\tau}) \\ \cos(\Lambda\_{0}\boldsymbol{\tau}) \\ \end{bmatrix} - G \begin{bmatrix} \mathbf{0} \\ \mathbf{1} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} \tag{5.17} \end{aligned} (5.17)$$

with four newly appearing nondimensional numbers

$$\Omega\_{\mathbb{R}} = t\_0 \omega\_{\mathbb{R}} = \frac{6\mu^{\diamond}}{p^{\diamond}} \left(\frac{\mathbb{R}}{\mathbb{C}}\right)^2 \omega\_{\mathbb{R}\prime} \qquad \qquad D\_{\mathbb{R}} = t\_0 \delta\_{\mathbb{R}} = \frac{6\mu^{\diamond}}{p^{\diamond}} \left(\frac{\mathbb{R}}{\mathbb{C}}\right)^2 \delta\_{\mathbb{R}\prime} \tag{5.18a}$$

$$M\_{\rm R} = \frac{p^{\diamond}}{36(\mu^{\diamond})^2 L} \left(\frac{\mathbb{C}}{\mathbb{R}}\right)^5 m\_{\rm R} \qquad\qquad\qquad G = \frac{36(\mu^{\diamond})^2}{(p^{\diamond})^2 R} \left(\frac{\mathbb{R}}{\mathbb{C}}\right)^5 \mathfrak{g}.\tag{5.18b}$$

While *Ω*R and *M*R are constructed and interpreted in the same way as *Ω*B and *M*B from Eq. (5.13), the characterization of the rotor requires two additional parameters. To this purpose, the viscous damping is reflected in *D*R according to Eq. (5.18a), while *G* in Eq. (5.18b) represents the gravitational loading of a horizontal rotor. Obviously, the discussed nondimensionalization applies analogously to Eq. (2.24) for the rigid rotor model, which is then independent of the frequency number *Ω*R.

### **5.2 Finite-Difference Discretization Schemes**

While the foil structure and the rotor are represented by discrete models that lead directly to systems of exclusively time-dependent ODEs, the continuous fluid model produces PDEs that depend on both time and space. Lacking analytical solutions, the common remedy in computational fluid dynamics (CFD) consists in using spatial discretization schemes to transform such PDEs into systems of coupled ODEs that deliver approximate solutions. Among the different classes of available approaches, e.g., finite-element (FEM), finite-volume (FVM), or finite-difference (FDM) methods, the selected FDM provides the greatest flexibility when facing non-standard PDEs.

Since the nondimensional domains that describe the lubrication gap are rectangular, the introduction of a corresponding structured computational grid is not complicated when specifying the number of grid points in the respective directions by *Nϕ* and *NZ*. However, referring to the overall fluid region to define the equidistant grid spacings

$$
\Delta \varphi = \frac{2\pi}{N\_{\varphi} - 1} \, \, \, \, \quad \qquad \qquad \Delta Z = \frac{1}{N\_{\Sigma} - 1} \, \, \, \, \} \, \, \, \tag{5.19}
$$

some additional care is indicated to ensure compatibility with the number of pads *N*P, with the number of bumps *N*B, and with the axial symmetry boundary conditions. More specifically, requiring grid lines that coincide exactly with each fixation gap, with each foil –sleeve contact line, and with the bearing's circumferential centerline, a closer look at the exemplary grids in Fig. 5.1 yields two discretization conditions3

$$\frac{\Delta\psi}{\Delta\varrho} = \frac{N\_{\varphi} - 1}{N\_{\text{P}} N\_{\text{B}}} \in \mathbb{N}, \qquad \qquad \frac{\frac{1}{2}}{\Delta Z} = \frac{N\_{\text{Z}} - 1}{2} \in \mathbb{N}.\tag{5.20}$$

Finally, having established the computational grid according to Eqs. (5.19) and (5.20), the grid origin is anchored to the nondimensional physical domain in such a way that the integer grid coordinates (*i*, *j*) correspond to discretized physical coordinates

$$
\varphi\_i = (i - 1)\Delta\varphi + \chi\_{1\prime} \quad i = 1, \ldots, N\_{\Phi\prime} \tag{5.21a}
$$

$$Z\_{j} = (j - 1)\Delta Z - \frac{1}{2}\,'\,\,\,j = 1\,\,\dots\,\,N\_{\mathbf{Z}}.\tag{5.21b}$$

Inserting Eq. (5.21a) into the nondimensional forms of Eqs. (2.13), (2.15), and (4.7), these field quantities transform into a discretized fluid film thickness *Hi* and into a discretized foil structure deformation *Qi* with discretized time derivatives *H*0 *i* and *Q*0 *i* .

3 Graphically, this implies that *∆ψ* is an integer multiple of *∆ϕ*, while 1/2 is an integer multiple of *∆Z*.

**(a)** Full-arc bearing with *N*P = 0 pads (described by the total number of *N*FDM = 27 internal nodes).

**(b)** Multi-pad bearing with *N*P = 3 pads (described by the total number of *N*FDM = 18 internal nodes).

**Figure 5.1:** Sketches of exemplary computational grids for the FDM with *Nϕ* × *NZ* = 10 × 7 grid points, which are either internal nodes (circles), boundary nodes (squares), or mirror nodes (dashed circles).

Considering once more the sketched grids in Fig. 5.1, it becomes obvious that only

$$N\_{\rm FDM} = \frac{1}{2} \left( N\_{\rm \varphi} - 1 - N\_{\rm P} \right) \left( N\_{\rm Z} - 1 \right) \tag{5.22}$$

grid points are internal nodes that lie in domains on which PDEs need to be solved. At each internal node, a discretized density *Di*,*j* and a discretized temperature *Ti*,*j* are introduced as variables, on which all other discretized thermodynamic properties such as the pressure *Pi*,*j* or the viscosity *Vi*,*j* depend via the given equations of state. Around the internal nodes, the remaining grid points are either boundary nodes, which store constant values given by Dirichlet boundary conditions, or mirror nodes, in which values from nodes of the other symmetric half of the domain are reflected.

The basic idea behind any FDM consists in the replacement of partial derivatives by corresponding difference quotients, which is close to the definition of a derivative and can be justified by reformulating a truncated Taylor series (see, e.g., [FPS20]). Evaluated on the introduced computational grid, this makes any spatial derivative at the position of a given internal node (*i*, *j*) an algebraic function of discrete values obtained from the node itself and from four neighboring nodes (*i* ± 1, *j*) and (*i*, *j* ± 1). Regarding the PDEs in Eqs. (5.6) and (5.7), the diverse nature of occurring derivatives as discussed in Section 3.4.1 leads to a distinction between terms of parabolic and terms of hyperbolic character, which affects the choice of suitable approximations. Following an approach that is typically used in CFD when dealing with any kind of diffusion equation, the second-order terms with inner and outer derivatives on the

right-hand side of Eq. (5.6) can be expressed by central differences (see, e.g., [FPS20])

$$\left. \frac{\partial}{\partial \boldsymbol{\varrho}} \left( \frac{D\boldsymbol{H}^{3}}{2\boldsymbol{V}} \frac{\partial \boldsymbol{P}}{\partial \boldsymbol{\varrho}} \right) \right|\_{i,j} \approx \left. \frac{\frac{D\_{i+\frac{\mathsf{I}\_{2}}{2},\mathsf{I}\_{2}^{\mathsf{I}\_{2}}}{2\mathsf{I}\_{i}+\mathsf{I}\_{2},j} \left(\mathsf{P}\_{i+1,j}-\mathsf{P}\_{i,j}\right) - \frac{D\_{i-\frac{\mathsf{I}\_{2}},\mathsf{I}\_{2}^{\mathsf{I}\_{2}}}{2\mathsf{I}\_{i}-\mathsf{I}\_{2},j} \left(\mathsf{P}\_{i,j}-\mathsf{P}\_{i-1,j}\right)}{\Delta \boldsymbol{\varrho}^{2}} \right|\_{i} \quad \text{(5.23a)}$$

$$\left. \frac{\partial}{\partial Z} \left( \frac{DH^3}{2V} \frac{\partial P}{\partial Z} \right) \right|\_{i,j} \approx \left. \frac{\frac{D\_{i,j+\rangle\_2}{2V} H\_i^3}{2V\_{i,j+\/\hbar\_2}} \left( P\_{i,j+1} - P\_{i,j} \right) - \frac{D\_{i,j-\/\hbar\_2} H\_i^3}{2V\_{i,j-\/\hbar\_2}} \left( P\_{i,j} - P\_{i,j-1} \right)}{\Delta Z^2}. \tag{5.25b}$$

Such schemes, in which the fractional indices *i* ± 1/2 or *j* ± 1/2 denote mean values of adjacent nodes, are also used for the terms in Eq. (5.7) without density as a factor. For the first-order derivatives in Eqs. (5.6) and (5.7), from which the derivatives of temperature are picked out as an illustrative example, numerical stability becomes a relevant issue that can be addressed by using upwind differences (see, e.g., [FPS20])

$$
\left| \frac{\partial T}{\partial \boldsymbol{\varphi}} \right|\_{\boldsymbol{i},\boldsymbol{j}} \approx \frac{T\_{\boldsymbol{i},\boldsymbol{j}} - T\_{\boldsymbol{i}-1,\boldsymbol{j}}}{\Delta \boldsymbol{\varphi}}, \qquad \qquad \left| \frac{\partial T}{\partial \boldsymbol{Z}} \right|\_{\boldsymbol{i},\boldsymbol{j}} \approx \begin{cases} \frac{T\_{\boldsymbol{i},\boldsymbol{j}+1} - T\_{\boldsymbol{i},\boldsymbol{j}}}{\Delta \boldsymbol{Z}} & \text{if } P\_{\boldsymbol{i},\boldsymbol{j}+1} \geq P\_{\boldsymbol{i},\boldsymbol{j}-1}, \\\frac{T\_{\boldsymbol{i},\boldsymbol{j}} - T\_{\boldsymbol{i},\boldsymbol{j}-1}}{\Delta \boldsymbol{Z}} & \text{otherwise.} \end{cases} \tag{5.24}
$$

In order to determine the direction of flow in Eq. (5.24), the same consideration as for the distinction between inflow and outflow edges in Eqs. (3.58) and (3.61) is applied. Accordingly, the circumferential velocity is assumed to be always positive, whereas the axial velocity must be judged by means of the discretized axial pressure gradient. Of course, the first-order derivatives in Eqs. (5.6) and (5.7) that involve *P*, *D*, or *H* are discretized in analogy to Eq. (5.24) even though not stated for the sake of brevity.

When the finite-difference stencil moves around to each of the *N*FDM internal nodes, the FDM-discretized Reynolds equation and temperature equation produce together a total number of 2*N*FDM first-order ODEs, which are reordered to the abstract forms

$$D'\_{i,j} = \text{FDM}\_{\mathcal{D}}\left\{ D\_{i,j}, D\_{i\pm 1,j}, D\_{i,j\pm 1\prime}, T\_{i,j\prime}, T\_{i\pm 1,j\prime}, T\_{i,j\pm 1\prime}, H\_{i\prime}, H\_{i\pm 1\prime}, H'\_{i\prime}, H'\_{i\pm 1} \right\},\tag{5.25a}$$

$$T'\_{i,j} = \text{FDM}\_{\Gamma} \left\{ D\_{i,j}, D\_{i \pm 1, j}, D\_{i, j \pm 1 \prime}, T\_{i,j}, T\_{i \pm 1, j \prime}, T\_{i, j \pm 1 \prime}, H\_{i \prime}, H\_{i \pm 1 \prime}, H'\_{i \prime}, H'\_{i \pm 1} \right\}. \tag{5.25b}$$

Realizing that peripheral nodes of the stencil may fall on boundary or mirror nodes, it becomes clear that these equations represent not only PDEs but even entire BVPs. Most remarkably, it is finally the very nature of the selected schemes that enables the central-difference parabolic terms to consider boundary conditions at all edges, while the upwind-difference hyperbolic terms do not even "see" outflow boundaries. Consequently, any solutions of the discretized ODEs in Eq. (5.25) are also considered approximate solutions of the original BVPs from Eq. (3.62) in nondimensional form. In order to assess the quality of the approximation, grid independence studies as shown exemplarily in Appendix B are conducted prior to performing actual analyses.

### **5.3 Nonlinear Dynamical System Representation**

Although the establishment of a dynamical system representation is a decisive step, the feasibility of which is an important merit of the present thesis, this operation does not cause major difficulties in view of the preparatory work already accomplished. It only remains to identify suitable variables that uniquely describe the system state when forming a smallest possible fluid –structure – rotor state vector (or phase vector)

$$\mathbf{S} = \left[ \begin{array}{cc} \mathbf{S}\_{\mathrm{F}}^{\mathsf{T}} & \mathbf{S}\_{\mathrm{S}}^{\mathsf{T}} & \mathbf{S}\_{\mathrm{R}}^{\mathsf{T}} \end{array} \right]^{\mathsf{T}} \in \Omega \subset \mathbb{R}^{N},\tag{5.26}$$

where the associated phase space (or state space) *Ω* is an *N*-dimensional manifold4 . Together with some evolution rule that characterizes the phase velocity according to

$$\frac{d\mathbf{S}}{d\tau} = \mathbf{R}\{\mathbf{S}, \Lambda\_0\}, \quad \mathbf{R} \colon \Omega \times \mathbb{R} \to \mathbb{R}^N,\tag{5.27}$$

this creates an autonomous dynamical system with the right-hand side vector field **R** depending solely on the current state and parameterized by the bearing number *Λ*0.

The fluid description5 goes back to Section 3.1.1, where density and temperature are selected as independent primary variables on which all other properties depend. After nondimensionalization and discretization, this results in the fluid state vector

$$\mathbf{S}\_{\mathbf{F}} = \begin{bmatrix} & \dots & D\_{i,j} & \dots & T\_{i,j} & \dots \end{bmatrix}^{\mathsf{T}} \in \mathbb{R}^{2\mathrm{N\_{\sf FDM}}}\_{\geq 0} = \mathbb{R}^{\begin{pmatrix} N\_{\sf F} - 1 - \mathrm{Np} \\ \geq 0 \end{pmatrix} \binom{N\_{\sf F} - 1}{}}\tag{5.28}$$

for which the evolution through time is known from Eq. (5.25) by first-order ODEs. In the case of an isothermal calculation that disregards the temperature equation, the discretized temperature is omitted and the fluid state vector is only half as large.

With regard to the equations of motion that describe the foil structure in Eq. (5.12), each of the second-order ODEs generates one position and one velocity state variable, to which the respective bristle variable is added to give the foil structure state vector

$$\mathbf{S}\_{\mathbf{S}} = \left[ \begin{array}{cccc} \dots & \mathbf{X}\_{\mathrm{n,m}} & \mathbf{X}\_{\mathrm{n,m}}' & \mathbf{Z}\_{\mathrm{n,m}} & \dots & \end{array} \right]^{\mathsf{T}} \in \mathbb{R}^{\mathsf{3NpN}\_{\mathbf{B}}}.\tag{5.29}$$

Here, the non-trivial evolutions through time of the second and third state variable for each bump can be obtained immediately from Eqs. (5.12) and (5.14), respectively. When the regularized Coulomb friction model is used instead of a bristle model, there is of course no need to include bristle variables in the foil structure state vector.

4 The restriction of **S** to a subset of **R***N* relates to the requirement of non-negative *Di*,*j* and *Ti*,*j* in Eq. (5.28).

5 Due to the symmetry of the problem, it is sufficient to consider fluid and structure of only one bearing.

Coming to the equations of motion of the rotor in Eq. (5.17), it is once more indicated to introduce position and velocity state variables for each of the second-order ODEs. Furthermore, it is usually purposeful and therefore required by Eq. (5.27) that all the ODEs are autonomized, so the resulting dynamical system is no longer explicitly time-dependent via the unbalance. Following a suggestion by Doedel et al. [Doe+02], this can easily be achieved by the introduction of two artificial nonlinear oscillators

$$O\_{\sin}^{\prime} = +\Lambda\_0 \mathcal{O}\_{\cos} - \mathcal{O}\_{\sin} \left( \mathcal{O}\_{\sin}^2 + \mathcal{O}\_{\cos}^2 - 1 \right) \tag{5.30a}$$

$$O\_{\rm cos}^{\prime} = -\Lambda\_0 O\_{\rm sin} - O\_{\rm cos} \left( O\_{\rm sin}^2 + O\_{\rm cos}^2 - 1 \right) \tag{5.30b}$$

and by replacing harmonic forcing in Eq. (5.17) by the asymptotically stable solutions

$$O\_{\sin} = \sin(\Lambda\_0 \tau), \qquad \qquad O\_{\cos} = \cos(\Lambda\_0 \tau). \tag{5.31}$$

Altogether, this yields the rotor state vector

$$\mathbf{S}\_{\mathsf{R}} = \left[ \begin{array}{ccccc} \ X & \ X' & Y & Y' & X\_{\mathsf{R}} & X'\_{\mathsf{R}} & Y\_{\mathsf{R}} & Y'\_{\mathsf{R}} & O\_{\sin} & O\_{\cos} \end{array} \right]^{\mathsf{T}} \in \mathbb{R}^{10},\tag{5.32}$$

of which the evolution through time is governed by Eqs. (5.17) and (5.30), respectively. If no unbalance is considered, obviously no autonomizing oscillators are required, just as there is no need for disk positions or velocities when dealing with rigid rotors.

Even though the most general form of Eq. (5.27) suggests that each component of **R** may depend explicitly on all of the *N* state variables, this is obviously not the case regarding the underlying equations. On the contrary, relevant coupling mechanisms can be broken down to the relatively simple interaction chart depicted in Fig. 5.2, which clarifies that there is no direct relation between the rotor and the foil structure. Most remarkably, the identification of physical quantities to be exchanged across domain interfaces allows for a modular design of the implemented research code, so that three submodels can be merged arbitrarily into one coupled overall model. Besides the most sophisticated modeling approaches on which this work is focused, the established collection of fluid –structure – rotor submodels includes, for instance, (visco-)elastic foundation models [LBS17b] and clamped rotor shaft models [LBS18]. Moreover, it is strongly emphasized that the described modular concept refers only to the assembly of governing equations, which are then solved simultaneously using a rigorous monolithic approach not to be confused with partitioned co-simulation. In this context, a recent contribution by Bou-Saïd et al. [Bou+20] gives a comparison of the two basic methods, which outlines the advantages of monolithic approaches in terms of precision and numerical stability, but also mentions possible limitations. Especially for less abstracted GFB geometries or even more complex fluid behavior, there would practically be no alternative but to choose a partitioned solution strategy.

**Figure 5.2:** Interaction chart illustrating the modular concept of the multiphysics model and underlying coupling mechanisms between the fluid submodel, the foil structure submodel, and the rotor submodel.

The present thesis is not intended to provide an exhaustive documentation of the implemented research code GFB++, but here are some insights into fundamental technical aspects and underlying ideas. Basically, there are two project-specific main requirements that affect the choice of programming languages and environments, which are a sufficient computational speed to deal with the relatively complex model and an object-oriented programming paradigm to facilitate a modular conception. Therefore, the core code that evaluates the model equations is programmed in C/C++ and completed by Unix shell scripts for sequential program runs and post-processing. Essentially, when inheriting fluid –structure– rotor submodels from abstract classes that principally ensure the accessibility of all coupling quantities according to Fig. 5.2, solution algorithms can treat them as black boxes without knowing internal details. Finally, some relevant program parts are actually based on AUTO 2000 [Doe+02], which is a widely used package for solution continuation and bifurcation analysis, as well as on the C++ libraries Eigen [GJ+10], odeint [AM11], and CoolProp [Bel+14]. On a standard desktop computer, the expected computational times typically range from several minutes up to only a few hours for comprehensive parameter studies.

# **6 Nonlinear Effects in GFB –Rotor Systems**

With the mathematical formalization of all model equations as one dynamical system, the associated theory provides a framework of analytical and numerical methods that are very extensively covered in the literature (see, e.g., [GH83, Kuz04, MV11]). As a starting point for analyzing GFB– rotor systems, a nominal set of parameters is established before investigating the general stability and bifurcation behavior. Hereafter, particular features of refrigerant lubrication and also some possibilities of vibration reduction are thoroughly examined in view of the overall performance. In doing so, this closing chapter not only demonstrates the merits of the adopted modeling approach, solution strategy, and implementation, but also guides through a representative selection of the most important nonlinear effects and phenomena.

## **6.1 Specification of Nominal Parameter Values**

The fluid–structure– rotor model and therefore the dynamical system in Eq. (5.27) depend on a large number of nondimensional parameters, for which physically meaningful values can be deduced from the corresponding dimensional parameters. As already mentioned at the end of Section 4.1 with regard to the properties of the foil structure, the present thesis is loosely based on GFB specifications published by Ruscitto et al. [RMG78], which are transferred here unchanged to three-pad bearings. For parameters with greater uncertainties, such as the heat transfer coefficient or the friction coefficient, moderate values that can reasonably be varied in both directions are chosen with some intuition. As further starting points for parameter variations, the ambient state corresponds initially to normal conditions and the standard rotor is as generic as those in preceding rotordynamic studies of Baum [Bau18] and others. A complete overview of nominal values for all dimensional parameters can be found in Table 6.2a, with exception of the rotational speed that is later varied systematically. It should of course be noted that, for instance, the nominal mass eccentricity is only relevant if an unbalance is considered at all and it can simply be ignored otherwise.

**Table 6.1:** Nominal parameter values used as reference for numerical analyses.

**(a)** Dimensional formulation (available precision varies between parameters).

**(b)** Nondimensional formulation (rounded to three significant digits at most).

Before the dimensional parameter values from Table 6.2a can actually be inserted into Eqs. (5.4), (5.9), (5.10), (5.13), (5.15), and (5.18) to perform the nondimensionalization, an evaluation of Eq. (3.7) must yield the dependent parameter value *µ* ◦ ≈ 11.6 µPa · s. Together with all the naturally dimensionless parameters, this finally results in the nominal set of nondimensional parameter values that Table 6.2b lists in its entirety. Reversely, the rotational speed being only expressed in the nondimensional form *Λ*0, it reveals purposeful to explicitly calculate1 the characteristic time scale *t*0 ≈ 0.1 ms from Eq. (5.8), thus allowing to perform the back-transformation *n*0 = *Λ*0/(2*πt*0). As a rule of thumb for the interpretation of all following nondimensional results, this implies that a bearing number of *Λ*0 = 0.3 roughly corresponds to 30 000 rpm.

## **6.2 Application of Dynamical Systems Theory**

In the following introduction to mathematical fundamentals of dynamical systems, it is intended to always keep a close reference to the actual GFB – rotor model and to immediately apply the theory to problems with mostly nominal parameter values. Even though many effects could also be shown with simpler modeling approaches, all of the results presented hereafter are obtained assuming non-ideal gas behavior under non-isothermal conditions and taking into account foil compliance and friction.

#### **6.2.1 Classification of Solution Types**

For any possible starting point **S**0 ∈ *Ω* in the phase space of the dynamical system, which is specified as initial condition **S**(*τ* = *τ*0) = **S**0 at an arbitrary initial time *τ*0, the solution *τ* 7→ **S**(*τ*) of Eq. (5.27) directly describes the arising transient process. Nevertheless, the primary interest in theoretical and practical rotor dynamics often consists only in the long-term system behavior that can be observed for *τ* → ±∞. In this context, the relevant analyses are usually concerned with limit sets toward which dynamical systems may evolve as time goes to positive or negative infinity. Topologically speaking, such limit sets *Ω*±∞ ⊆ *Ω* are subsets of the phase space and can be identified as points, curves, manifolds, or non-trivial fractal structures. More specifically, a limit set is called attractor *Ω*+∞ if it is reached asymptotically for *τ* → +∞ after starting from a surrounding region known as basin of attraction, while it is called repeller *Ω*−∞ if it has that property for *τ* → −∞ (see, e.g., [MV11]).

1 The found numerical value for the characteristic time scale remains valid throughout this entire chapter since *p* ◦ and *µ* ◦ are basically constants for a given fluid and neither *R* nor *C* are to be varied anywhere.

**Figure 6.1:** Plots of transient rotor disk trajectories, which are obtained by numerical time integration and give a first impression of solution types toward which the dynamical system may evolve in the long term.

In a corresponding numerical experiment, the first 200 ms of transient rotor drops are simulated by numerically integrating Eq. (5.27) over the time interval *τ* ∈ [0, 2000] using a fourth-order Rosenbrock method for systems of stiff ODEs (see, e.g., [HW96]). To create a simple starting point **S**0, the fluid is initialized with ambient properties, while the foil structure and the centrally positioned rotor are undeformed and at rest. For the results, the dimension of the phase space being far too large to be visualized, the plots in Fig. 6.1 show projections2 of the disk displacements *X*R(*τ*) and *Y*R(*τ*), which may evolve in completely different ways under slightly varied circumstances.

Beginning with a rotational speed of *Λ*0 = 0.2 and setting the unbalance to *e*R = 0 as shown in Fig. 6.1a, the dynamical system and therefore the plotted displacements exhibit a transient evolution (dashed in green) that tends asymptotically toward a constant final state (green dot) and thereby reaches an attractive fixed-point solution. In Fig. 6.1b, when increasing the speed to *Λ*0 = 0.3 while still ignoring the unbalance, the transient motion (dashed in orange) now guides the dynamical system toward an attractive periodic solution or limit cycle (in orange), ending in self-excited vibrations. Finally, the addition of a mass eccentricity *e*R = 0.15 with unaltered speed in Fig. 6.1c results in a quasi-periodic motion pattern (in blue), which is due to the occurrence of two incommensurable fundamental frequencies from self-excitation and unbalance. In fact, the presence of an unbalance does not necessarily induce quasi-periodicity, but can also lead to simply-periodic, multiply-periodic, or chaotic behavior [Boy11]. Since the present thesis is mainly concerned with fixed-point and periodic solutions, the reader is referred to the thesis of Becker [Bec20] on the topic of quasi-periodicity.

2 Alternative projections of other phase variables including densities, temperatures, displacements, and velocities could of course be plotted in direct analogy and would obviously lead to the same conclusions.

#### **6.2.2 Stability and Bifurcation Analysis**

By pursuing the approach of brute-force time integration as in the previous section, it would be possible to get a more or less random overview of the solution behavior by just varying the rotational speed in small steps from one simulation to another. Regrettably, such a direct strategy turns out to be not only computationally inefficient, but principally unsuitable for several reasons when dealing with dynamical systems. First of all, time simulations can only give evidence for attracting limit sets *Ω*+∞, even though this shortcoming may be cured by additionally inverting the time axis *τ*. Most importantly, however, it is commonly observed in nonlinear dynamical systems that there are two or more coexisting limit sets, each of which can then only be found by time integration when starting by chance within the associated basin of attraction. In this perspective, more systematic investigations may be conducted by making use of numerical continuation methods, which are essentially based on the observation that the continuous variation of a parameter constructs a family of related solutions. The following gives a brief explanation of the underlying theory and methodology based on some preliminary numerical results and without claiming to be exhaustive.

For the initialization of a specific solution branch, a small initial bearing number *Λ*F1 0 is chosen such that there is a fixed-point solution **S** F1 of Eq. (5.27) to fulfill (with *i* = 1)

$$\mathbf{R}\left\{\mathbf{S}^{\mathrm{Fi}},\Lambda\_0^{\mathrm{Fi}}\right\}=\mathbf{0}.\tag{6.1}$$

This is achieved either by trying if time integration runs into a non-evolving state or, more targeted, using an under-relaxed3 Newton –Raphson method (see, e.g., [FPS20]). From a practical point of view, finding a very first system state that solves Eq. (6.1) constitutes the most tricky part of each analysis since the Newton –Raphson method requires a good initial guess, which is more difficult the smaller the bearing number. Once a solution is known together with the corresponding bearing number, however, a predictor–corrector method can deliver a sequence of further fixed-point solutions

$$\left( \left( \mathbf{S}^{\mathrm{Fi}} , \Lambda\_0^{\mathrm{Fi}} \right) \right)\_{i \in \mathbb{N}} = \left( \left( \mathbf{S}^{\mathrm{F1}} , \Lambda\_0^{\mathrm{F1}} \right) , \left( \mathbf{S}^{\mathrm{F2}} , \Lambda\_0^{\mathrm{F2}} \right) , \dots \right) \tag{6.2}$$

the existence of which is assured by the implicit function theorem (see, e.g., [Kuz04]). For additional hints on the underlying procedure of pseudo-arc-length continuation, the interested reader is also referred to the documentation of AUTO 2000 [Doe+02].

Especially for practical applications, not only the existence and identification but also the stability of solutions is an important concern, which decides on whether some

3 In order to improve the chance of convergence while not unnecessarily slowing down the calculation, the initial relaxation factor of 0.1 is increased successively during the final iterations until 1.0 is reached.

predicted behavior can actually be observed and is not of purely academic interest. In all of the following discussions, stability is understood in the sense of Lyapunov and therefore refers to perturbations in the initial conditions of a known solution. Denoting a perturbed solution by **S**∗(*τ*), this leads to the mathematical definition4 that the solution of interest **S**(*τ*) is called Lyapunov-stable if (see, e.g., [Hag78])

$$\begin{aligned} \forall \varepsilon > 0, \quad \exists \delta(\varepsilon) > 0 \colon \quad \left\| \mathbf{S}\_{\*} (\boldsymbol{\tau} = \boldsymbol{\tau}\_{0}) - \mathbf{S} (\boldsymbol{\tau} = \boldsymbol{\tau}\_{0}) \right\| < \delta(\varepsilon) \\ \implies \left\| \mathbf{S}\_{\*} (\boldsymbol{\tau}) - \mathbf{S} (\boldsymbol{\tau}) \right\| < \varepsilon, \quad \forall \boldsymbol{\tau} \ge \boldsymbol{\tau}\_{0}. \end{aligned} \tag{6.3}$$

For any arbitrarily small neighborhood5 *ε* > 0 around the solution **S**(*τ*), this requires that all trajectories of perturbed solutions **S**∗(*τ*) that start at *τ* = *τ*0 from a sufficiently small neighborhood *δ*(*ε*) > 0 remain within the *ε*-neighborhood for all times *τ* ≥ *τ*0. For practical purposes, this definition of stability needs to be accompanied by some concept of attraction, which is introduced mathematically by (see, e.g., [Hag78])

$$\exists \delta > 0 \colon \quad \left\| \mathbf{S}\_{\*} (\boldsymbol{\tau} = \boldsymbol{\tau}\_{0}) - \mathbf{S} (\boldsymbol{\tau} = \boldsymbol{\tau}\_{0}) \right\| < \delta \implies \lim\_{\boldsymbol{\tau} \to +\infty} \left\| \mathbf{S}\_{\*} (\boldsymbol{\tau}) - \mathbf{S} (\boldsymbol{\tau}) \right\| = 0. \tag{6.4}$$

Any solution that satisfies both Eq. (6.3) and Eq. (6.4) is called asymptotically stable and thereby represents an attractor *Ω*+∞ of the dynamical system (see, e.g., [NB95]). Completing these definitions, a solution is called unstable if it is not Lyapunov-stable.

In order to evaluate Eqs. (6.3) and (6.4) for the fixed-point solutions **S** F*i* from Eq. (6.2), it is permissible by virtue of the Hartman –Grobman theorem to linearize Eq. (5.27) by simply calculating the Jacobians **J** F*i* **R** for all hyperbolic equilibria (see, e.g., [MV11]). As an efficient criterion, denoting by *λ* F*i n* the *N* complex eigenvalues of **J** F*i* **R** such that

$$\det\left(\mathbf{J}\_{\mathbf{R}}^{\mathrm{Fi}} - \lambda\_n^{\mathrm{Fi}}\mathbf{I}\right) = 0, \quad n = 1, \ldots, N,\tag{6.5}$$

it comes that a fixed-point solution **S** F*i* is asymptotically stable if (see, e.g., [MV11])

$$\max\_{n=1,\ldots,N} \Re \left(\lambda\_n^{\text{Fi}}\right) < 0. \tag{6.6}$$

In the contrary case, the solution of interest is principally assumed to be unstable even though there are undecidable cases in which the largest real part is just zero.

A convenient way to make the behavior of a dynamical system visually accessible consists in the establishment of bifurcation diagrams, in which the progressions of selected solution-dependent quantities are plotted against a bifurcation parameter. With regard to rotating machinery, the displacements of the rotor are certainly among

4 In definitions, some quantities such as *ε* and *δ* are introduced independently of the actual nomenclature.

5 The notion of neighborhood is defined in terms of a suitable norm, which may be the Euclidean norm.

the most natural variables to be observed when assessing the overall performance. To this effect, the bifurcation diagram in the upper part of Fig. 6.2a traces the extrema min*τ Y*R(*τ*) and max*τ Y*R(*τ*), which gives the envelope of vertical rotor disk motions, as curves that are locally parameterized by the nondimensional rotational speed *Λ*0. For any of the fixed-point solutions in Eq. (6.2), the disk remains at a specific place for all times and the represented envelope therefore degenerates into a single curve. In this particular plot, the nominal configuration according to Table 6.1 (in orange) is opposed to cases without (in green) and with increased rotor damping (in blue), which leaves the fixed-point solution branch in the middle completely unaffected6 . As the bearing number grows, the rotor approaches asymptotically some position near the bearing center and thereby defines characteristic locus curves in Fig. 6.2b that are highlighted (in yellow) for the disk (top row) and the journals (bottom row). Moreover, all the plots in Fig. 6.2 contain information about the stability of solutions, which means that asymptotically stable fixed-point solutions according to Eq. (6.6) correspond to solid lines, whereas unstable solutions are shown as dashed lines. Obviously, there seems to be a critical bearing number *Λ*HB 0 for fixed-point solutions to become unstable, which is shifted to higher speeds the more the rotor is damped.

When following the complex eigenvalues of the Jacobian in a closer vicinity of *Λ*HB 0 , it is found that a pair of complex conjugate eigenvalues crosses the imaginary axis, which is a characteristic finding during the occurrence of a Hopf bifurcation (HB). Computationally, such HB points can efficiently be detected by introducing specific test functions to be evaluated as part of numerical continuation (see, e.g., [Doe+02]). In addition to the loss of stability for fixed-point solutions, a HB also marks the birth of a branch of periodic solutions **S** P*i* (*τ*), all of which satisfy Eq. (5.27) according to

$$\frac{\mathbf{d}}{\mathbf{d}\tau}\mathbf{S}^{\mathrm{Pi}}(\tau) = \mathbf{R}\left\{\mathbf{S}^{\mathrm{Pi}}(\tau), \Lambda\_0^{\mathrm{Pi}}\right\}, \quad \mathbf{S}^{\mathrm{Pi}}(\tau) = \mathbf{S}^{\mathrm{Pi}}(\tau + \tau^{\mathrm{Pi}}).\tag{6.7}$$

For periodic solutions of autonomous dynamical systems, the period *τ* P*i* in Eq. (6.7) is a priori unknown and must thus be found by an appropriate solution strategy. Detailed descriptions of the orthogonal collocation method used here can be found in the book by Kuznetsov [Kuz04], which also presents alternative shooting methods. Having calculated for *Λ*P1 0 ≈ *Λ*HB 0 the first periodic solution **S** P1(*τ*) with period *τ* P1 , a predictor–corrector method similar to that from Eq. (6.2) can provide a sequence

$$\left( \left( \mathbf{S}^{\mathrm{Fi}}(\tau), \Lambda\_0^{\mathrm{Pi}}, \tau^{\mathrm{Pi}} \right) \right)\_{i \in \mathbb{N}} = \left( \left( \mathbf{S}^{\mathrm{Pl}}(\tau), \Lambda\_0^{\mathrm{Pl}}, \tau^{\mathrm{Pl}} \right)\_{\prime} \left( \mathbf{S}^{\mathrm{Pl}2}(\tau), \Lambda\_0^{\mathrm{Pl}2}, \tau^{\mathrm{Pl}} \right), \dots \right), \tag{6.8}$$

which is constructed by numerical continuation until a defined end point is reached or until convergence of the orthogonal collocation method fails (see, e.g., [Doe+02]).

6 For better distinction, the three fixed-point solution branches are slightly shifted relative to each other.

**(a)** Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom).

**(b)** Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row).

**Figure 6.2:** Results of numerical continuation and bifurcation analysis for different damping coefficients (solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations/bifurcating solutions). The assessment of stability according to Eqs. (6.3) and (6.4) is slightly more difficult for the periodic solutions **S** P*i* (*τ*) from Eq. (6.8) than it is for fixed-point solutions. Using again the idea of linearization to judge the evolution of an initial perturbation, the Jacobian **J** P*i* **R** (*τ*) comes out as a time-dependent matrix with periodic coefficients. For dealing with the systems of linear ODEs with periodic coefficients that arise, the mathematical foundation is given by classical Floquet theory (see, e.g., [Hag78]). When calculating the fundamental matrix solution **M**P*i* **R** (*τ*) that is related to **J** P*i* **R** (*τ*) by

$$\frac{\mathbf{d}}{\mathbf{d}\tau} \mathbf{M}\_{\mathbf{R}}^{\mathrm{Pi}}(\tau) = \mathbf{J}\_{\mathbf{R}}^{\mathrm{Pi}}(\tau) \mathbf{M}\_{\mathbf{R}}^{\mathrm{Pi}}(\tau), \quad \mathbf{M}\_{\mathbf{R}}^{\mathrm{Pi}}(\tau=0) = \mathbf{I},\tag{6.9}$$

an evaluation after just one period yields the monodromy matrix **M**P*i* **R** (*τ* = *τ* P*i* ). The *N* complex eigenvalues *λ* P*i n* of the monodromy matrix, which are also referred to as the Floquet multipliers of the periodic solution, can ultimately be determined by

$$\det\left[\mathbf{M}\_{\mathbf{R}}^{\mathrm{Pi}}(\boldsymbol{\tau}=\boldsymbol{\tau}^{\mathrm{Pi}})-\lambda\_{n}^{\mathrm{Pi}}\mathbf{I}\right]=\mathbf{0},\quad n=1,\ldots,\mathrm{N}.\tag{6.10}$$

As a characteristic of autonomous dynamical systems, there is always one special Floquet multiplier *λ* P*i* 1 = 1 that represents tangential perturbations along an orbit, which can be visualized graphically by means of a corresponding Poincaré map. Based on the remaining (*N* − 1) Floquet multipliers, the Andronov–Witt theorem states that a periodic solution **S** P*i* (*τ*) is asymptotically stable if7 (see, e.g., [MV11])

$$\max\_{n=2,\ldots,N} \left| \lambda\_n^{\text{Pi}} \right| < 1. \tag{6.11}$$

Similar to Eq. (6.6), the periodic solution is considered as being unstable otherwise, knowing that there may be again cases that cannot be decided by a linearized theory.

When finally coming back to the bifurcation diagram in the upper part of Fig. 6.2a, upper and lower envelopes of vertical disk motions are obviously no longer identical for the periodic solutions that branch from the respective HB points (circle symbols). Also, the exemplary trajectories of the disk (top row) and the journals (bottom row) plotted in Fig. 6.2b now correspond to closed curves around each of the HB points, which illustrate the orbit8 shapes in relation to the nominal clearance (dotted line). Due to the compliance of the foil structure, it is clear that the journals can travel to eccentricities *e*(*τ*) ≥ 1, which are even amplified for disk motions by the elastic shaft. By convention, asymptotic stability according to Eq. (6.11) is indicated by different line types for stable (solid lines) and for unstable (dashed lines) periodic solutions.

7 Here, it is pointed out that | . . . | denotes the absolute value √ (<2 (. . .) + =2 (. . .)) of a complex number.

8 Strictly speaking, the two-dimensional closed curves are not the actual orbits but only their projections.

As an additional information, the bifurcation diagram in the lower part of Fig. 6.2a visualizes all corresponding periods as nondimensional subsynchronous frequencies

$$A\_{\rm sub}^{\rm Pi} = \frac{2\pi}{\pi^{\rm Pi}}\tag{6.12}$$

in relation to *Λ*0. This implies *Λ*P*i* sub/*Λ*0 ≈ 0.5 for typical half-frequency whirling, whereas the frequency ratio would become *Λ*P*i* sub/*Λ*0 ≈ 1 for synchronous vibrations.