faculty of mathematics and natural sciences

## A computational fluid

## dynamics analysis on air bearing designs

### Bachelor Project Mathematics

June 2015 Student: I.W. Tan

First supervisor: Prof.dr.ir. R.W.C.P. Verstappen

**Abstract**

This report treats the results of a research project on simulating the air flow in aero-
static thrust air bearings with the Open Source Computational Fluid Dynamics (CFD)
software package OpenFOAM. The type of air bearing considered in this research
is a simple type circular bearing with eight orifices that are equally distributed over
a concentric circle. In this report we first present some background information on
air bearings, treat the fluid dynamics that is involved and consider the previously
conducted research on computational models of air bearings. Next, we will discuss
the simulations that we performed where we varied the choice of solvers, boundary
conditions, turbulence models and the viscosity parameter. Lastly, we will analyze the
results and give suggestions for future work on this subject. We find that manually
*altering the viscosity via the Sutherland Coefficient A** _{s}* provides a stable initial flow,
from which we can incrementally lower the viscosity back to values closer to the realis-
tic value. Furthermore we find that the solvers sonicFoam and rhoCentralFoam show
similar behaviour on our meshes and that sonicFoam and rhoSimpleFoam in Open-
FOAM 1.7 perform differently from sonicFoam and rhoSimpleFoam in OpenFOAM
2.3.

*Keywords: aerostatic thrust air bearing, OpenFOAM, CFD*

**Contents**

**1 Introduction** **1**

**2 Background** **3**

2.1 Flat aerostatic thrust air bearings . . . 3

2.2 Fluid dynamics . . . 4

2.2.1 General conservation principles . . . 4

2.2.2 Supersonic flow . . . 5

**3 Model** **9**
3.1 OpenFoam . . . 9

3.2 Mesh . . . 10

3.3 Solvers . . . 10

3.3.1 rhoSimpleFoam . . . 10

3.3.2 sonicFoam . . . 14

3.3.3 rhoCentralFoam . . . 16

3.4 Turbulence models . . . 18

3.4.1 *The Realizable k-ε model . . . .* 19

3.4.2 *k-ω SST . . . 21*

3.5 Geometry and boundary conditions . . . 23

3.6 Schemes and Transport model . . . 24

3.7 Postprocessing . . . 26

**4 Results** **27**
4.1 Adapting the mesh . . . 27

4.2 Changing several settings . . . 34

4.2.1 Initial pressure . . . 34

4.2.2 Adjusting the time-step size . . . 35

4.3 Adjusting the viscosity . . . 36

4.4 Comparison to OpenFOAM 1.7 and a different geometry . . . 42

**5 Conclusion** **47**

**6 Future Work** **51**

**A Additional results** **I**

**Chapter 1** **Introduction**

This report treats the results of a research project on numerical simulations of aero- static thrust air bearings as commissioned by the company Schut Geometrical Metrol- ogy. Air bearings are bearings that use a thin film of air as a lubricant to reduce the friction between two moving surfaces. The company wishes to use the results from these simulations to improve the bearing design. Ultimately, the bearing has to satisfy certain stability properties required for the implementation in high precision measurement devices.

Figure 1.1: Three types of aerostatic circular thrust air bearings; 1.

multiple orifices 2. single orifice 3. porous medium

In Figure 1.1, three types of flat aerostatic air bearings are illustrated. In the past, the company has carried out several successful simulations on air bearings of type 2, which contains a single inlet orifice in the center[6]. However, applying the same modelling procedure to the geometry of the bearing of type 1 did not yield any realistic results. Therefore, the goal of this project was to find out how a bearing of type 1 can be simulated correctly. The simulations were carried out in the Open

Source Computational Fluid Dynamics (CFD) software package OpenFOAM 2.3.0.

This software provides features for applications ranging from complex turbulent flow involving chemical reactions to solid dynamics and electromagnetics[2].

When simulating fluid flow, it is important to tailor the numerical model to the flow dynamics involved. By making the right assumptions and simplifications, a lot of computation time can be saved without losing too much accuracy. After studying some of the available literature, we found that as the air flow exits the nozzle at the bottom of the bearing, it is subject to an extreme decrease in flow area. This can result in locally supersonic flow, which requires the use of algorithms for compressible fluids. This is presented in Chapter 2, along with a description of the bearing and the governing equations.

In reality, inevitable disturbances, such as disturbances in the load that is exerted on the bearing, will affect the height of the air gap. However, in this project we investigated a simplified situation in which the gap height is fixed. We have tried multiple strategies to obtain a stable and realistic simulation, such as changing the boundary conditions, trying different solvers and manually altering the viscosity. The specific settings and procedures of the simulations are explained in Chapter 3.

The results from our simulations are presented in Chapter 4 together with an analysis of the results in Chapter 5. We find that running a simulation with an altered viscosity is an effective strategy to obtain a stable initial flow for lower and more realistic viscosities. Furthermore, we find that the solvers sonicFoam and rhoSimpleFoam in OpenFOAM 1.7 perform differently from sonicFoam and rhoSimpleFoam in OpenFOAM 2.3. Lastly, we find that the solvers rhoCentralFoam and sonicFoam show similar results in nearly all configurations, despite their different nature.

Finally, we present some suggestions for future work in Chapter 6.

**Chapter 2** **Background**

**2.1** **Flat aerostatic thrust air bearings**

Air bearings exist in different types, such as aerodynamical air bearings, which use their movement to create a thin high pressure air film between the surfaces, or aerostatic air bearings, in which the thin air film is externally pressurized. Moreover, there exist many other variable properties such as shape and material.

Figure 2.1: A schematic representation of the cross-section of an aero- static air bearing with a single orifice.

Source: https://www.mech.kuleuven.be/apt/intro.en.html

In Figure 2.1, the cross-section of an aerostatic air bearing with a single orifice
is shown. The shape and size of the air bearings considered in this research are
comparable to a hockey puck. An external compressor connected to the chamber
*inside the bearing, applies a supply pressure p*_{s}, which is higher than the atmospheric
*back pressure p*_{a}. The pressure difference induces an air flow which exits the chamber
via the nozzle and flows into the air gap via the orifice. Finally, as it passes the edge
of the bearing, it can escape into the open air. The air bearing considered in this

research has eight orifices distributed evenly over a concentric circle, comparable to type 1 in Figure 1.1.

**2.2** **Fluid dynamics**

**2.2.1** **General conservation principles**

Because the smallest length scale in this set up, the gap height under the bearing, is 4µm long, we can regard the gas as a continuum rather than as a collection of discrete particles [15]. This means that we can use arbitrary control volumes to derive conservation equations for mass, momentum and scalar properties such as temperature and energy. By specifying the inflow, outflow, and sources for an arbitrary control volume we obtain equations 2.1, 2.2 and 2.6 for mass, momentum and energy respectively[8].

