FOAM1.7performdifferentlyfromsonicFoamandrhoSimpleFoaminOpenFOAM airbearings,treatthefluiddynamicsthatisinvolvedandconsiderthepreviously softwarepackageOpenFOAM.Thetypeofairbearingconsideredinthisresearch

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 tensorTin 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

2grad 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²− 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_sand 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_ato p_son 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, 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 ¹₈ 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ⁿ⁻¹, ρⁿ⁻¹ and pⁿ⁻¹

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⁰

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 nth 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ⁿ⁻¹ andρⁿ⁻¹ 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^m−1 and ρ^m−1 from the previous outer iteration are used instead. These 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ⁱu^m∗_i,P+X

l

A^u_lⁱu^m∗_i,l = Q^m−1_u

i − ∆Ω

δp^m−1 δxi



. (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ⁱu^m∗_i,l A^u_Pⁱ −∆Ω

A^u_Pⁱ

δp^m−1 δxi



. (3.2)

Because the velocity field obtained from equation 3.2 together with the old density valuesρⁿ⁻¹ do not satisfy continuity, substituting the new interpolated cell face mass fluxes ˙m^∗_l into the continuity equation,

(ρ^m−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ρ⁰and the normal velocity correction v_n⁰ such that the corrected mass flux for the mth outer iteration and lth cell face can be expressed as,

˙

m^m_l = ρ_l^mv_n,l^m S_l= (ρ^m−1+ ρ⁰)l(v_n^m∗+ v_n⁰)lS_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⁰_l = (ρ^m−1S v_n⁰)l+ (v_n^m∗Sρ⁰)l+ (ρ⁰v_n⁰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−1S v_n⁰)l= − (ρ^m−1S∆Ω)l

 1 A^u_P



l

δp⁰ δ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ρ⁰in the second term of equation 3.5 by

ρ⁰=

∂ ρ

∂ p

‹

T

p⁰= C_ρp⁰. (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ρ⁰)l = C_ρm˙^∗

ρ^m−1



l

p⁰_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⁰_l = − (ρ^m−1S∆Ω)l

 1 A^u_P



l

δp⁰ δn

‹

l

+ C_ρm˙^∗ ρ^m−1



l

p⁰_l. (3.10)

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

ρ⁰_P∆Ω

∆t +X

l

˙

m⁰_l+ Q^∗_m= 0. (3.11)

If we substitute equation 3.7 forρ⁰_P and equation 3.10 for ˙m⁰_l into equation 3.11, we can write the pressure correction equation as,

A_Pp⁰_P+X

l

A_lp⁰_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+ αPp⁰. (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⁰into equation 3.2 gives,

u⁰_i,P= − P

l A^u_lⁱu⁰_i,l

A^u_Pⁱ −∆Ω A^u_Pⁱ

δp⁰ δxi

‹

P

, (3.14)

from which the corrected velocity u^mis 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ⁿ⁻¹, ρⁿ⁻¹ and pⁿ⁻¹

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⁰

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 nth 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ρⁿ⁻¹, this estimate u^m∗ is then used to calculate an estimate of the mass fluxes through the cell faces ˙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⁰, a correction on the densityρ⁰and a correction on the cell face normal velocity v_n⁰, 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 +Sf and−Sf 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± = Sf · uf±

p_f± = ρf±RT_f± c_f± = ÆγRTf±

T_exp = µ[(∇u)^T −²₃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) − ∇ · (Texp) = 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



, (3.18)

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

∂ (ρcvT)

∂ 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 = 〈ui〉 + u⁰_i, (3.20)

p= 〈p〉 + p⁰, (3.21)

S_{i j} = 〈S_{i j〉 +}S⁰_{i j}. (3.22) Here〈ui〉, 〈p〉 and 〈S_{i j}〉 denote the time averages defined as,

〈ui〉 = 1 T

Z t+T t

u_i(t⁰) dt⁰, (3.23)

for〈ui〉 and defined in a similar way for 〈p〉 and 〈S_{i j}〉. The variables u⁰_i, p⁰ andS⁰_{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⁰_iu⁰_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ρ⁰ and u⁰_i in the continuity equation and the triple correlation of ρ⁰, u⁰_i and u⁰_j 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,

u_i = 〈ui〉F+ u⁰⁰_i, (3.25) where 〈 〉F denotes the Favre-average given by,

〈ui〉F = 1 ρ lim

T→∞

Z t+T t

ρ(x, τ)ui(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⁰⁰_iu⁰⁰_j, (3.27) where u⁰⁰_i and u⁰⁰_j denote the fluctuating parts of the Favre decomposition.

3.4.1 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} ≡ 〈uiu_j〉 −²₃kδi j, is linearly related to the mean rate-of-strain tensor S_{i j} via the turbulent eddy viscosity νT as,

〈uiu_j〉 −2

3kδi j = −νT

∂ 〈Ui〉

∂ xj

+∂ 〈Uj〉

∂ xi



. (3.28)

The realizable k-ε model then calculates the Reynolds stress by the relation, νT = C_µk²

ε , (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 ¹

0+AskU∗

ε

Ω˜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₀ and A_s are defined as in table 3.3.

A₀ = 4.0 W = S_{i j}S_jkS_ki S˜³ A_s =p

6 cosφ S˜ = ÆSi jS_{i j}

φ = 1

3cos⁻¹(p

6W) Si j = 1 2

∂ uj

∂ xi

+∂ ui

∂ xj



Table 3.3: The coefficients forA₀ 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) + ∂

∂ xj

(ρkuj) = ∂

∂ xj



µ + µt

σk

‹ ∂ k

∂ xj



+ Pk+ Pb− ρε − YM+ Sk, (3.30)

∂

∂ t(ρε) + ∂

∂ xj

(ρεuj) = ∂

∂ xj



µ + µt

σ_ε

‹ ∂ ε

∂ xj



+ ρ C1Sε − ρ C2

ε² k+pνε + C1εε

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_b is the generation of turbulence kinetic energy due to buoyancy, Y_m represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, S_kand S_εare user-defined source terms and C_3erepresents how muchε 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₁ = max

•

0.43, η η + 5

˜

σ_ε = 1.2 η = Æ2Si jS_{i j}k

ε σk = 1.0

C_1ε = 1.44 C₂ = 1.9

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

3.4.2 k- ω 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₁k

max(a1ω, ΩF2), (3.32)

where,

Ω = Æ2Wi jW_{i j}, W_{i j}= 1

2

∂ ui

∂ xj

−∂ uj

∂ xi

 .

The transport equations for k andω are,

∂ (ρk)

∂ t +∂ (ρujk)

∂ xj

= P − β^∗ρωk + ∂

∂ xj



(µ + σkµt)∂ k

∂ xj



, (3.33)

∂ (ρω)

∂ t +∂ (ρujω)

∂ xj

= γ νt

P− βρω²+ ∂

∂ xj



(µ + σ_ωµt)∂ ω

∂ xj

 +

2(1 − F1)ρσ_ω2 ω

∂ k

∂ xj

∂ ω

∂ xj

. (3.34)

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

∂ ui

∂ xj

, (3.35)

τi j = µt



2S_{i j}−2 3

∂ uk

∂ xk

δi j

‹

−2

3ρkδi j, (3.36)

S_{i j} =1 2

∂ uj

∂ xi

+ ∂ ui

∂ xj



. (3.37)

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

F₁ = tanh ((arg1)⁴) F₂ = tanh ((arg2)²) arg₁ = min

– max

‚ p

k

β^∗ωd,500ν d²ω

Œ

,4ρσ_ω2k C D_kωd²

™

arg₂ = max

‚ 2

pk

β^∗ωd,500ν d²ω

Œ

C D_kω = max



2ρσ_ω2 1 ω ∂ k

∂ xj

∂ ω∂ xj

, 10⁻²⁰

‹

α1 = 0.31

k = 0.41 β^∗ = 0.09

γ1 = β1

β^∗−σ_ω1k²

pβ^∗ γ2 = β2

β^∗ −σ_ω2k² pβ^∗

σk1 = 0.85 σk2 = 1.0

σ_ω1 = 0.5 σ_ω2 = 0.856

β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₁. For example forβ we have,

β = F1β1+ (1 − F1)β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⁻⁵ 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⁻⁶)

zeroGradient zeroGradient compressible::

kqRWallFunction value: 1×10⁻¹²

ε 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_sp 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⁻⁵s with a final maximum velocity of 635 m s⁻¹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⁻⁵s) v_max(m s⁻¹) 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 a

pressure 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⁻⁵ 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⁻⁵s) v_max(m s⁻¹) 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

rhoSimpleFoam

showed pressure depressions around the orifice again. This time however, the pressure depressions were not static but they were moving away from the orifice, similarly to a sound wave. Furthermore, the depressions were divided into four directions, and there were more visible depressions than in the previous simulations, as illustrated in Figure 4.10a. Lastly, the velocity profile in the nozzle showed a higher velocity near the walls and corners, as can be seen in Figure 4.11b.