*∂*

*∂ t*
Z

*Ω*

*ρ dΩ +*
Z

*S*

* ρv · n dS = 0,* (2.1)

Equation 2.1 represents the conservation of mass, in which *ρ represents the*
* density, v the velocity vector and n the normal vector on the surface S of the control*
volume

*Ω. It states that the amount of mass in a control volume can only change by*a mass flow through the boundary of that control volume.

*∂*

*∂ t*
Z

*Ω*

* ρv dΩ +*
Z

*S*

* ρv v · n dS =*
Z

*S*

T*_{· n dS +}*
Z

*Ω*

* ρb dΩ,* (2.2)

Equation 2.2 represents the conservation of momentum, which is derived from Newton’s second law of motion,

d(mv)

*dt* =X

* f*. (2.3)

The tensorT* in equation 2.2 represents the stress tensor and b represents the body*
forces. The equation states that the amount of momentum in a control volume can
change by advection, surface forces (for example pressure, shear stress or surface
tension) or body forces (for example gravity or the Coriolis force).

In this research we will approximate air as a Newtonian fluid. Therefore, the stress tensor can be written as,

T= −
*p*+2

3* µ div v*

I*+ 2µ*D, (2.4)

*where p is the static pressure,µ is the dynamic viscosity,*Iis the unit tensor andDis
the deformation tensor[8],

D=1

2**grad v + (grad v)*** ^{T}* . (2.5)

Equation 2.1 and 2.2 together are known as the Navier-Stokes equations.

*∂*

*∂ t*
Z

*Ω*

*ρh dΩ +*
Z

*S*

* ρhv · n dS =*
Z

*S*

*kgrad T · n dS +*
Z

*Ω*

* (v · grad p +* S

**: grad v**) dΩ+ *∂*

*∂ t*
Z

*Ω*

*p* dΩ. (2.6)

*Lastly, equation 2.6 represents the conservation of energy, in which h is the enthalpy*
*per unit mass, k is the thermal conductivity defined as,*

*k= µc*_{p}

Pr, (2.7)

*with c*_{p} *the specific heat at constant pressure and Pr the Prandtl number, T is the*
absolute temperature (K) andSis the viscous part of the stress tensorS_{=}T* _{+ p}*I

_{[8].}

**2.2.2** **Supersonic flow**

In a lot of the previously conducted research on air bearings the presence of supersonic
flow (Ma*> 1) and an undesirable pressure depression attributed to a shock wave is*
mentioned[5, 6, 14]. Since our configuration is comparable to the configurations in
those articles, it is plausible that we will encounter supersonic flow too.

**Shock waves**

An important component of supersonic flow is the possible occurrence of shock waves.

A shock wave is a thin region in the fluid across which the flow properties change drastically. Because the highest velocity at which information on a fluid’s properties travels is the speed of sound, this information cannot travel upstream in supersonic flow. For the scope of this project it is sufficient to acknowledge the existence of shock waves, in order to analyze the results of our model. More information on shock waves can be found in[4] among many other text books.

**Pressure depression**

The undesirable pressure distribution mentioned earlier, is located close to the inlet of the air bearing. This depression has a negative effect on the load capacity and has

been studied extensively both theoretically and experimentally. Some researchers [20] attribute this depression to the transition from laminar to turbulent flow, but this contradicts the majority of articles on the subject, which show that the depression exists due to shock waves in the bearing clearance[5, 7, 14]. In general, the latter is accepted.

This explanation states that the air bearing can be considered as a converging- diverging duct, since the cross-sectional flow area at the entrance of the gap is smaller than the flow area in the nozzle. From the area-velocity relation

*dA*

*A* = (Ma^{2}− 1)*dv*

*v* , (2.8)

we can see that the nature of the flow determines the effect of a change in cross-
sectional flow area. For subsonic flow (Ma*< 1), a decrease in area will accelerate the*
flow and an increase in area will decelerate the flow. For supersonic flow (Ma*> 1) a*
decrease in area will decelerate the flow and an increase in area will accelerate the
flow[4]. This explains how a small adjustment in supply pressure can induce a large
change in flow behaviour, as can be seen from figure Figure 2.2. If the difference
*between the supply pressure p*_{s}*and the atmospheric pressure p*_{a} is large enough, the
flow may become choked at the entrance of the gap and become supersonic.

Figure 2.2: A converging-diverging duct and the possible pressure distributions and velocity profiles.

Source: http://turbo.mech.iwate-u.ac.jp/Fel/turbomachines/

stanford/images/LR_5.3.gif

*In Figure 2.2 the effect of the ratio of p*_{a}*to p*_{s}on the flow in a converging diverging
duct is illustrated. The pressure distributions that we encountered most often are
*of the type d, where the pressure difference is large enough to induce supersonic*
flow and a compression shock wave exists in the diverging duct that is the bearing
clearance.

Because the transition from the inlet to the gap is at a right angle, the flow separates from the boundary at high velocities. The resulting separation bubble as illustrated in Figure 2.3 narrows the flow area even more. For a more detailed description of the fluid flow and the pressure depression often encountered in air bearings we refer to[7] and [14].

Figure 2.3: Qualitative representation of the shock/boundary layer in- teraction in the inlet region [7]

**Chapter 3** **Model**

**3.1** **OpenFoam**

OpenFOAM is a free open source software package for applications ranging from complex turbulent flow involving chemical reactions to solid dynamics and electro- magnetics[2]. In addition to the official software package that is managed by the OpenFOAM Foundation, there are numerous programmers that write extensions for the existing OpenFOAM software themselves. The most of these unofficial extensions are public too, but they have not been checked by the OpenFOAM Foundation. Often, these extensions contribute to new OpenFOAM versions after being processed by the programmers of the OpenFOAM Foundation. The software is being used for commercial as well as academic purposes and a lot of information can be found on internet fora where users discuss their problems and experiences.

To perform a simulation in OpenFOAM, the user must create a case directory with at least the following sub-directories,

1. a time directory,

2. the "constant" directory, 3. the "system" directory.

*The time directory, for example named 0 if t*_{0} = 0, contains the boundary and
initial conditions for at least all variables that will be used in the calculations. The
constant directory contains files that define the mesh as well as files in which the
thermophysical and turbulence settings are specified. The system directory contains
the files in which the discretization schemes, linear system solvers, parallelisation
settings, solver settings and control settings are specified.

**3.2** **Mesh**

The mesh that we used was provided by the commissioning company and was created
in Mathematica[19]. It exploits the symmetry of the bearing by covering only ^{1}_{8} of
the actual bearing geometry, which saves a lot of computation time. The mesh is
made of hexahedral cells exclusively because that is best suited for flow simulations
of very thin layers, such as our bearing clearance. Furthermore, by starting from a
three dimensional structured hexahedral mesh it is less difficult to adopt the explicit
mesh description format that OpenFOAM requires. The hexahedral mesh was fitted
to the geometry of the bearing by transformations and removals of redundant cells
without disturbing the structure and ordering of the mesh.

**3.3** **Solvers**

We have used three solvers from the compressible solvers that are available in Open- FOAM 2.3.0, namely

### sonicFoam

,### rhoSimpleFoam

and### rhoCentralFoam

. We will describe these three solvers below.**3.3.1** **rhoSimpleFoam**

The solver

### rhoSimpleFoam

is a steady-state Semi-Implicit Method for Pressure Linked Equations (SIMPLE) solver for laminar or turbulent Reynolds-averaged Navier- Stokes (RANS) flow of compressible fluids[2]. It uses the SIMPLE algorithm to solve the coupled equations for velocity and pressure. The SIMPLE algorithm is an iterative algorithm proven to be fairly efficient in solving steady state problems by the use of a predictor-corrector loop and under-relaxation factors[8]. This algorithm is also the basis for many other iterative Navier-Stokes algorithms such as SIMPLEC, SIMPLER, PISO and PIMPLE. In Figure 3.1 a flow chart of the algorithm is illustrated.*Start with u** ^{n−1}*,

*ρ*

^{n−1}*and p*

^{n−1}### k = 0 m = 0

Linearized momen-
tum equation *→ u*^{m∗}

Calculate mass fluxes ˙*m*^{m∗}* _{l}* at
cell faces from

*ρ*

^{m}^{−1}

*and v*

_{n}

^{m}^{∗}

Pressure correction*→ p*^{0}

Correct mass fluxes at cell faces

### k++

*Calculate u*^{m}*from p*^{m}

### m++

Calculate other variables

Proceed to next time step

### n++

Outer
ifk _{≤} loop

nNonOrthogonal
Correctors
if m *<*

nOuter Correctors

Figure 3.1: A flowchart of the SIMPLE algorithm for the *n*th time step

*To calculate the values for the nth time step, the algorithm first enters a loop*
to calculate the pressure field and the velocity field, using the variable fields from
*the previous time step as initial values. The iterations of this loop are called outer*
*iterationsto distinguish them from inner iterations which are the iterations performed*
on linear systems with fixed coefficients[8]. These outer iterations continue as long
as the pressure and velocity fields do not satisfy the specified residual tolerance and
as long as the number of performed iterations

### m

does not exceed the specified value of### nOuterCorrectors

_{[11].}

Lastly, all other variable fields, such as temperature and turbulent kinetic energy are calculated. After that, the algorithm proceeds to the next time step.

We will now follow the derivations and steps as presented in[8] to explain how the outer loop operates.

The outer loop consists of three tasks. First, in case

### m

= 0, the pressure and density*fields p*

^{n}^{−1}and

*ρ*

^{n}^{−1}from the previous time step are used to solve the segregated linearized momentum equations. Else, in case

### m

*> 0, the pressure and density*

*fields p*

*and*

^{m−1}*ρ*

*from the previous outer iteration are used instead. These*

^{m−1}*linearized momentum equations give a new estimate u*

^{m}^{∗}for the velocity field, where the asterisk is used to indicate that this estimate does not necessarily satisfy the continuity equation. The discretized segregated linearized momentum equation for

*the ith component of the velocity estimate u*

^{m}^{∗}

*at the point P that is solved in this*first step of the outer iteration can be written as,

*A*^{u}_{P}^{i}*u*^{m∗}* _{i,P}*+X

*l*

*A*^{u}_{l}^{i}*u*^{m∗}_{i,l}*= Q*^{m−1}_{u}

*i* *− ∆Ω*

*δp*^{m−1}*δx**i*

. (3.1)

*Here m is the outer iteration counter and∆Ω denotes the control volume centered*
around a cell face, which is not necessarily the same as the volume of the control vol-
*ume around point P. The coefficients A*_{P}*and A** _{l}* depend on the chosen discretization
schemes and represent the terms in the momentum equation that contain unknown

*variable values, in this case u*

^{m}

_{i}^{∗}

*. The subscript l denotes the index of the neighbouring*

*cells. The coefficient Q*

^{m−1}

_{u}*i* contains all terms of the momentum equation which do
not contain unknown variable values and thus are presumed known, except for the
pressure gradient which is written explicitly in the last term in this equation [8].

The pressure gradient is written explicitly here, because later on we want to relate a correction of the velocity and density to a correction of the pressure.

Equation 3.1 can be rewritten as,

*u*^{m∗}* _{i,P}* =

*Q*

^{m−1}

_{u}*i* −P

*l* *A*^{u}_{l}^{i}*u*^{m∗}_{i,l}*A*^{u}_{P}* ^{i}* −

*∆Ω*

*A*^{u}_{P}^{i}

*δp*^{m−1}*δx**i*

. (3.2)

Because the velocity field obtained from equation 3.2 together with the old density
values*ρ*^{n}^{−1} do not satisfy continuity, substituting the new interpolated cell face mass
fluxes ˙*m*^{∗}* _{l}* into the continuity equation,

*(ρ*^{m}^{−1}*− ρ*^{n}^{−1}*) ∆Ω*

*∆t* +X

*l*

˙

*m*^{∗}_{l}*= Q*^{∗}* _{m}*, (3.3)

*results in an imbalance Q*^{∗}* _{m}* [8]. To enforce the conservation of mass it is necessary to
impose a correction step.

This correction is the second task of the outer loop. Since we are dealing with
compressible flow, the mass flux varies not only with the velocity but with the density
as well. We must therefore calculate the density correction*ρ*^{0}and the normal velocity
*correction v*_{n}^{0} *such that the corrected mass flux for the mth outer iteration and lth*
cell face can be expressed as,

˙

*m*^{m}_{l}*= ρ*_{l}^{m}*v*_{n,l}^{m}*S*_{l}*= (ρ*^{m−1}*+ ρ*^{0})*l**(v*_{n}^{m∗}*+ v*_{n}^{0})*l**S** _{l}*, (3.4)

*where S*

_{l}*denotes the cell face area between the cell of point P and the lth neigbouring*cell. The mass flux correction in terms of density and velocity can thus be expressed as,

˙

*m*^{0}_{l}*= (ρ*^{m}^{−1}*S v*_{n}^{0})*l**+ (v*_{n}^{m}^{∗}*Sρ*^{0})*l**+ (ρ*^{0}*v*_{n}^{0}*S*)*l*. (3.5)
Though methods exist to calculate the last term in equation 3.5[8], we will neglect
it here because it is of second order in terms of corrections. For a colocated variable
arrangement, a mesh in which the velocity field and pressure field share the same
grid, the first term of equation 3.5 can be approximated by,

*(ρ*^{m}^{−1}*S v*_{n}^{0})*l**= − (ρ*^{m}^{−1}*S∆Ω)**l*

1
*A*^{u}_{P}

*l*

*δp*^{0}
*δn*

*l*

, (3.6)

*where the overbar denotes an interpolated value and n is the coordinate of the*
outward facing normal to the cell face[8]. If we assume the temperature to be fixed
for each outer iteration, we can approximate*ρ*^{0}in the second term of equation 3.5 by

*ρ*^{0}=

*∂ ρ*

*∂ p*

*T*

*p*^{0}*= C*_{ρ}*p*^{0}. (3.7)

Since we assumed that the fluid is a perfect gas, we can use the equation of state to
*determine the coefficient C** _{ρ}* by,

*C** _{ρ}*= 1

*RT*. (3.8)

The second term of equation 3.5 can be written as,
*(v*_{n}^{m}^{∗}*Sρ*^{0})*l* =* C*_{ρ}*m*˙^{∗}

*ρ*^{m−1}

*l*

*p*^{0}* _{l}*. (3.9)

If we now substitute equation 3.6 and 3.9 into equation 3.5 and neglect the last term, we get an expression for the mass flux correction in terms of known variable values and the pressure correction,

˙

*m*^{0}_{l}*= − (ρ*^{m−1}*S∆Ω)**l*

1
*A*^{u}_{P}

*l*

*δp*^{0}
*δn*

*l*

+* C*_{ρ}*m*˙^{∗}
*ρ*^{m}^{−1}

*l*

*p*^{0}* _{l}*. (3.10)

For the mass fluxes to satisfy the continuity equation, the corrected mass fluxes and the corrected densities should satisfy,

*ρ*^{0}_{P}*∆Ω*

*∆t* +X

*l*

˙

*m*^{0}_{l}*+ Q*^{∗}* _{m}*= 0. (3.11)

If we substitute equation 3.7 for*ρ*^{0}* _{P}* and equation 3.10 for ˙

*m*

^{0}

*into equation 3.11, we can write the pressure correction equation as,*

_{l}*A*_{P}*p*^{0}* _{P}*+X

*l*

*A*_{l}*p*^{0}_{l}*= −Q*^{∗}* _{m}*. (3.12)

After the resulting system of equations has been solved, only the last step of the outer loop remains. First, the corrected pressure is calculated. It has been found that adding only a portion of the pressure correction improves the convergence[8].

This is done by choosing an under-relaxation factor*α**P*∈ [0, 1] ⊂ R which acts like a
pseudo-time step, such that,

*p*^{m}*= p*^{m}^{−1}*+ α**P**p*^{0}. (3.13)

This under relaxation factor is not used in the final outer iteration because the results
of each time step must satisfy continuity. After this pressure correction, substituting
*equation 3.13 and the expression u*^{m}_{i}*= u*^{m}_{i}^{∗}*+ u*^{0}into equation 3.2 gives,

*u*^{0}* _{i,P}*= −
P

*l* *A*^{u}_{l}^{i}*u*^{0}_{i,l}

*A*^{u}_{P}* ^{i}* −

*∆Ω*

*A*

^{u}

_{P}

^{i}*δp*^{0}
*δx**i*

*P*

, (3.14)

*from which the corrected velocity u** ^{m}*is solved.

Via the file

### fvSolution

the user can control several settings for the SIMPLE algorithm. For example, the maximum number of outer iterations can be set through the variable### nOuterCorrectors

. Also, it is possible to solve the pressure correction equation and calculate the mass fluxes at the boundary cell faces multiple times by setting the variable### nNonOrthogonalCorrectors

larger than zero, to account for possible grid non-orthogonality. The under-relaxation factors and residual tolerances can also be specified there.**3.3.2** **sonicFoam**

The solver

### sonicFoam

in OpenFOAM version 2.3 uses the PIMPLE algorithm, which is a hybrid between the SIMPLE and PISO algorithms. The PISO algorithm is very sim- ilar to the SIMPLE algorithm, which we described above, except for two differences.It executes the outer loop only once for each time step and can therefore not use under-relaxation factors, but on the other hand the pressure correction and resulting velocity field are calculated more than once per time step, which we will call the PISO loop from now on. The PIMPLE algorithm combines SIMPLE and PISO such that both the outer loop and the PISO loop can be executed multiple times per time step. For steady flows this allows larger time steps, as with the SIMPLE algorithm, since we can use the under-relaxation factor from the SIMPLE algorithm which acts as a pseudo time step.

*Start with u*^{n}^{−1}, *ρ*^{n}^{−1} *and p*^{n}^{−1}

### k = 0

,### l = 0

,### m = 0

Linearized momen- tum equation*→ u*

^{m}^{∗}Calculate mass fluxes ˙

*m*

^{m}

_{l}^{∗}at cell faces from

*ρ*

^{m−1}*and v*

_{n}

^{m∗}Pressure correction*→ p*^{0}

Correct mass fluxes at cell faces

### k++

*Calculate u*^{m}*from p*^{m}

### l++

,### m++

Calculate other variables

Proceed to next time step

### n++

Outer PISO loop

loop ifk ≤

nNonOrthogonal Correctors

ifl *<*

nCorrectors
ifm *<*

nOuter Correctors

Figure 3.2: A flowchart of the PIMPLE algorithm for the *n*th time step

In Figure 3.2 a flow chart of the PIMPLE algorithm is shown. The results from the
*n−1th time step are used as initial values for the nth time step as the algorithm enters*
the outer loop. First the discretized linearized segregated momentum equations are

*solved to obtain an estimate u*^{m}^{∗}, where again the asterisk is used to indicate that this
value does not necessarily satisfy the law of conservation of mass. Together with*ρ** ^{n−1}*,

*this estimate u*

*is then used to calculate an estimate of the mass fluxes through the cell faces ˙*

^{m∗}*m*

^{m}

_{l}^{∗}. After that, the flux estimates are used to solve the imbalance of the continuity equation and subsequently to solve the pressure correction equation.

*This results in a correction on the pressure p*^{0}, a correction on the density*ρ*^{0}and a
*correction on the cell face normal velocity v*_{n}^{0}, by which the mass fluxes at the cell faces
can be corrected. Depending on the variable

### nNonOrthogonalCorrectors

, these two steps can be executed multiple times to account for non-orthogonality. After that,*the velocity cell center value estimates u*

^{m}

_{P}^{∗}are corrected and depending on the value of

### nCorrectors

the algorithm can perform another PISO iteration by returning to step 2 of the outer loop. When the PISO loop is finished, the algorithm can perform another outer iteration, depending on the value of the variable### nOuterCorrectors

. We see that the PIMPLE algorithm is very similar to the SIMPLE algorithm. For a more detailed description of the PISO loop we refer to[8, 17].**3.3.3** **rhoCentralFoam**

The solver

### rhoCentralFoam

is a density based solver for high speed compressible flow[12]. Similarly to the SIMPLE based algorithms it is also suitable for implicit time-stepping. But whereas SIMPLE based algorithms use a pressure correction to allow this, the### rhoCentralFoam

solver uses an operator-splitting approach involving a diffusion correction to allow this. This implies that in### rhoCentralFoam

for the velocity and energy first an estimate based on the inviscid fluxes alone is calculated, which is then corrected by a diffusion correction. We will now present a more detailed description based on the article written by the programmers who have implemented### rhoCentralFoam

in OpenFOAM[9].**First, the values of u and** *ρ of the previous time step are used to compute the*
density-weighted field ˆ**u*** = uρ. By using the Van Leer limiter and the value T from*
the previous time step, the cell face values ˆ

**u**

_{f}_{±},

*ρ*

*f*±

*and T*

_{f}_{±}for the convective terms are then computed by interpolation for inward and outward flux separately. Here,

*the subscripts f*

**+ and f − denote the flow in the direction +S***f*and

**−S**

*f*respectively.

**The vector S*** _{f}* is defined as the face area vector normal to the boundary face surface,
pointing out of the owner cell with a magnitude equal to the area of the boundary
face. After this, the remaining variables are calculated via the expressions as denoted
in table 3.1.

**u**_{f}_{±} = **u**ˆ_{f}_{±}
*ρ**f*±

*φ**f*± **= S***f* **· u***f*±

*p*_{f}_{±} *= ρ**f*±*RT*_{f}_{±} *c*_{f}_{±} *= ÆγRT**f*±

**T**_{exp} **= µ[(∇u)*** ^{T}* −

^{2}

_{3}tr

**(∇u)I]**

Table 3.1: The variable expressions for the cell face values by inward and outward flux [9].

Using these variables and updated*µ and k, the continuity equation is solved for ρ.*

Next, the inviscid convective momentum equation,

**∂ ˆu**

*∂ t*

*I*

* + ∇ · [uˆu] + ∇p = 0,* (3.15)

is solved for ˆ**u. This step is also known as the inviscid momentum prediction. By**
using this new* ρ and ˆu we then update u after which it is corrected by solving the*
diffusive velocity correction equation,

**∂ (ρu)**

*∂ t*

*V*

* − ∇ · (µ∇u) − ∇ · (T*exp) = 0, (3.16)

**for u, where T**

_{exp}, the part of the viscous stress tensorScontaining inter-component

**coupling, is expressed explicitly in terms of the old value for u and the Laplacian term**

* ∇ · (µ∇u), the part of the viscous stress tensor without inter-component coupling is*
treated implicitly.

Now, only the energy, temperature and pressure have to be calculated. This is done by first solving the equation,

*∂ ˆE*

*∂ t*

*I*

* + ∇ · [u(ˆE + p)] + ∇ · (T · u) = 0,* (3.17)

for ˆ*E. We can then compute the temperature T from the equation,*
*T* = 1

*c*_{v}

*E*ˆ
*ρ* −**|u|**^{2}

2

, (3.18)

with our updated* ρ, ˆE and u. After that, the temperature T is corrected by solving*
the diffusive temperature correction,

*∂ (ρc**v**T*)

*∂ t*

*V*

*− ∇ · (k∇T ) = 0,* (3.19)

*for T . Lastly the pressure p is updated by p= ρRT.*

**3.4** **Turbulence models**

To take into account that the flow is subject to turbulence, we used Reynolds-Averaged Navier-Stokes (RANS) turbulence models. For these models the Reynolds-Averaged Navier-Stokes equations are used which are slightly different from the Navier-Stokes equations we presented in chapter 2. For incompressible flow, the equations are obtained by decomposing the expressions for the velocity, pressure and viscous stress tensor into a mean component and a fluctuating component such that,

*u*_{i}*= 〈u**i**〉 + u*^{0}* _{i}*, (3.20)

*p= 〈p〉 + p*^{0}, (3.21)

S* _{i j}* = 〈S

_{i j}_{〉 +}S

^{0}

*. (3.22) Here*

_{i j}*〈u*

*i*

*〉, 〈p〉 and 〈*S

*〉 denote the time averages defined as,*

_{i j}*〈u**i*〉 = 1
*T*

Z *t**+T*
*t*

*u*_{i}*(t*^{0}*) dt*^{0}, (3.23)

for*〈u**i**〉 and defined in a similar way for 〈p〉 and 〈*S_{i j}*〉. The variables u*^{0}_{i}*, p*^{0} andS^{0}* _{i j}*
denote the fluctuating parts. This technique is called the Reynolds Decomposition
[15]. When we substitute these expressions into the Navier-Stokes equations from
chapter 2 and then average the equations, an additional term arises in the momen-
tum equations. This term contains the time average of the product of the velocity
fluctuations, known as the Reynolds Stress,

*R*_{i j}*= 〈u*^{0}_{i}*u*^{0}* _{j}*〉. (3.24)

Because of the new unknown Reynolds Stress there are not enough equations to determine all variables and therefore the problem is not closed anymore. Turbulence models provide closure to this problem by using additional equations based on exper- iments and assumptions to calculate the Reynolds Stress.

For compressible flow however, we have to take into account that the density is
variable as well. By applying the Reynolds Decomposition to the density and following
the same procedure as before we would obtain a more difficult problem due to the
arising correlation of*ρ*^{0} *and u*^{0}* _{i}* in the continuity equation and the triple correlation
of

*ρ*

^{0}

*, u*

^{0}

_{i}*and u*

^{0}

*in the momentum equation [18]. This problem is simplified by using mean velocity components based on the mass-averaged values opposed to the time-averaged values, also known as Favre-Averaging. In this case the velocity is expressed as,*

_{j}*u*_{i}*= 〈u**i*〉F*+ u*^{00}* _{i}*, (3.25)
where 〈 〉F denotes the Favre-average given by,

*〈u**i*〉*F* = 1
*ρ* lim

*T*→∞

Z *t**+T*
*t*

*ρ(x, τ)u**i**(x, τ)dτ,* (3.26)
and where*ρ denotes the time average. When we now substitute the Favre decom-*
positions of the velocity, specific internal energy, temperature and specific enthalpy,
together with the Reynolds decompositions of the pressure, density and heat flux
vector into the conservation equations, we obtain the Favre-averaged equations[18].

The Favre-averaged continuity and momentum equations differ from the laminar equations only by the Favre-Averaged Reynolds-Stress tensor,

*τ**i j* *= ρu*^{00}_{i}*u*^{00}* _{j}*, (3.27)

*where u*

^{00}

_{i}*and u*

^{00}

*denote the fluctuating parts of the Favre decomposition.*

_{j}**3.4.1** **The Realizable k-** **ε model**

**The Realizable k-**

**ε model**

*The realizable k-ε model uses the Turbulence Viscosity Hypothesis coined by Boussi-*
nesq in 1877, to provide closure to turbulence problem[15]. This hypothesis is anal-
ogous to the assumption of the viscous stress relation for a Newtonian fluid, which
we imposed earlier by employing equation 2.4. The Turbulence Viscosity Hypothesis
*states that the Reynolds-stress-anisotropy tensor a*_{i j}*defined as a*_{i j}*≡ 〈u**i**u** _{j}*〉 −

^{2}

_{3}

*kδ*

*i j*, is

*linearly related to the mean rate-of-strain tensor S*

*via the turbulent eddy viscosity*

_{i j}*ν*

*T*as,

*〈u**i**u** _{j}*〉 −2

3*kδ**i j* *= −ν**T*

*∂ 〈U**i*〉

*∂ x**j*

+*∂ 〈U**j*〉

*∂ x**i*

. (3.28)

*The realizable k-ε model then calculates the Reynolds stress by the relation,*
*ν**T* *= C*_{µ}*k*^{2}

*ε* , (3.29)

*where k andε represent the turbulent kinetic energy and the turbulent dissipation*
rate respectively, which are calculated from two transport equations for each time
*step. In contrary to the standard k-ε model where C** _{µ}* is a constant, in the Realizable

*k-ε model C*

*is a variable and defined by the expressions as listed in Table 3.2, where*

_{µ}*Ω*

*i j*is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular velocity

*ω*

*k*[16].

*C** _{µ}* =

_{A}^{1}

0*+A**s**kU∗*

*ε*

*Ω*˜*i j* *= Ω**i j**− 2ε**i jk**ω**k*

*U*^{∗} ≡q

*S*_{i j}*S*_{i j}*+ ˜Ω**i j**Ω*˜*i j* *Ω**i j* *= Ω**i j**− ε**i jk**ω**k*

Table 3.2: The expressions for *C** _{µ}*.

*Furthermore, the coefficients A*_{0} *and A** _{s}* are defined as in table 3.3.

*A*_{0} = 4.0 *W* = *S*_{i j}*S*_{jk}*S*_{ki}*S*˜^{3}
*A** _{s}* =p

6 cos*φ* *S*˜ *= ÆS**i j**S*_{i j}

*φ =* 1

3cos^{−1}(p

*6W) S**i j* = 1
2

*∂ u**j*

*∂ x**i*

+*∂ u**i*

*∂ x**j*

Table 3.3: The coefficients for*A*_{0} and *A** _{s}* in

*C*

*[16].*

_{µ}*The two additional equations of the k-ε model to close the problem are the*
*transport equations for the variables k andε [3],*

*∂*

*∂ t(ρk) +* *∂*

*∂ x**j*

*(ρku**j*) = *∂*

*∂ x**j*

*µ +* *µ**t*

*σ**k*

*∂ k*

*∂ x**j*

*+ P**k**+ P**b**− ρε − Y**M**+ S**k*, (3.30)

*∂*

*∂ t(ρε) +* *∂*

*∂ x**j*

*(ρεu**j*) = *∂*

*∂ x**j*

*µ +* *µ**t*

*σ*_{ε}

*∂ ε*

*∂ x**j*

*+ ρ C*1*Sε − ρ C*2

*ε*^{2}
*k*+*pνε*
*+ C*1*ε**ε*

*kC*_{3}_{ε}*P*_{b}*+ S** _{ε}*, (3.31)

where *µ**t* *= ρν**t**, P** _{k}* represents the generation of turbulence kinetic energy due to

*the mean velocity gradients, P*

*is the generation of turbulence kinetic energy due to*

_{b}*buoyancy, Y*

*represents the contribution of the fluctuating dilatation in compressible*

_{m}*turbulence to the overall dissipation rate, S*

_{k}*and S*

*are user-defined source terms and*

_{ε}*C*

*represents how much*

_{3e}*ε is affected by the buoyancy [10]. The other coefficients*for the transport equations 3.30 and 3.31 can be found in table 3.4.

*C*_{1} = max

0.43, *η*
*η + 5*

*σ** _{ε}* = 1.2

*η*

*= Æ2S*

*i j*

*S*

_{i j}*k*

*ε* *σ**k* = 1.0

*C*_{1ε} = 1.44 *C*_{2} = 1.9

Table 3.4: The coefficients for the transport equations of *k* and * _{ε}* [16].

**3.4.2** *k-* **ω SST**

**ω SST**

*The Shear Stress Transport (SST) k-ω turbulence model combines the k-ω model*
*with the k-ε model such that k-ω is used in the inner region of the boundary layer*
*and the k-ε in the free stream region. We will present the standard k-ω SST model as*
described by NASA in*[1], which is based on the article [13] by F.R. Menter in which*
he presented the model.

The turbulent viscosity is calculated as follows,

*ν**t*= *a*_{1}*k*

max*(a*1*ω, ΩF*2), (3.32)

where,

*Ω = Æ2W**i j**W** _{i j}*,

*W*

*= 1*

_{i j}2

*∂ u**i*

*∂ x**j*

−*∂ u**j*

*∂ x**i*

.

*The transport equations for k andω are,*

*∂ (ρk)*

*∂ t* +*∂ (ρu**j**k*)

*∂ x**j*

*= P − β*^{∗}*ρωk +* *∂*

*∂ x**j*

*(µ + σ**k**µ**t*)*∂ k*

*∂ x**j*

, (3.33)

*∂ (ρω)*

*∂ t* +*∂ (ρu**j**ω)*

*∂ x**j*

= *γ*
*ν**t*

*P− βρω*^{2}+ *∂*

*∂ x**j*

*(µ + σ*_{ω}*µ**t*)*∂ ω*

*∂ x**j*

+

2*(1 − F*1)*ρσ*_{ω2}*ω*

*∂ k*

*∂ x**j*

*∂ ω*

*∂ x**j*

. (3.34)

*With the following expression for P,*
*P* *= τ**i j*

*∂ u**i*

*∂ x**j*

, (3.35)

*τ**i j* *= µ**t*

*2S** _{i j}*−2
3

*∂ u**k*

*∂ x**k*

*δ**i j*

−2

3*ρkδ**i j*, (3.36)

*S** _{i j}* =1
2

*∂ u**j*

*∂ x**i*

+ *∂ u**i*

*∂ x**j*

. (3.37)

The remaining variable expressions and constants can be found in table 3.5.

*F*_{1} = tanh ((arg1)^{4}) *F*_{2} = tanh ((arg2)^{2})
arg_{1} = min

max

p

*k*

*β*^{∗}*ωd*,500*ν*
*d*^{2}*ω*

,4*ρσ*_{ω2}*k*
*C D*_{kω}*d*^{2}

arg_{2} = max

2

p*k*

*β*^{∗}*ωd*,500*ν*
*d*^{2}*ω*

*C D** _{kω}* = max

2*ρσ** _{ω2}* 1

*ω*

*∂ k*

*∂ x**j*

*∂ ω∂ x**j*

, 10^{−20}

*α*1 = 0.31

*k* = 0.41 *β*^{∗} = 0.09

*γ*1 = *β*1

*β*^{∗}−*σ*_{ω1}*k*^{2}

p*β*^{∗} *γ*2 = *β*2

*β*^{∗} −*σ*_{ω2}*k*^{2}
p*β*^{∗}

*σ**k1* = 0.85 *σ**k2* = 1.0

*σ** _{ω1}* = 0.5

*σ*

*= 0.856*

_{ω2}*β*1 = 0.075 *β*2 = 0.0828

Table 3.5: Expressions and constants for the standard SST *k*-* _{ω}*model
as presented by NASA [1].

The constants in table 3.5 that appear with subscript 1 and 2 are used in the equations
*with the blending function F*_{1}. For example for*β we have,*

*β = F*1*β*1*+ (1 − F*1*)β*2. (3.38)

**3.5** **Geometry and boundary conditions**

*y*

*x*

Figure 3.3: The geometry we meshed and its position within the bearing (not on scale).

Number of orifices Bearing clearance (µm) Orifice radius (µm) Bearing radius (mm)

8 4 50 25

Table 3.6: The specifications and dimensions of the modeled air bearing

In the simulations we have carried out, we considered a vertically and horizontally stationary bearing floating at a constant height of 4µm. As illustrated in Figure 3.3, the mesh surrounds only one-eighth of the bearing with dimensions as listed in table 3.6. Because of the symmetric shape of the bearing this is sufficient to determine the flow in the entire bearing clearance.

The mesh has four types of boundaries, namely the inlet at the top of the nozzle, the outlet at the edge of the bearing, the two symmetry planes on the clockwise and counterclockwise sides and the solid walls consisting of the surface above which the bearing floats and the bearing itself.

As suggested in most of the literature on pressure driven flow simulations, we
have used an initial flow consisting of a fluid at rest with a uniform pressure. During
the first 5× 10^{−5} s, we gradually increased the pressure difference by enhancing
the supply pressure to 6 bar and reducing the atmospheric back pressure to 1 bar.

By gradually introducing this pressure difference, it is more likely that the flow will develop stably. At first we used an initial uniform pressure of 1.7 bar, but later we changed it to 5 bar for two reasons. Firstly, based on results of simulations on the grooved bearing in earlier research[6], we expected that a uniform pressure of 5 bar

is closer to the final converged situation than a uniform pressure of 1.7 bar. Secondly, since most of the instabilities we encountered originated near the inlet, we expected it to be better to keep the inlet pressure more steady.

The majority of the simulations were executed with the boundary conditions as shown in table 3.7, but we have also tried others as described in Appendix A. For example for the inlet velocity, we have also used

### fixedValue

and### zeroGradient

as well, and for the inlet pressure we have tried### totalPressure

too.Variable Inlet Outlet Symmetry walls Solid walls

*u* pressureInlet-

Velocity

zeroGradient slip fixedValue (0,0,0)

*p* fixedValue

(Table)

fixedValue (Table)

zeroGradient zeroGradient

*T* fixedValue

(293)

zeroGradient zeroGradient zeroGradient

*k* fixedValue

(1× 10^{−6})

zeroGradient zeroGradient compressible::

kqRWallFunction
value: 1×10^{−12}

*ε* fixedValue

(1)

zeroGradient zeroGradient compressible::

epsilonWallFunction value: 500000

*ω* fixedValue

(1)

zeroGradient zeroGradient compressible::

omegaWallFunction
value: 500000
*α**t* zeroGradient zeroGradient zeroGradient compressible::

alphatWallFunction value: 0

*µ**t* zeroGradient zeroGradient zeroGradient compressible::

mutkWallFunction value: 0

Table 3.7: The boundary conditions for all patches.

**3.6** **Schemes and Transport model**

In order to perform a simulation, the user first has to define the thermophysical properties in the file

### thermophysicalProperties

. In listing 3.1 an overview of the settings that we have used for the thermophysical properties of the fluid can be found.t h e r m o T y p e {

20 t y p e h e P s i T h e r m o ; m i x t u r e p u r e M i x t u r e ; t r a n s p o r t s u t h e r l a n d ; s p e c i e s p e c i e ; t h e r m o e C o n s t ;

25 e q u a t i o n O f S t a t e p e r f e c t G a s ;

e n e r g y s e n s i b l e I n t e r n a l E n e r g y ; }

30 m i x t u r e {

s p e c i e {

n M o l e s 1;

35 m o l W e i g h t 2 8 . 9 6 ; }

t h e r m o d y n a m i c s {

40 Cp 1 0 0 4 . 5 ; // s p e c i f i c h ea t Cp

Hf 2 . 5 4 4 e + 0 6 ; // h e a t of f u s i o n Cv 7 1 7 ; // s p e c i f i c h e a t Cv }

45 t r a n s p o r t {

As 1 . 4 5 8 e - 0 6 ; // s u t h e r l a n d c o e f f i c i e n t Ts 1 1 0 . 4 ; // s u t h e r l a n d t e m p e r a t u r e }

50 }

Listing 3.1: thermophysicalProperties line 18 to 50.

The settings we used, imply that we relate the pressure to the temperature
through the equation of state of a perfect gas and that we model the thermodynamical
*behaviour by a constant specific heat c** _{p}*-model with evaluation of the internal energy

*eand entropy s without considering the heat of formation. Furthermore, we have*used the Sutherland transport model, which relates the dynamic viscosity

*µ and the*

*temperature T via the expression,*

*µ =* *A** _{s}*p

*T*1+

*T*

_{s}*T*

, (3.39)

*with the parameters A*_{s}*and T** _{s}* set as shown in listing 3.1.

**3.7** **Postprocessing**

To process the results of our simulations, we have used the open source data analysis and visualisation application paraView combined with the OpenFOAM reader module as supplied by OpenFOAM. With paraView we could for example investigate the pressure distribution and velocity profiles throughout the entire mesh. With these tools we inspected whether the results we received were realistic and whether they were converging or not. We also used it extensively to determine our strategy to improve future simulations.

**Chapter 4** **Results**

**4.1** **Adapting the mesh**

**The first Mesh**

(a) A close-up of the nozzle. (b) A close-up of the cross-section of the
bearing with the *y, z*-plane.

Figure 4.1: The first mesh of the bearing with eight orifices.

For the first simulation, we used a mesh with a relatively uniform cell distribution, containing 814.464 cells of which 99.62% is hexahedral. Figure 4.1 contains two pictures of the mesh. The nozzle is octagonal instead of round, to make the meshing procedure easier.

For this simulation we used the settings that led to good results in the past when Schut performed simulations of a grooved air bearing. For those simulations they used the solver

### sonicFoam

*together with the realizable k− ε turbulence model. The*grooved bearing that they modeled is comparable to our bearing in dimensions, but the design is different. It has only three orifices instead of eight, it contains recesses, also known as pockets, around the orifices and it has anchor shaped distribution

grooves from the orifices to about three-quarter of the bearing, as can be seen in Figure 4.2. Furthermore, in their simulations, only bearing clearances of 9µm and larger were considered.

(a) A close-up of the air chamber above the nozzle and the recess around the orifice.

(b) An overview of the en- tire mesh.

Figure 4.2: The mesh of one-third of the grooved bearing, used by the company in earlier research.

The first simulation on the first mesh with the settings as described above crashed
*at t*= 201.60×10^{−5}s with a final maximum velocity of 635 m s^{−1}and a final maximum
pressure of 17.4 bar. In Figure 4.3 it can be seen that there is a pressure accumulation
moving downward through the nozzle, which then collapses back towards the inlet
after which the simulation crashed. The pressure in the nozzle reaches values of over
18 bar. Beyond the orifice however, very little happens. Only in the region near the

outlet a little activity can be seen.

#ID Solver Turbulence Model Progress Latest Time (×10^{−5}s) *v*_{max}(m s^{−1}) *p*_{max}(bar)

5 sonicFoam *Realizable k**− ε* Crash 201.60 635 17.4

Table 4.1: The results of the first simulation on the first mesh.

Figure 4.3: The pressure distribution over time in the *y, z*-plane. Simu-
lation #5.

**Second Mesh**

The second mesh that we used is a finer mesh, see Figure 4.4. This mesh contains 7.756.848 cells, of which 99.85% is hexahedral. This mesh has a lot more cells in the nozzle than the previous mesh, which can be seen by comparing Figure 4.1a with Figure 4.4a.

(a) A close-up of the top of the nozzle. (b) A close-up of a cross-section of the
mesh with the *y, z*-plane.

Figure 4.4: The second mesh of the bearing with eight orifices.

The two simulations that we ran with the solver

### rhoSimpleFoam

both show apressure distribution with multiple stationary circular pressure depressions around the orifice. These depressions range from strong to invisible as you move away from the orifice, as can be seen in the top-view in Figure 4.5a and in a cross-section in Figure 4.5b.

(a) A top-view of the nozzle and the sur- rounding pressure depressions.

(b) The pressure in the *y, z*-plane.

(c) The normalised velocity vectors in the transition from the orifice to the gap.

(d) The velocity magnitude in the inlet.

Figure 4.5: The results of the second mesh. Simulation #6.

Furthermore, the flow exits the orifice in multiple strong streams between which the flow is reversed, as indicated by the yellow arrows that represent the normalised velocity vectors in Figure 4.6.

The results of the simulation with

### sonicFoam

barely show any movement nor change in pressure, except for in small regions near the inlet and outlet. In Fig- ure 4.7 and Figure 4.8 one can see that the fluid underneath the bearing is nearly stationary and that the pressure remains at the initial value of 1.7 bar almost every- where. Because this did not seem to change with time we pauzed simulation #8 at*t*= 6.3 × 10^{−5} s.

Figure 4.6: A top-view of the velocity profile. The yellow arrows repre- sent the normalised velocity vectors. Simulation #6.

#ID Solver Turbulence Model Progress Latest Time (×10^{−5}s) *v*_{max}(m s^{−1}) *p*_{max}(bar)

6 rhoSimpleFoam *Realizable k**− ε* Crash 1.6 702 2.14

7 rhoSimpleFoam Laminar Crash 1.5 646 2.06

8 sonicFoam *Realizable k**− ε* Pauzed 6.3 365 6.02

Table 4.2: The results of the second mesh.

(a) The pressure distribution in a cross-
section with the *y, z*-plane.

(b) A top view of the pressure distribution close to the outlet.

Figure 4.7: The results of the second mesh. Simulation #8.

(a) The velocity in a cross-section with the
*y, z*-plane.

(b) A top view of the velocity close to the outlet.

Figure 4.8: The results of the second mesh. Simulation #8.

**Third Mesh**

The final mesh that we have used for all simulations on this bearing from this point on, contains 1.150.044 cells of which 99.74% is hexahedral. Because our geometry contains a lot of thin regions, we could reduce the amount of cells without losing a lot of accuracy. This is done by elongating the cells in the flow direction in regions where the flow is mostly laminar. To do so, we have used four gradients as visible in Figure 4.9. Firstly, there is a vertical gradient in the nozzle, since the flow in the top of the nozzle is irrelevant to the company compared to the flow in the bearing clearance.

(a)A close-up of the nozzle. (b) A cross-section of the nozzle with the
*y, z*-plane.

Figure 4.9: The third mesh of the bearing with eight orifices.

*Secondly, there are gradients in the bearing clearance in the x and y direction,*
because we expect the flow to laminarize as it gets farther away from the orifice.

*Lastly, there is a vertical gradient in the z direction in the bearing clearance to make*

(a) The top view of the final pressure dis- tribution.

(b) The pressure versus the *y* coordinate.

Figure 4.10: The pressure in simulation #10.

(a) The normalised velocity vectors in the transition from the orifice to the gap.

(b) The velocity magnitude in the inlet.

Figure 4.11: The velocity in simulation #10.

the cells near the two solid walls more fine. Furthermore, this mesh has a square orifice to make the meshing procedure easier.

We executed two simulations of which the behaviour is almost similar to the results of the simulations we performed with the second mesh. Similarly to simulation #8, the simulation with

### sonicFoam

barely initiated a flow so we pauzed the simulation.The simulation with