R Gourdon

M Jean

Marc Lancelin

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

H Charles

Gudrun Mme

Bornette

Safia Ait

Chalal Insa

M Lagarde

M Didier

M Johannes Kellendonk

M Laboune

Jean-Yves Buffiere
email: yves.buffiere@insa-lyon.fr

Mega

M Philippe Boisse

Scso Scso

M Obadia

Lionel Sec

Viviane Polsinelli

J Y Toussaint

This study is the result of collaborative work between the LaMCoS laboratory (Lyon) and the LIEBHERR AEROSPACE company (Toulouse).

Résumé

Introduction

La réduction des masses et volumes embarqués est une question centrale dans le domaine des transports, en particulier dans le secteur aéronautique. Il est à l'heure actuelle possible d'améliorer les performances des turbomachines, notamment en augmentant les vitesses de rotation et en diminuant le diamètre des arbres. De cette fac ¸on, il est alors possible de concevoir des turbomachines qui seront à la fois plus compactes et plus efficaces.

Les machines à cycle air font parties des turbomachines à haute vitesse et sont un des éléments centraux du système de contrôle de l'environnement de l'avion (Figure 1). Ce système est celui qui est responsable du contrôle du refroidissement, du chauffage et de la mise en pression au sein de l'appareil. De nos jours, la plupart des systèmes de contrôle de l'environnement installés sur les avions civils et militaires utilisent des paliers à feuilles.

FIGURE 1: Air Conditioning Pack

La solution la plus appropriée au fonctionnement de ces turbomachines dans un environnement réfrigerant consiste en la lubrification des paliers à feuilles directement par ce même gaz (Figure 2). Une approche similaire [HOW 06, HOW 08a] a été étudiée avec des gaz inertes à haute pression. Toutefois, certains problèmes persistent lors de la mise en place de paliers à feuilles au sein de nouveaux systèmes, notamment lorsqu'il s'agit d'un milieu réfrigerant. Les paliers à feuilles ont montré de nombreuses qualités. Ils peuvent fonctionner à des vitesses de rotations très élevées, ne nécessitent que peu de maintenance, ont de bonnes propriétés dynamiques, sont compactes et n'ont pas de système de lubrification externe. En comparaison avec des paliers lubrifiés à l'huile, les paliers à feuilles sont capables de travailler à des températures plus extrêmes. (Depuis les milieux cryogéniques jusqu'à 700 ˚C [DEL 97]) En revanche, la mauvaise gestion des aspects thermiques des paliers à feuilles est l'un des deux mécanismes de base qui limitent les performances des paliers à feuilles [ DYK 06].

De nos jours, la gestion des aspects thermiques est toujours un aspect central dans l'étude des paliers, considérée comment étant aussi important que la stabilité dynamique du rotor ou encore la capacité de charge [LEE 10]. La génération de chaleur dans le fluide et l'élévation de température associée devient conséquente lorsque la vitesse de rotation est élevée car les pertes de puissance sont en majorité transformées en chaleur.[ HOR 08] Une élévation de la temperature moyenne du palier à feuille ou encore une forte surpression locale peuvent être néfastes pour les performances du palier, dégradant notamment la capacité de charge, l'épaisseur de film minimale, etc. [MA 96] Cette élévation de température modifie également certaines propriétés du lubrifiant (telles que la viscosité moléculaire). Ceci peut entrainer des déformations au sein même du palier et de son fourreau. De plus, dans le cas particulier des paliers à feuilles, la pression de vapeur saturante dépend fortement de la température locale du lubrifiant.

De récentes études concernant les paliers à feuilles qui utilisaient des modèles THD 3D, et pour lesquelles les prédictions étaient en accord avec les données expérimentales, ont montrées que les gradients de pression axiaux augmentaient avec la vitesse de rotation de l'arbre et qu'ils ne pouvaient être négligés à très haute vitesse : l'utilisation d'approximation par des modèles 2D entraine la perte d'informations en rapport avec les bords du palier mais également une surrestimation des valeurs maximales des paramètres [MIC 07]. Dobrica et Fillon [DOB 06] ont confirmé que les algorithmes 3D étaient indispensables à une bonne prédiction des phénomènes thermiques dans la direction axiale. De plus, l'étude des déformations élastiques axiales nécessite la modélisation des variations de température dans les trois directions du palier [FIL 08, MAO 08, SIN 08]. L'étude de l'utilisation des gaz réfrigérants jouant le rôle de lubrifiant est justement l'un des cas où la modélisation 3D aide à déterminer l'existence éventuelle d'une phase liquide.

Théorie

Equation d'état

L'hypothèse de base de la théorie des gaz parfaits est que la densité est proportionnelle à la pression, ces deux paramètres étant reliés par le biais d'un coefficient qui est lui même dépendent des propriétés intrinsèques du gaz ainsi que de sa température. Cependant, lorsque la pression est suffisamment proche de la pression de vapeur saturante, qui est la pression pour laquelle la transition liquide/vapeur se produit, cette hypothèse n'est plus valable. Un écoulement diphasique a été observé expérimentalement dans des conditions spécifiques, à savoir quand du réfrigérant est introduit dans le système de lubrification d'un palier à feuilles. Dans ces conditions, la loi des gaz parfait n'est plus valable.

En conséquence, nous utilisons une équation d'état capable de décrire la variation de densité en fonction de la pression et de la température, mais également la transition liquide/vapeur. Si de nombreuses équations d'état sont capables de décrire ce comportement non linéaire, peu d'entre elles satisfont nos critères. Premièrement, l'équation d'état doit pouvoir décrire le comportement du fluide avec suffisamment de précision pour ne pas nuire à la celle du modèle dans son ensemble. De plus, cette équation d'état n'est qu'une partie d'un modèle beaucoup plus complexe (THD) et se doit d'être relativement simple afin de ne pas trop complexifier le modèle.

En fait, la densité et la pression sont fortement couplées et l'équation de Reynolds iii généralisée [DOW 62] ainsi que l'équation d'état doivent être résolues simultanément. Enfin, ce modèle est destiné à être adapté à différents types de lubrifiant et doit donc être facilement modifiable.

C'est en prenant en compte l'ensemble de ces critères que nous avons choisi le type d'équation d'état le mieux adapté à ce modèle, un modèle cubique qui consiste en une équation polynomiale du troisième degré où la densité joue le rôle d'inconnue et où la pression est l'un des paramètres de l'équation. Nous avons choisi l' équation modifiée de ].

La comparaison de ce modèle avec la base de données thermodynamiques de référence "REFPROP" (NIST) montre que la pression et la densié du fluide (qui dépend de la température) peuvent être prédits grâce à cette équation d'état avec une erreur inférieure à 2% dans la phase vapeur et moins de 5% d'erreur dans la phase liquide tant que le fluide n'est pas dans des conditions proches de celles de la transition liquide/vapeur. La seule imprécision provient des résultats donnés à proximité de la zone de transition. Afin d'éliminer cette imprécision, nous utilisons la formule de Clapeyron pour finalement prédire la pression de vapeur saturante avec une erreur inférieure à 1%. Finalement, l'erreur maximum concernant la valeur de la densité chute donc à 2% dans la phase vapeur et à 5% pour la phase liquide.

La formule de Clapeyron décrit le comportement thermodynamique du fluide, et plus précisément la valeur de la pression de transition de l'état vapeur à l'état liquide. Grâce à l'approximation de Dupré, seuls trois couples (température/pression de vapeur saturante associée) sont nécessaires afin de prédire la pression de vapeur saturante à n'importe quelle température.

En conclusion, cette équation d'état peut être résolue en connaissant les paramètres caractéristiques suivants (qui sont des données entrées par l'utililsateur) : masse molaire, température critique, pression critique, et les trois couples (température/pression de vapeur saturante associée).

Transition liquide/vapeur

Nous avons également intégré un modèle de mélange liquide/vapeur à notre modèle [ ODY 03]. Une transition sinusodale, fonction de la pression locale, de la pression de vapeur saturante et de la vitesse du son minimale dans le mélange permet de décrire les valeurs des paramètres équivalent du mélange. Grâce à ce modèle nous pouvons calculer un domaine de pressions allant de la pression de vapeur saturante jusqu'à sa limite supérieure (calculée) où les paramètres vont être ceux du mélange (densité, viscosité, conductivité thermique) avec une transition continue de l'état purement vapeur à l'état purement liquide.

Equation de Reynolds généralisée

Le champ de pression dans les paliers à feuilles s'obtient en résolvant l'équation de Reynolds généralisée en régime permanant, pour des fluides turbulents, compressibles et avec une viscosité qui varie à travers l'épaisseur de film. Les conditions aux limites imposées pour la pression sont : la valeur de la pression connue à chaque extrémité du palier et au niveau de la gorge d'alimentation éventuelle. Concernant la transition liquide/vapeur, il y a une variation des valeurs physiques telles que la densité et la viscosité mais notre modèle assure la continuité de ces valeurs.

De plus, nous tenons compte des effets de la turbulence dans cette équation de Reynolds généralisée par l'intermédiaire de la viscosité équivalente qui est la somme des viscosités moléculaire et tourbillonnaire (due à l'écoulement turbulent). Ici aussi, si les conditions de fonctionnement conduisent à un écoulement laminaire du fluid, alors la viscosité tourbillonnaire est égale à zéro et la viscosité équivalente devient égale à la seule viscosité moléculaire.

Equation de l'énergie

On résout l'équation de l'énergie en film mince tridimensionnelle en régime permanent [FRE 90] afin d'obtenir le champ de température et de viscosité moléculaire associée. Ce modèle tient compte des variations de température dans les trois directions de l'écoulement : circonférentielle, axiale et radiale. La direction radiale (celle de l'épaisseur de film) est fondamentale pour comprendre le processus de génération de chaleur au sein du film lubrifiant. Ce sont les contraintes de cisaillement selon cette direction qui sont responsables de la majeure partie des pertes par dissipation.

Le gradient axial de pression a également un effet sur le champ de température : quand le gradient de pression est négatif, une partie du lubrifiant est aspiré par les côtés du palier et contribue ainsi au processus de refroidissement. Toutefois, le gradient de temperature axial est néfaste puisqu'il engendre des contraintes supplémentaires dans la structure du palier ce qui est potentiellement dangereux, notamment dans le cas des feuilles de métal des paliers à feuilles.

Modèle de turbulence

Nous avons choisi le modèle de viscosité tourbillonnaire tridimensionnelle dans lequel l'impact de la turbulence est calculé grâce à une version modifiée de la formule empirique de Reichardt. Les coefficients empiriques sont basés sur les mesures des contraintes de cisaillement et des vitesses dans un conduit [REI 51] et adaptés aux données expérimentales [NG 64] v Ce modèle donne de bons résultats prédictifs pour les paliers et fourni un modèle 3D de la viscosité tourbillonnaire. En cela, il tient compte de la turbulence dans les équations de Reynolds et de l'énergie.

Méthodes numériques

La théorie de la méthode des différences finies est basée sur des principes simples mais elle est capable de résoudre des problèmes THD complexes avec précision et efficacité. De ce fait, de nombreuses études récentes utilisent cette méthode afin de résoudre des problèmes hydrodynamiques ou THD [CHA 02, BOU 03, COS 03, BRU 05, AHM 10]. De plus, la méthode des différences finies est une méthode très efficace pour des géométries simples telles que celle des paliers à feuilles.

L'équation d'état est résolue en utilisant un algorithme de Gauss-Seidel et une méthode de discrétisation par différences finies centrées. Les tests que nous avons effectués ont montré que la relaxation pouvait être utilisée pour avoir une accélération de la convergence dans certains cas simples ou bien pour raffiner la résolution dans des cas plus critiques (Par exemple, une forte excentration de l'arbre combinée à une vitesse élevée). L'erreur relative locale sur la pression est calculée en chaque point et à chaque itération. L'équation de Reynolds généralisée est considérée résolue lorsque l'erreur maximum est inférieure à un critère prédéfini. L'équation de l'énergie en film mince est résolue en utilisant un algorithme de Gauss-Seidel et une méthode de discrétisation par différences finies décentrées. Afin de prendre en compte l'orientation de l'écoulement au sein du palier, il est nécessaire d'utiliser un schéma décentré pour obtenir la convergence dans le cas où des dérivées de la température apparaissent, selon que le noeud en un point va être influencé par tel ou tel voisin. Afin de s'affranchir d'un certain manque de précision de l'algorithme en certains points isolés (notamment lorsque des recirculations apparaissent), nous utilisons un double critère de convergence qui combine à la fois les erreurs locales et globales. Cet algorithme est suffisamment robuste pour donner la solution de l'équation de l'énergie en film mince à de très hautes vitesses (100 000 tr/min) et des excentricités relatives élevées (0.9) ce qui, combinées, sont des conditions de fonctionnement très sévères comparées aux applications industrielles actuelles.

Résultats et discussions

Contexte et hypothèses

Nous avons étudié les performances d'un type de palier à feuilles ayant des propriétés et fonctionnant dans des conditions à même de reproduire un cas diphasique. L'un des objectifs principaux de cette étude est de décrire le comportement des pavi ramètres locaux du fluide et la manière dont ils sont couplés. C'est la raison pour laquelle nous avons dans cette partie fait des hypothèses concernant le lubrifiant à la fois simples mais réalistes. Cette étape est indispensable pour obtenir une bonne compréhension des phénomènes complexes qui ont lieu au sein du palier. Parmi les hypothèses principales que nous avons effectuées, nous avons :

-Un palier rigide. Il s'agit d'une hypothèse forte pour les paliers à feuilles, mais même si notre modèle est capable de prendre en compte les déformations des feuilles, ce phénomène interagit fortement avec le fluide et ne peut donc être pris en compte dans cette première étude.

-Une température du rotor imposée. Nous avons limité cette étude à l'influence de l'écoulement diphasique sur les paramètres du fluide puisque notre but est d'analyser le couplage entre la pression, la densité, la viscosité. Afin d'obtenir une bonne compréhension de la plupart des intéractions, l'excentricité relative est fixée à 0.9 (ce qui est une condition de fonctionnement sévère). La vitesse varie de 20 000 à 100 000 tr/mn ce qui, une fois de plus, est une limite élevée comparée à la plupart des applications industrielles. Cette étendue de l'étude permet une bonne compréhension générale du comportement du fluide.

Etude de la pression

Concernant la pression à faible vitesse (20 000 tr/mn), nous retrouvons un profil classique avec un pic de pression aux alentours de 180 ˚et un très léger effet due à l'écoulement diphasique. Toutefois, à des vitesses de rotation plus élevées, les effets combinés de la turbulence et de l'écoulement diphasique conduisent à un profil de pression inhabituel.

Les effets de la turbulence sont pris en compte dans l'équation de Reynolds grâce à la viscosité équivalente qui a tendance à accrotre la pression dans des conditions données, comparé à un écoulement turbulent. Les effets de l'écoulement diphasique sont plus subtiles. Puisque la pression de vapeur saturante est donnée par la température locale, la transition est déterminée de fac ¸on claire. Dans la zone de mélange, le pic de pression est plus faible que celui obtenu dans un écoulement monophasique (dans des conditions vii équivalentes) puisque, dans la zone de mélange, la pression a tendance à être proche de la pression de vapeur saturante tandis que la densité augmente fortement.

Il faut également noter que la localisation circonférentielle du maximum de pression se décale dans le sens de la rotation quand la vitesse augmente.

Etude de la densité

La turbulence n'a pas d'impact direct sur la densité, puisque le terme de viscosité n'apparait pas dans l'équation d'état. Mais puisque la densité est fonction de la pression, les paramètres de la turbulence et la densité sont également fortement couplés. Le couplage direct avec la pression est clairement visible sur le profil de densité. Dans une zone autour du minimum d'épaisseur de film, qui s'étend jusqu'à un sixième de la circonférence environ, la pression est au dessus de la pression de vapeur saturante et le mélange se forme.

Une densité équivalente prend en compte la proportion de liquide et, de ce fait, génère un important pic au niveau du minimum d'épaisseur de film. La densité est également fonction de la température puisque la pression de vapeur saturante est fortement influencée par la température. C'est pourquoi un modèle THD 3D est nécessaire dans le cas de l'étude d'écoulements diphasiques.

Etude de la viscosité

Les variations de viscosité moléculaire dues à la température ou à la densité sont du même ordre de grandeur (ou moins) que sa valeur initiale, calculée dans l'environnement du palier. Comparé à la viscosité moléculaire, la viscosité tourbillonnaire peut augmenter jusqu'à un ordre de grandeur supérieur. Nous en concluons que la viscosité équivalente en régime turbulent est principalement déterminée par la viscosité tourbillonnaire. Cependant, il est important de calculer la viscosité dynamique avec précision puisque ses plus importantes variations sont situées aux endroits critiques (la zone de mélange autour du minimum d'épaisseur de film).

Etude de la température

Le champ de température est tridimensionnel mais la plupart des effets sont visibles sur le profil de température à mi-longueur. Les résultats suivants sont ceux d'un test à 100 000 tr/mn. La même tendance est observée pour toutes les vitesses de rotation. De 0 ˚à 60 ˚une augmentation progressive de la température est constatée, jusqu'à ce qu'une discontinuité dans la pente survienne, du fait du courant inverse généré par le gradient de pression. (A proximité de la feuille supérieure, l'écoulement de Poiseuille domine largement celui de Couette).

Ensuite, le profile de température augmente depuis cette position jusqu'au début de viii la zone de mélange de manière standard. Dans la zone de mélange (de 140 ˚à 210 ˚) il y a une baisse de la température due à la proportion de liquide que le fluide contient. De 210 ˚à 270 ˚un gradient de pression axial permet au lubrifiant provenant de l'extérieur du palier de pénétrer à l'intérieur de ce dernier, générant ainsi une baisse soudaine de la température. Enfin, de 270 ˚à 360 ˚le processus de génération de chaleur due aux contraintes de cisaillement continue, alors que le gradient de pression axial diminue.

Conclusions

Au cours de cette étude, nous avons montré l'importance d'une description précise des paramètres du film, dont les variations influencent le comportement du palier. Parmi les principales théories qui ont été utilisées ici, nous pouvons citer :

-Un lubrifiant compressible, décrit par une loi de comportement non-linéaire appropriée, au voisinage de la transition liquide/vapeur.

-La transition liquide/vapeur et le calcul des paramètres équivalents du mélange.

-Un écoulement turbulent pour les paliers à feuilles à haute vitesse, avec un modèle de viscosité tourbillonnaire.

-Un comportement tridimensionnel pour la viscosité, en particulier concernant les variations selon l'épaisseur de film. (dépendant de la température)

-Un comportement tridimensionnel de la température, en particulier concernant les variations selon l'épaisseur de film afin d'être cohérent avec le modèle de viscosité (moléculaire et tourbillonnaire), mais également dans la direction axiale, afin de prendre en compte les gradients de température potentiels, qui modifient considérablement le champ de température tridimensionnel du palier.

ix

Table des matières

Introduction

The reduction of on-board masses and volumes is a central issue regarding transports, particularly in aeronautics. It means either more passengers or lower energetic cost, which is a major issue in regards to the increase in oil prices. Improvements in turbomachinery are possible, thanks to higher rotation speeds and smaller shaft diameters. As a result, turbomachines can be more compact and more efficient at the same time [LEL 07].

The Environmental Control System (ECS) is the structure in an aircraft that manages cooling, heating and pressurizing the air using an air conditioning pack (Figure 3). The main process consists in bleeding the air from the engines to cool it to the suitable pressure and temperature. Then this air is distributed to the various cabin zones at the right temperature and pressure to manage passenger's comfort and aircraft safety. (Figure 4) The purpose of the air conditioning pack is to decrease the temperature and remove the water contained in the hot air bled from the engines. The main target of the operation carried out by the air cooling pack is to put the air into a device where its volume will be increased and accordingly its temperature will decrease. The following explanations describe how this is achieved : Two heat exchangers decrease, in two steps, the temperature of this air. Between theses two steps, an Air Cycle Machine (ACM) (Figure 5 and Figure 6) in its compressor stage compresses the air. Then in its turbine stage, the air expands thus lowering the air temperature to a very low level. To avoid ice built up, some operations are added to the previous steps described here-above to remove water. This is done in a condenser which condenses the water vapor included in the air and a water extractor which removes this water outside. (Gross, Barnett and Sylver that would become Honeywell) and it has made major progress since then. Reliability of many high-speed turbo-machines has been improved over tenfold compared to former turbo-machines which were using rolling element bearings. The use of GFBs in air-system turbo-machinery for civilian aeronautics started at the end of the 1970's [START_REF]The Application of Foil Air Bearing Turbomachinery in Aircraft Environmental Control Systems[END_REF]. Moreover, using GFBs in micro gas turbines reduces initial and running costs [OHK 05].

Introduction

The following paragraph is a general description of the main components of a GFB (Figure 7) : The housing is supposed to be a rigid component where the bearing itself is located. The bearing is made up of several compliant surfaces (bumps and top foil) and a high-speed rotating shaft. The stiffness and damping of the bump-strip layer(s) are very important since it influences the whole bearing behavior. The wedge created by the top foil and the rotating shaft generates a hydrodynamic lift inside the lubricant film. The lubricant is required in order to avoid a metal/metal friction when the shaft is rotating. However, the lubricant is still subject to high shear stress which causes power losses by heat generation. Whereas the major drawbacks are :

-Limited load capacity -Thermal management -Design difficulties because of strong fluid-structure interactions Looking at the above-mentioned characteristics, it is clear that GFBs are very well suited for advanced turbomachinery applications such as turbochargers, auxiliary power units, and gas turbine engines [DEL 00b], since low viscosity lubricants allow very high rotation speeds with small friction torque in the fluid film.

Compared to oil bearings, GFBs are also able to work within a larger temperature panel (from cryogenic environment up to 700˚C [DEL 97]) because gas viscosity variations are smaller than oil viscosity ones, and because oil cannot work at these very high temperatures. Besides, if the lubricating gas and the one in the environment are the same -as it happens most of the time in the air -a separate lubricant feeding system becomes unnecessary. Hence, GFBs solve both lubricant supply and bearing environment pollution problems at the same time.

In some cases, GFBs can be used in an environment different from air, like refrigerants in the ECS. Then, the first possibility would be to work with air lubricated GFBs, but inside a refrigerant environment. This way GFBs would behave as usual, but it would also require an extra air lubrication system to avoid contamination. However, contamination cannot be allowed and an extra lubrication system would be really detrimental regarding the efficiency and the compactness of the turbomachinery.

Therefore, a more realistic solution in refrigerant environment is refrigerant lubricated GFBs, exactly as it is done with air lubricated GFBs in air environment. A similar approach has already been studied by Howard [HOW 06, HOW 08a] for power generation in space where the processed fluid is a high-pressure liquefied inert gas such as argon, krypton, or a helium-xenon mixture. In this context, GFBs are a potential solution since if they could use the inert gas as a lubricant.

Over the last 20 years, a significant number of studies have shown GFBs were the best options for a consequent range of applications, such as oil-free turbomachinery and particularly small gas turbines [

Limiting assumptions

The development of multidisciplinary and complex theoretical and numerical models more often than not implies some simplifications when looking for a solution, not only because a model is always based on assumptions regarding the physical world, but also because it helps the understanding of the various effects on a physical behavior. In the present work, some strong assumptions are made for both reasons.

Bearing geometry

First, we choose to focus on the lubricant behavior itself, and not the whole bearing. The present work is strongly related to GFB since we are studying turbo-machinery running with these GFB in very specific contexts. However, we focus on the lubricant behavior only, in the operating conditions of foil bearings. The Thermo-HydroDynamic (THD) approach describes lubricant characteristics such as pressure, density, viscosity, and temperature. This is a possibly restrictive assumption, but we believe that a thorough understanding of the fluid behavior by itself is a necessary condition to understanding the overall GFB performance. The development of a model for the structure deformation of a GFB has already been carried out at Liebherr-Aerospace FRANCE (LTS) and a coupling is feasible (even if not straightforward) between our model and the structure deformation model. It will be presented in future publication(s).

Thermal boundary conditions

A similar approach is adopted for the THD study. LTS have already developed a nodal model for thermal transfers in a GFB in which the real geometry of the different components is taken into account (top foil, bumps, housing, hollow shaft). Therefore, the fluid model is developed under more simple assumptions regarding boundary conditions between the shaft and the fluid and top foil and the fluid. Again, the coupling with the structure thermal transfers in the GFB will be presented in further study.

It is essential to underline here that these two assumptions will strongly modify the results compared to the reality of a GFB. However, the first goal of this work is not to accurately predict the behavior of the GFB as a whole, but mostly to improve the comprehension and the prediction of the lubricant itself. This will also bring a better understanding of the behavior of the entire GFB but we are aware that only the coupling with the structure deformation and thermal transfer models will lead to accurate results regarding global parameters such as load capacity and power losses.

Thermal transfers and phase transition

In this study, we consider phase transition issues in the lubricant film. The phenomena of phase transition is linked to the variation of the enthalpy of the fluid (which is the total energy of the fluid). However, we choose to neglect the enthalpy variation as a first approximation since only a small fraction of the refrigerating gas is liquefied in the process.

The second reason why we make this hypothesis is the following : as we have already explained above, this fluid model will be implemented into a global thermal transfer model for the whole bearing. Only then will we have a good estimation of the energy transfers between the bearing and the lubricant film. Before this step, using the enthalpy for energy transfers in the fluid would be a better approximation but probably not accurate enough to describe the reality of energy transfers inside the whole bearing since the thermal boundary conditions (and therefore the energy transfers) between the bearing and the lubricant film would not be accurate.

Inertia effects

Finally, inertia effects are neglected in the derivation of the Generalized Reynolds Equation since the modified Reynolds number, which characterize the influence of these inertia effects, reaches the value of 1 locally in the worst running conditions, but is more often close to 0.1 or smaller. With increasing working speeds in turbo-machines, we are aware that we are close to the limit where the inertia effects could not be neglected. That is why we will try to integrate this aspect of the theory in further studies.

Chapitre 1 General considerations about gas foil bearings

In this chapter, we sum up the major works available in literature referring to gas foil bearings. 1. General considerations about gas foil bearings

GFBs characteristics

The first 0 • gas bearings were plain, rigid bearings but the main principle is the same as for GFBs : The wedge created by the static housing and the eccentric rotating shaft generates hydrodynamic lift when the two surfaces are moving relatively to each other. Plain journal bearings have good static characteristics such as a high load capacity but they suffer stability problems.

Tilting pad bearings cope with instability problems but they are limited by the necessity to work with small clearance ratios that create geometrical issues in relation with misalignment and thermal expansion. For all these reasons, deformable bore bearings can be more suitable when higher clearance ratio is needed. [GRA 04a] Several flexible bore bearing types already exist but we will only focus on GFBs in this study, whose capacity to work under extreme conditions has already been underlined in the introduction (particularly concerning very low or high temperatures). One of the main benefits of flexible bearings compared to rigid bearings is that they work at a greater radial clearance for a given pressure field in the lubricant, thus generating less shear stress and so reducing power loss. However, they are limited to relatively small load capacity since the peak pressure in the lubricant film is smaller. [HES 83].

GFBs are able to cope with the aforementioned problems. First, friction loss due to the shaft rotation is very small since it is lubricated with gas instead of liquid. In contrast to rigid bearings, foils are pushed away from the rotor with increasing hydrodynamic pressure. This design allows a good adaptation since the clearance ratio remains about the same value to variable running conditions. Moreover, the stability missing in rigid bearings is partly provided by GFBs thanks to Coulomb damping due to foil relative sliding [AGR 97].

The early development of hydrodynamic lubrication

Journal bearings have been used since the antiquity in order to have smooth machine functioning. However, it is only since the XIX th century that scientists have been studying lubrication theory, which plays a key role in the bearing domain. For instance in 1854, Hirn finds that laws describing lubricated contacts are different from those describing dry contacts (Coulomb and Amontons laws) and that torque directly depends on fluid viscosity. He also notices that for relatively small load, water or air are really efficient lubricants. [HIR 54] Concerning bearing functioning comprehension, Petrov sheds the light on the existence of a fluid film that keeps the shaft and the housing separated from each other [PET 83]. He shows that viscosity is the only oil characteristic which plays a significant role in bearing friction. To describe the fluid film, he proposes a famous law to calculate the friction force as a function of speed, contact surface and lubricant viscosity. Moreover, he carries out an interesting global thermal analysis which links the friction force to the temperature.

At the same time, Beauchamp Tower is the first to experimentally study the relation between hydrodynamic pressure generation and journal bearing load capacity [TOW 83].

Effect of the temperature on GFBs

The effect of temperature on the lubricant film will be further discussed. What we want to underline in this section is the influence of the temperature on the bearing itself (housing, foils, etc.) and on its global characteristics. Temperature increases lubricant viscosity in case of compressible fluids, but it might deteriorate the ability of the foil to support the load (temperature gradient might create deformations) [HOW 99].

These phenomena have a combined effect on both global stiffness and load capacity. Below a given temperature, the bearing global stiffness slowly decreases with increasing temperature whereas the trend becomes significant above this temperature. It appears that load capacity follows the same trend. [GRA 04a] Moreover, GFB foils surface treatments are rather complex and use very thin layers of specific material. It is very likely that these material properties deteriorate with high temperatures and/or strong temperature gradients. Therefore, the study of GFB thermal behavior is crucial if one wants to give accurate predictions when GFBs are working under extreme thermal conditions such as cryogenic environment, close to heat sources, very high rotational speeds, etc.

Predictive GFB models

In 1953, Blok and van Rossum [BLO 53] introduce the concept of GFBs. They point out that the GFB film thickness, larger than the rigid gas bearing one, can improve operational reliability and provide a solution for problems related to thermal expansion of the shaft and the housing. Very valuable experiments on GFBs are carried out by Licht [LIC 70], whose conclusion is that GFBs are recommended for high-speed turbomachinery.

Field experience has proved, since the late 1960s, that GFBs are far more reliable than ball bearings previously used in ACMs. Since then, GFBs have been used in almost 1. General considerations about gas foil bearings every new ACM, in both civil and military aircraft [AGR 97]. Agrawal also sums up the aforementioned assets of GFBs such as high rotational speed and extreme temperature capabilities. However, the use of GFB for turbomachinery applications requires accuracy in modeling capabilities.

Predictive GFB Models of Heshmat et al. [HES 83] first presents analyses of bump type GFBs and details the bearing static load performance. The predictive model couples the gas film hydrodynamic pressure generation to a local deflection of the support bumps. In this simple model, the top foil is altogether neglected and the elastic displacement is proportional to the local pressure difference through a structural compliance coefficient which depends on the bump material, thickness and geometric configuration. Peng and Carpino [PEN 93,PEN 94] presents finite difference formulations to calculate the linearized stiffness and damping force coefficients of GFBs. The model gives the bump stiffness which is used in order to make calculation on bearing global parameters. In the model of the underlying foil structure, a perfectly extensible foil is placed on top of the corrugated bumps. The model includes the equivalent viscous damping of dry-friction between the bumps and the bearing housing. As the dry-friction coefficient increases, the direct damping coefficients significantly increase.

Yu et al. [YU 05

] publish an article in which they detail a generalized solution of an elasto-aerodynamic lubrication problem, where they study both static and dynamic deformations of foils and gas film behavior simultaneously. But recent studies have shown that structural damping as well as dry friction are responsible for non-linear dynamic behavior. [BAR 11] Chapitre 2

Hydrodynamic study and two-phase flow

For an accurate description of GFB lubricant behavior, we need to consider fluid compressibility. Most of the studies use the ideal gas law because it is an accurate and easy way to deal with fluid compressibility in most running conditions. The main assumption of the ideal gas law is that density is proportional to pressure. However, when the pressure is close enough to the vapor pressure, this assumption becomes irrelevant because of the vapor/liquid transition. Furthermore, density is also strongly influenced by temperature. Internal experiments at Liebherr-Aerospace FRANCE on new refrigerant-lubricated compressor designs have shown that under specific operating conditions, a mixture of vapor and liquid appears in the compressor, instead of a single-phase vapor flow. One of the main goal of this study is to accurately describe the refrigerating gas acting as a lubricant under single and two-phase flow. The theoretical work of Reynolds [REY 86] are generally considered to be the hydrodynamic lubrication theory basis since they lead to the famous "Reynolds equation". The Reynolds equation is derived from the mass and momentum conserving equations under thin film conditions and it governs the hydrodynamic pressure generation and the lubricant flow.

After Reynolds' major breakthrough in the hydrodynamic lubrication theory, scientific research carried on focusing on the development of analytical models since computers did not exist yet. Solving the Reynolds equation and its boundary conditions could not be performed analytically. Therefore, simplifying assumptions had to be done.

In 1904, Sommerfeld [SOM 04] proposed a long journal bearing approximation theory that consists in deriving an analytical solution to the Reynolds equation that neglects the effect of the bearing edges (infinitely long bearing approximation) as well as the occurrence of the film rupture.b Neglecting this phenomenon causes the long bearing approximation to predict unrealistic load carrying capacities. Film break-up modeling then becomes the subject of numerous works. We will not study this phenomena which happens only with liquid lubricants.

In 1914, Gumbel [GUM 14] improves Sommerfeld's solution by presenting the half-Sommerfeld theory. In this approach, the negative pressures at the divergent portion of the gap obtained with the former theory are turned to zero (reference pressure) instead. By setting pressure to zero, Gumbel tries to describe the previously observed film break-up phenomenon. Swift [SWI 32] further improves the film rupture model by imposing a pressure gradient continuity at the rupture boundary. These conditions are those that give the most accurate solutions and they have been subject to numerous experimental studies.

The short bearing approximation

The problem with the long bearing approximation originates from the fact that it is only valid for length to diameter (L/D) ratios above 4, whereas most of the bearings used in the industry today are under this limit, because of misalignment problems. However, an other approximation, called "short bearing approximation", is often used when the L/D ratio is below 1/4.

Hydrodynamic study and two-phase flow

This approximation is based on the fact that small L/D ratio means negligible circumferential pressure gradients (in comparison with axial pressure gradients). With such an approximation, a simplified analytical solution can be found for academical work. The short bearing approximation unmistakably helps solving the problem more easily but it reduces the accuracy of the model as well.

The need for numerical methods

As we have seen before, the hydrodynamic lubrication problem must be approximate in order to obtain a purely analytical solution. In some cases, these approximations are justified and give satisfactory results. However, complex phenomena such as lubricant film rupture or thermal effects cannot be accurately described under these assumptions : The model then becomes over simplified. With the development of computers, scientists start using computer-based numerical methods in order to described the aforementioned neglected phenomena.

The first computer-based models are developed by Raimondi and Boyd [RAI 58] and Pinkus [PIN 58]. After these pioneer works, authors try to describe reality as accurately as possible, keeping in mind that fewer simplifications generally means higher calculation costs. Progress in computer-based calculations finally allow authors to add thermal and elastic behavior to the hydrodynamic problem itself.

The Reynolds model

The Reynolds model is based on the Reynolds equation under the Swift-Stieber conditions (both pressure and pressure gradient are set to zero at the rupture boundary), solved with an algorithm developed by Christopherson [CHR 41]. This algorithm uses a finite difference discretization and the Gauss-Seidel iterative process. At the nodes localized inside the film-break up region, the pressure is set to the reference pressure in each iteration until convergence is achieved and reformation is set either at the maximum film thickness angular position or at the grooves. This model has been very popular over the past decades and is still used in some works today as it is a simple numerical model able to give satisfactory results concerning the pressure profile in plain journal bearings [KHO 86,BON 86].

It is also used for tilting-pad bearings [KNI 90, FIL 97] and recently, several authors used the Reynolds conditions in their models [BRI 03, MAJ 04]. Unfortunately, this model generally does not give good results concerning the flow rate or the mass conservation in a two-phase lubricant, the critical point being the liquid-gas interface.

A Generalized Reynolds Equation (GRE)

When Dowson, in 1962, publishes a paper introducing a generalized form of Reynolds equation, he explains that his work comes from a lack of a general theory in thin-film lubrication when fluid parameters are subject to strong cross-film variation [DOW 62]. For this reason, Dowson comes up with its GRE which accounts for transport property variations, density and dynamic viscosity, as a function of pressure and temperature.

Thermodynamic equation of state (EoS)

An EoS is a general relation between three local fluid parameters : local pressure, temperature and density. The ideal gas law which is among the most current EoS gives a linear relation between pressure and density and shows really good agreement with experimental data when the conditions are far enough from the liquid/vapor transition.

The Van der Waals description is a semi-empirical relation which has the advantage to give more accurate results with a relatively simple relation : The Van der Waals EoS describes both vapor/liquid transition and supercritical state. The Van der Waals EoS is based on two main physical phenomena : the intrinsic volume of the molecule of gas B w and the interaction between the molecules of gas, which is called the "cohesive pressure" of the gas A w [BON 64] :

p = RT M ρ -B w - A w M ρ 2 (2.1)
More generally, the so called "cubic" EoS category groups together all the EoS which consists in a 3 rd order polynomial form : Density is the unknown and the pressure is one of the parameters. The Van der Walls EoS belongs to this category and was in fact the first of its kind. Several EoS were found as different ways to describe fluid behavior but some of them are more suitable to describe refrigerating gas in our model. That is why we chose the modified Peng-Robinson equation [PEN 76] which gave the best results in term of accuracy and suitability to the numerical model in our comparative study.

The Viriel development is an other category of EoS in which the ratio of the volume over the temperature is a function of a polynomial expansion in power of 1/V . The Van der Waals equation also belongs to this category as a second-order development.

Thermodynamic model for two-phase flow -A nonlinear equation of state 2.2.1 Symbol description

From now on, a given parameter with the subscripts • x , • y , • z will correspond to the given parameter for the respective direction x, y, z, the subscript • 0 will correspond to the reference value in the surrounding environment, the subscripts • L , • V will correspond respectively to the liquid and vapor phase and the subscripts • S , • T will correspond respectively to the value at the shaft and the top foil. Also, the exponent • t will correspond to the turbulent regime , the exponent • * will correspond to the sum of the laminar and turbulent parameter values, or "equivalent" parameter.

Choosing the equation of state

In bearings, for compressible gas, the pressure (given by solving the Reynolds equation) is related to density, temperature and viscosity. This relation is often described by the ideal gas law. The main assumption of the ideal gas theory is that density remains dependent on pressure and temperature.

However, when pressure is close enough to the vapor pressure, which is value for which the vapor/liquid transition occurs, this assumption becomes invalid. Two-phase flow has been observed experimentally under specific conditions when refrigerant is introduced into GFBs' lubrication system. When working with air-lubricated bearings, the ideal gas law is used in its valid domain. The ideal gas law gives a linear relation between pressure and density.

Therefore we use a non-linear EoS able to describe the density variation as a function of pressure and temperature, as well as the vapor/liquid transition. If a lot of EoS are able to describe non-linear behavior, only a few match our requirements. First of all, the EoS has to describe the behavior with enough accuracy so that it does not deteriorate the precision of the whole model. Besides, this EoS is only one part of a far more complex model (THD) and should be simple enough not to increase the complexity of the model. In fact, density and pressure are strongly coupled and both the GRE and the EoS have to be solved simultaneously. Viscosity is also linked to the temperature and to the fraction of liquid in the fluid. As a result, it also depends on this modifications Finally, this THD model is bound to be adapted to different lubricant types. As a result, parameters used in this model should be easily modifiable.

A modified Peng-Robinson equation of state

Considering all these criteria, the most convenient type of equation for this model is the "cubic" EoS which consists in a 3 rd order polynomial form : Density is the unknown and the pressure is one of the parameters. We choose a modified Peng-Robinson EoS [PEN 76], whose dimensional formulation is written as follow :

p = RT M ρ -B - a ρ M ρ 2 + 2MB ρ -B 2 (2.2)
With : R = 8.314462175 (2.3)

B = 0.0777960739039 • RT c P c (2.4) m ρ = 0.378893 + 1.4897153ω ρ -0.17131848ω 2 ρ + 0.0196554ω 3 ρ (2.5) a ρ = 0.457135528921 • R 2 T 2 c P c 1 + m ρ 1 - T T c 1/2 (2.6)
In order to solve the EoS, we put it in a 3 rd order polynomial form written as follow :

pB 3 + RT B 2 -a ρ B ρ M 3 + -3 • pB 2 -2 • BRT + a ρ ρ M 2 + (pB -RT ) ρ M + p = 0 (2.7)
Where p is the pressure, R the ideal gas constant, T the local temperature, B and m ρ are EoS constant coefficients (for a given fluid) whereas a ρ is an EoS temperaturedependent coefficient. The molar mass M, the critical pressure P c , the critical temperature T c and the a-centric factor ω ρ are among the fluid characteristics which can be found in Refprop or any other fluid database.

Hydrodynamic study and two-phase flow

Theoretical model validation

This type of EoS gives a more accurate description of a real gas physical behavior compared to the ideal gas law. It takes into account the volume of gas molecules and their local interactions (see section 2.1.6). There are a lot of EoS using this approach that would have been able to describe the refrigerant behavior. The aim of this study is not to create a specific EoS for the numerical model and a comparison between what seemed to be the most relevant laws allowed us to chose the Modified Peng-Robinson EoS to implement our numerical model. Figure 2.1 shows the comparison between the EoS in our model with our thermodynamic reference "REFPROP" (NIST). The symbol REFPROP denotes a reference database of thermodynamic properties which will be our reference of comparison to validate our results. We plotted the density versus pressure for several temperature. On this figure, we can see that the pressure and temperature-dependent fluid density can be predicted using this EoS. For tempera+ture which does not exceed the critical temperature,we see that density is predicted with less than 2% of error in vapor phase (low pressure) and less than 5% of error in liquid phase (high pressure) as long as the fluid is not close to the vapor/liquid transition (discontinuity for density). Above the critical temperature, the discontinuity does not exist and the fluid density is predicted with less than 2% error.

The only lack of accuracy of the EoS prediction is due to bad prediction of the vapor pressure, for which the transition occurs. In order to cope with this, we use the Clapeyron formula to predict the vapor/liquid transition with an error smaller than 1%, lowering our maximum error down to 2% for vapor phase density prediction and 5% for liquid phase.

The Clapeyron formula

The Clapeyron formula describes the thermodynamic behavior of a gas, and more precisely its transition from a vapor to a liquid phase, giving a relation between the temperature and the value of the vapor pressure. Thanks to the Dupré approximation, one only needs three couples (T i , p sat (T i )) (for a given gas) to predict the vapor pressure at any given temperature [LOB 01] :

ln(p sat (T )) = C 1 - C 2 T -C 3 ln(T ) (2.8)
C 1 ,C 2 ,C 3 are three phenomenological constants (for a given gas) which values can be obtained thanks to three couples (T i , p sat (T i )). These three couples can easily be found in the literature for a large panel of fluid, such as refrigerants, or using REFPROP. In conclusion, this EoS can be solved with a few characteristic numbers given as an input by the user : molar mass, critical temperature,critical pressure and three couples (T i , p sat (T i )).

The vapor/liquid transition issue

The non-linear model can describe a two-phase flow in the lubricant. It means that for a given set of running conditions, the pressure in low-pressure zones might be under the vapor pressure whereas in some high-pressure zones it might exceeds this limit. In this case, some of the lubricant might change from vapor to liquid phase. This is one of the major numerical difficulty in this study : The strong variations in density is difficult to simulate since it generates strong perturbations in the iterative numerical process.

In order to take compressibility into account, a lubricant model uses an EoS to describe the density variations with pressure and temperature. Most of the time, the ideal gas law is enough to describe the compressibility. But a fluid particle infinitely close to the shaft has the same velocity as the shaft surface (under no-slip assumption). Therefore, when the rotational speed is high enough, this particle will go through the high-pressure zone thousands of times per second. In the mixture zone, There is formation of small liquid particle in suspension in the gas. We apply a transitional law to link local pressure and density and compute the mixture equivalent local properties (density, molecular viscosity, thermal conductivity).

As we will see in the next section, we choose to use a vapor/liquid transition model which is not directly linked to the enthalpy calculation in order to compute the local fraction of liquid in the lubricant in two-phase flow scenario. This is a first approximation that we choose to make because it simplifies the THD model and the numerical process in severe running conditions. The model we choose has shown good efficiency in previous works [ ODY 03]. However, the full thermodynamic calculation (enthalpy) is an aspect of the problem we would like to study in further works in order to refine the model and also in order to adapt it to transient calculations.

A vapor/liquid transition model

We used a vapor/liquid mixture model [ODY 03] to describe two-phase flow. A sinusoidal transition, function of the local pressure, the vapor pressure and the minimal speed of sound a min in the mixture accurately describes mixture equivalent parameters. The subscript . V is for the value of a given parameter in the vapor phase, . L in the liquid phase. Within the transition zone, the value of the given parameter is between the theoretical values X V and X L and depends on the proportion of liquid in the mixture. This transition occurs for three parameters : density ρ, thermal conductivity k and molecular viscosity µ. Thanks to this model we can calculate a range of pressure from the vapor pressure p sat to its upper limit p sat + ∆p for which the fluid equivalent parameters are calculated as follows :

a V = γp ρ (2.9) a min ≈ 2 • a V ρ V ρ L (2.10) ∆ρ = ρ L -ρ V (2.11) ∆µ = µ L -µ V (2.12) ∆k = k L -k V (2.13) ∆p = π 2 a 2 min ∆ ρ (2.14) ρ = ρ V + ∆ρ sin p -p sat ∆ρa 2 min (2.15) µ = µ V + ∆µ sin p -p sat ∆ρa 2 min (2.16) k = k V + ∆k sin p -p sat ∆ρa 2 min (2.17)
Behavior of refrigerant in laminar isothermal biphasic running conditions a V is the minimal speed of sound in the vapor phase [BOH 88] calculated from pressure, density and the adiabatic index γ (Assumed to be 7/5 as it is often for diatomic molecules from kinetic theory) for an ideal isentropic gas a first approximation. Then we compute the minimal speed of sound in the mixture from the vapor over liquid density ratio. This approximation is only possible when ρ L ρ V , a condition which is always satisfied within our running conditions range, where liquid density is at least two orders of magnitude higher than the vapor density. In this model, viscosity and thermal conductivity are described with the same sinusoidal law as density in the mixture.

Behavior of refrigerant in laminar isothermal biphasic running conditions 2.3.1 Dimensionless parameters

In this section, we sum up the different characteristic dimensions and the related dimensionless parameters used in the Reynolds equation :

Parameter name

Symbol Dimensionless form Circumferential coordinate

x x = x • 1/R a Radial coordinate y ỹ = y • 1/C Axial coordinate z z = z • 1/L Thickness h h = h • 1/C Pressure p p = p • 1/p 0 Density ρ ρ = ρ • 1/ρ 0 Dynamic molecular viscosity µ μ = µ • 1/µ 0 Circumferential velocity component u ũ = u • 1/(R a Ω) Radial velocity component v ṽ = v • L/(R a ΩC) Axial velocity component w w = w • 1/(R a Ω)

The Reynolds equation

In this first approach, the pressure field in GFBs is obtained solving the steady-state isothermal Reynolds equation for laminar, compressible fluids [REY 86]. It can be written in its dimensionless form as follows :

∂ ∂ x ρh 3 12μ ∂ p ∂ x + R a L 2 ∂ ∂z ρh 3 12μ ∂ p ∂z = Λ ∂ ρh ∂ x (2.18) FIGURE 2

.2: Simplification in the model

With the bearing number Λ being :

Λ = µ 0 R 2 a Ω p 0 C 2 (2.19)
And the circumferential and axial velocity components :

ũ = 1 Λ 1 2μ ∂ p ∂ x ỹ( ỹ -h) + h - ỹ h (2.20) w = R a L 1 Λ 1 2μ ∂ p ∂z ỹ( ỹ -h) (2.21)
Besides, even if the pressure is constant along the radial direction, the velocity radial component is not zero. In fact, one can extract it from the mass conservation equation written in its dimensionless form as follows [BON 86] :

∂ ρ ũ ∂ x + ∂ ρ ṽ ∂ ỹ + R a L ∂ ρ w ∂z = 0 (2.22)
2. 4 From the isothermal laminar model to the THD turbulent model

A Generalized Reynolds Equation

Viscosity and density terms are replaced by integral coefficients that takes into account the cross-film variations. Coefficients are separated in two groups : The first group only contains coefficients without density-derivative terms in the radial direction and the other group contains coefficients with density-derivative terms in the radial direction.

Most of the time, density can be considered constant across the film and, as a result, the terms in the second group are set to zero in the GRE. More generally, if all radial variations are negligible in the GRE, one finds the classic Reynolds equation formulation.

The GRE [DOW 62] accounts for turbulent flow in the equivalent-viscosity terms. The equivalent viscosity combines the molecular viscosity, an intrinsic fluid property, and the "eddy" viscosity, which is due to turbulent flow whenever it occurs. If the flow is laminar, then the "eddy" viscosity is set to zero and the equivalent viscosity will be equal to the molecular viscosity. (For more details, refer to the chapter dedicated to turbulent flow) The GRE can be written in its dimensionless form as follows :

∂ ∂ x ( F2 + G1 ) ∂ p ∂ x + R a L 2 ∂ ∂z ( F2 + G1 ) ∂ p ∂z = Λ h∂ ρS ∂ x - ∂ ∂ x F3 + G2 F0 + ( hV ) S (2.23)
With the dimensionless viscosity integral coefficients being :

2. Hydrodynamic study and two-phase flow

F0 = h 0 d ỹ μ * (2.24) F2 = h 0 ρ ỹ( ỹ -ỹ) μ * d ỹ (2.25) F3 = h 0 ρ ỹ μ * d ỹ (2.26) G1 = h 0 ỹ∂ ρ ∂ ỹ ỹ 0 ỹ μ * d ỹ -ỹ ỹ 0 d ỹ μ * d ỹ (2.27) G2 = h 0 ỹ∂ ρ ∂ ỹ ỹ 0 d ỹ μ * d ỹ (2.28)
The cross-film viscosity variations are accounted for integral coefficients whereas two-phase flow transition theory ensure local parameter continuous variations. The two-phase flow phenomena is treated with a single GRE which is identical with either single or two-phase flow occurs.

Dimensionless parameters

In this section, we sum up the different characteristic dimensions and the related dimensionless parameters used in the GRE :

Parameter name

Symbol Dimensionless form Viscosity integral coefficients

F 0 F0 = F 0 • µ 0 /C Viscosity integral coefficients F 2 F2 = F 2 • µ 0 /(ρ 0 C 3 ) Viscosity integral coefficients F 3 F3 = F 3 • µ 0 /(ρ 0 C 2 ) Viscosity integral coefficients G 1 G1 = G 1 • µ 0 /(ρ 0 C 3 ) Viscosity integral coefficients G 2 G2 = G 2 • µ 0 /(ρ 0 C 2 ) TABLE 2
.2: Dimensionless fluid parameters for the GRE These coefficients come from the integration of fluid properties across the film thickness. Compared to the classic Reynolds equation, the analogy can be made thanks to the dimensions of the parameters as shown in (Table 2.2)

Model improvements

To understand the advantages of the GRE compared to the Reynolds equation, one has to go back to the demonstration of these equations. If the inertia forces are neglected and the thin-film hypothesis is applied to the Navier equations, they can be written in their dimensionless form as follows :

∂ p ∂ x = ∂ ∂ ỹ μ * ∂ ũ ∂ ỹ (2.29) ∂ p ∂ ỹ ≈ 0 (2.30) ∂ p ∂z = ∂ ∂ ỹ μ * ∂ w ∂ ỹ (2.31)
From that point, we integrate this set of equation in order to find an analytical formulation of the velocity field. Then we use this formulation in the continuity equation and integrated it across the film thickness. If we assume that local parameters are independent from the radial coordinate, then it leads to the Reynolds equation. If not, it leads to the GRE with the corresponding fluid velocity components written in their dimensionless form as follows :

ũ = 1 Λ ∂ p ∂ x ỹ 0 ỹ μ * d ỹ + 1 F0 - 1 Λ ỹ ∂ p ∂ x ỹ 0 d ỹ μ * (2.32) w = R a L 1 Λ ∂ p ∂z ỹ 0 ỹ μ * d ỹ - ỹ∂ p ∂z ỹ 0 d ỹ μ * (2.33)
The radial velocity component can be extract from the mass conservation equation (2.22) as before.

The 3D-THD study makes the variable viscosity along the radial direction mandatory. We would like to underline the strong theoretical link between the GRE and the energy transfer in the film. The combined effects of turbulent flow and temperature gradient across the film thickness is taken into account into the integral coefficients in the GRE and modifies pressure and density fields. They have a direct impact on global parameters, such as load capacity and power loss.

Thanks to the theory, the strong link between the GRE, the 3D thin-film energy equation, a temperature-dependent EoS and an appropriate 3D turbulent model clearly appears.

Conclusion on the hydrodynamic study and twophase flow

In this chapter, we set the basis of a two-phase flow, 3D, turbulent THD model. We implemented an adequate EoS able to describe two-phase flow parameters, with a smooth transition from vapor to liquid phase (mixture). We also implemented a GRE, able to account for two-phase flow as well as local parameters cross-film variations, due to THD and turbulent phenomena.

The basis of the hydrodynamic lubrication model are now integrated to the model. However, in order to obtain realistic values of the pressure and density fields, we need to have an efficient model of other local parameters. Temperature variations directly influences density and viscosity in the fluid, as well as bearing structure material properties. Turbulence affects viscosity and thermal conductivity. Therefore, the next step is to accurately predict the influence of these phenomena.

Chapitre 3 ThermoHydroDynamic model for GFBs

Lubricant local temperature has an influence on bearing parameters (density, viscosity, vapor pressure) and heat transfers impact on the bearing structure itself. A THD model enables better predictions of both local and global parameters. Moreover, the two-phase flow strongly depends on local temperature variations. Adequate thermal management is necessary when designing high temperature applications using GFBs, such as gas turbine engines [DYK 04]. Hori and Kato [HOR 08] explain that heat generation in the lubricant film and the related temperature rise are the most important factors in journal bearings, particularly for high rotational speeds because power loss is almost only converted into heat.

A rise in the average bearing temperature as well as the local peak pressure can be very detrimental regarding bearing properties such as load capacity, minimum film thickness, etc. [MA 96] Not only does the rise in temperature modifies some of the lubricant properties, but it also generates deformations of the shaft, the foils and the housing. Furthermore, in the particular case of two-phase flow, the vapor pressure strongly depends on the lubricant local temperature.

In the early developments of bearings theoretical studies, thermal effects are first neglected or at least too much simplified in journal bearing analyses because of the complexity of the phenomena involved. Then authors start studying heat effects in lubricated film [ KIN 33,HAG 44] and show more detailed thermal models are mandatory in order to avoid the discrepancies observed between experimental data and the first theoretical predictions.

Lubricant film thermal study

Global thermal model

This is the simplest approach of thermal phenomena, based on the assumption that the lubricant temperature and viscosity are constant in the lubricant (constant in space, and in time since we are uner the steady-state assumption). Thanks to a global energy balance, one calculates the effective average temperature and viscosity [SWI 37]. In this particular case, the thermal aspect is even more important for high-speed GFBs since the amount of heat produced by viscous dissipation cannot be neglected [BOU 04].

The isoviscous assumption is legitimate as a first approach or when the model includes other theoretical aspects which, combined to a THD approach, would require excessive computational costs [TAN 00, STA 01, TAL 02, YU 02, JOH 04, BOE 04, IGH 08]. However, this study involves complex lubricant behavior related to local temperature, which has to be calculated accurately.

Evolution to a local thermal model

In order to cope with the lake of efficiency of global thermal methods, authors add complexity to there models. The first step is to consider that temperature and viscosity are not constant along the circumferential direction [COP 49] which is the flow main direction. The next step is to consider that, because of the temperature gradient between the top foil and the shaft, temperature and viscosity are functions of the radial coordinates as well [ZIE 57].

Viscosity variations along the bearing circumference do not modify the Reynolds equation itself whereas considering viscosity variation across the film thickness has a big theoretical impact. It implies a modified Reynolds equation called "Generalized Reynolds Equation" (or GRE), which integral terms for the viscosity result from the integration of the continuity equation across the film thickness.

The first Thermohydrodynamic (THD) analysis

With the publication of "A generalized Reynolds equation for fluid-film lubrication" [DOW 62] comes the first THD model, which is still fundamental for many authors. In this model, the GRE and the energy equation are solved in a coupled iterative mode, because they are strongly non-linear and coupled through the lubricant local viscosity.

During the 1970's, several authors point out the need for model including radial temperature variations [MCC 70, SEI 72]. In practice, most of the following works consist either in improving boundary conditions for heat transfers between the fluid, the shaft and the housing or taking into account more phenomena in the model (like film break-up, turbulence, inertia effects, etc.).

After that, authors implement more complex THD analysis where lubricant properties are studies in severe conditions and/or unusual environments, such as cryogenic gas.

[BRA 87]

Possible simplifications inside the lubricant film

Constant fluid film parameters Setting fluid parameters as constants, such as viscosity, density or specific heat allows major simplifications in bearing analysis. However, when working at high rotational speed, it has been shown [CHU 04] that variable density, viscosity and specific heat were mandatory in order to predict bearing load and frictional power loss accurately. Sharma et al. [SHA 03] find in their TEHD study of a hybrid journal bearing system that viscosity variations due to the temperature rise have significant influence on both static and dynamic performance on the system as well.

The axial temperature gradient When referring to the literature, the axial thermal gradient is generally small compared to other significant values (radial and circumferencial temperature gradients) and most of the time a 2-D thermal analysis gives realistic results for standard journal bearings. 2-D thermal analysis is also used when dealing with rather complex problem such as thermal transient behavior [FIL 97, KUC 00, TOU 03], or when the purpose of the study is to give a rapid estimation of bearing parameters [JAN 04, WAN 11b].

However, recent studies on GFBs [SAN 09] using 3-D thermal models, in which predictions are in good agreement with test data, show that the axial thermal gradient increases with shaft rotational speed and can no longer be neglected at very high speeds : Not only does the use of 2-D thermal approximation loose information concerning bearing edges but it overestimates values as well [MIC 07]. The study of refrigerating gas acting as a lubricant is one of the case were the use of 3-D model helps to accurately localize the possible liquid phase.

The Couette approximation The Couette approximation consists in neglecting all the convective heat transfer terms which are related to the pressure in the energy equation. Considering only the drag-driven flow as being responsible for heat transfer considerably simplify the coupling between the GRE and the energy equation. Even though the simplification might be a good way to reduce computer cost, the results obtained with this method (like the temperature profile for instance) are not realistic enough to be considered satisfactory [SAL 01].

Cross-film variations

The "standard" Reynolds equation is the result of the integration of continuity across the film thickness, accounting for the thin film hypothesis. The GRE for its part additionally takes into account the local viscosity variation across the film thickness through integral coefficients. The THD problem, whose purpose (among others) is to solve both GRE and energy equation in the lubricant film iteratively (where the radial coordinates appears explicitly), cannot ignore the cross-film dimension.

After the first 1-D thermal models have revealed their lake of accuracy, authors tried different methods to create a model that will describe cross-film phenomena accurately without being time-consuming. It would be too long to detail every single technique employed in order to approximate cross-film value but the trend clearly shows that not considering the cross-film dimension in a THD analysis would dramatically deteriorate results accuracy [SIN 08].

ThermoHydroDynamic model for GFBs

The common principle of several methods is to approximate local viscosity and temperature profile across the film thickness as accurately as possible but with a minimum of points thanks to polynomial approximations. Vijayaraghavan [VIJ 95] explains that since the coefficients in the GRE are functions of integral of viscosity across the film, it is very important to determine the viscosity profile. He uses the Legendre method, which is still widespread [SHY 04, MOR 07, SHY 06, SHY 08].

Several authors also use a cross-film averaged energy equation thanks to temperature profile polynomial approximations, whose order may vary from one to another [FAT 06, STE 09].

Lubricant mixing in the groove

Since "standard" journal bearing groove dimension is small compared to the circumferential dimension, one might think there is no need to treat the groove area with particular caution. However, previous works have revealed that reality is a little bit more complex.

The first problem that may occur when one tries to solve the THD problem in the groove area is that the thin film assumption used in the GRE is no longer acceptable, since the dimensionless film thickness reaches the same order of magnitude than other dimensionless length. Therefore, the THD problem in this area should be treated differently [BRI 09].

Most of the authors choose to perform a global heat and mass balance at the boundaries of the groove in order to estimate the average temperature of the fluid at the leading edge of the groove (or inlet section temperature). Some of them use uniform temperature profiles but others choose a non uniform temperature profile, determined by energy conservation and boundary conditions, thus preventing discontinuities at the interface.

It has also been shown that for particular types of regime, like turbulent flows in high speed journal bearing, the effect of the mixing at the groove should always be a part of the analysis.

Shaft thermal study

The shaft has a specificity compared to other parts : a given point on a journal bearing inner housing surface always stay at the same position compared to the minimum film thickness localization at steady state, whereas the order of magnitude of the shaft rotating speed suggests that there is a circular symmetry concerning thermal parameters in the shaft [DOW 66, VIJ 95, ?]. Experiments have confirmed this theory : the circumferential variation of the inner housing temperature is significant while the one of the shaft can be ignored.

Even if the shaft temperature is set as constant, one still need to determine a realistic value for it. In some case, a prescribed temperature can be set if there is a particular shaft environment or if the value comes from experiments. Other authors simply set the shaft temperature to a value equal to the average housing surface temperature.

Of course, classic heat transfer calculations are possible to determine the temperature at the film/shaft interface [PIE 00]. Lee and kim [LEE 10] carry out a numerical model in which 3D energy equation is used inside the lubricant and 1D energy equation is used in the shaft (Also used by Banwait [BAN 98]). The analysis shows that the shaft temperature increases with load and rotational speed and decreasing bearing clearance. The shaft temperature profile is parabolic with its maximum being a function of the boundary conditions.

GFBs are intrinsically subject to deformation since it is one of their major assets. However, adequate models of GFBs should consider several aspect regarding elastic deformations inside the bearing due to temperature and shaft rotation [SIM 08]. For example, Michaud et al. [MIC 07] reveal that not considering elastic deformation inside a journal bearing can lead to underestimations regarding the temperature of the bearing components.

Housing thermal study

The housing-film interface

The treatment of the housing-film interface is much more difficult than the shaftfilm one, because of the circumferential coordinate dependence. The most common assumptions for heat transfers between the lubricant and the housing are "isothermal" and "adiabatic", and they are limiting cases. These two assumptions are unrealistic and cannot lead to a good predictions for the whole bearing, because they do not describe heat transfers accurately. However, they are very useful in order to study one part of the bearing. In our case, we will use the adiabatic condition as a first approximation at the interface between the lubricant and the housing. When this work will be coupled with a more realistic structural thermal model for GFB, this assumption will no longer be needed.

As we said before, classic heat transfer calculations are also possible, which means setting heat flux continuity at the housing-film interface and using the Fourier equation inside the housing body. Experimental works confirm that there is a non negligible radial heat flux through the housing body. The radial heat flux means that it is pertinent to consider heat flux continuity at the housing-film interface, at least for a wide range of operating conditions.

Wang et al. [WAN 02

] also proposed to consider the shaft-film-housing system as a whole. What it means is that the interfaces between the lubricant and the solid bodies are treated as internal regions in the model. Therefore, no specific conditions are needed at the interface.

Heat conduction inside the housing

Experimental work has shown that there is a non negligible part of the energy dissipated inside the fluid film which is transfered to the housing. For instance Ma and Taylor [MA 96] show in experimentations on plain journal bearings that the temperature differences between the bearing liner and the back surface were considerable in the vicinity of maximum temperature location for high speed and load operating conditions. That is the reason why many authors use the Laplace equation in one, two or even three dimensions to describe heat conduction inside the housing.

The most common case is a 2-D Laplace equation approach, where temperature is supposed to be independent of the axial direction and where the thermal study is then limited to the bearing mid-plane. The wide use of the 2-D model is justified by experimental works that show the axial temperature gradient is negligible in most of the case and so a 3-D approach is only necessary in specific cases such as non symmetric geometry for misaligned bearings.

In GFBs, the structure is far more complex than in plain journal bearings and heat flux can generate important thermal structural effects [TAL 11]. That is why 3D THD study is more important in GFB model compared to plain journal bearing. As the film is the major heat generation location, it affects the whole structure thermal behavior. Structural dynamic coefficient are also influenced by these effects, either directly by heat flux in the lubricant film, or indirectly by the shaft temperature [SAN 11].

The housing-environment interface

At the outer housing surface, free convection is the most common boundary condition used.

3-D energy equation models

As we have seen in previous paragraphs, simplifications are not always possible when dealing with rather complex phenomena inherent to THD problems. Numerical methods and computer performance improvements allow scientists to implement

Relations between bearing parameters and thermal behavior

When studying bearing behavior, several trends can be found. Theoretical and experimental analysis allow authors to predict temperature variations inside the bearing versus a given global parameter (A parameter that reflects the bearing global behavior, such as load capacity, power loss) [SHY 12]. In this chapter, we tried to sum up the most important trends :

-Thermal effects in a highly loaded plain journal bearing strongly influence minimum film thickness, maximum pressure and power loss [BOU 04].

-In GFBs, the temperature inside the lubricant film increases with the shaft rotational speed and the static load. In both cases, the temperature rise is the result of a higher lubricant velocity gradient across the film, which creates more shear stress. However, it is the shaft speed, more than the static load, which influences the temperature rise in the film [ZHA 00, RAD 04, SAN 09, FEN 09]. Besides, Lee and kim [LEE 10] find that the rotor and the bearing sleeve temperatures can be approximate by parabolic functions of the shaft speed.

-The load capacity and the torque on GFBs significantly increase if temperature variations are taken into account [PEN 06,FEN 09].

-The maximum film temperature in GFBs is located just downstream the maximum gas pressure along the circumferential direction [FEN 09, SAN 09], and closer to the top foil than the shaft surface along the radial direction [FEN 09].

-A direct relation exists between the radial clearance and the thermal behavior of a GFB. Radil et al. [RAD 02] show that, for a given GFB, an optimum radial clearance exists, for which the GFB produce its maximum load capacity coeffi-cients. This last relation is crucial since smaller radial clearance (compared to the optimum) can cause thermal runaway.

Thermal management

Thermal management is the capacity to control the bearing thermal behavior in the field. It helps to prevent GFBs from being subject to thermal seizure as well as maintaining a desired load capacity and stability. In 2006, Dykas [DYK 06] finds that inadequate GFB thermal management is one of the two basic mechanisms that limit GFB performance. Nowadays, thermal management still remains a key issue regarding GFBs and is considered by Lee [LEE 10] to be "as important as rotordynamic stability and load capacity" .

There are several way to improve heat evacuation and prevent bearing from failure. Several authors have been studying cooling systems that carry away the heat generated inside the bearing, for example with cooling streams inside a hollow shaft (inner flow stream) or between the bearing housing and the top foil (outer flow stream) [SAN 10]. Radil et al. [RAD 06,RAD 11] show that the outer flow stream cooling method is the most effective method compared to direct or non direct inner flow stream cooling method since it can efficiently lower both bearing internal temperature (which can reduce materials properties if too high) and thermal gradient (which creates stress inside the bearing). However, GFB geometry requires complex model of the structure [LEE 11] In 2006 Seghir-Ouali et al. [SEG 06] carried out an experimental identification technique for the convective heat transfer coefficient inside a rotating cylinder with an axial airflow in order to predict heat transfer more accurately.

Another way to act on heat transfers inside the bearing is to improve its design, like the bearing profile (working on a optimum clearance, or a specific conical geometry for instance), in order to obtain better heat evacuation and load capacity [MIH 10].

The need for a 3D ThermoHydroDynamic model

In the gas film, temperature has a direct or indirect influence on almost every important local parameter : viscosity, conductivity, density, pressure. As a result, it not only influences the lubricant flow but it also modifies the bearing global parameters such as load capacity and power loss, as well as the bearing structure.

The need for a 3D ThermoHydroDynamic model Nowadays, the large majority of the studies take into account the local temperature variation in order to improve the reliability of their model. However, we will see in more details in this chapter why a 3D-THD model is mandatory when describing two-phase flow lubricant behavior.

General considerations

GFBs are likely to work under very specific conditions, such as very high rotational speed (100 000 R.P.M.). This can result in high shear stress area, which means high power loss and heat generation. As we already said in the literature review, experimental data show that the heat generation is one of the main issues to deal with in order to improve GFB reliability and performance.

Specific considerations related to two-phase flow

In our case, the gas lubricant film is likely to contain a zone where the lubricant is a mixture of vapor and some small liquid particles. In this zone, the fluid properties as well as the power loss will be different from the vapor-phase characteristics. However the vapor pressure, which is settle the thermodynamical limit between the vapor and the mixture zone in the fluid film, strongly depends on the local temperature.

As a result, even if a global thermal balance gave a good estimation in the temperature rise in the bearing, its lack of precision would lead to great error when dealing with two-phase flow. Neglecting the local temperature variation is not possible and a THD model is mandatory.

A 3D THD model

In the case of GFBs, a THD model is obviously taking into account the variation in the circumferential direction which is the main direction of the fluid flow. It also has to take into account the radial direction variations because of the strong temperature gradients across the film thickness. These two directions (circumferential and radial) are integrated to most of the THD model. However, we saw in the literature review that most of the recent models also consider the axial direction variations. Again, strong axial temperature gradients are likely to occur because of the inner flow generated when pressure drops the ambient pressure, locally cooling down the fluid film and creating a "natural" lubricant supply zone.

In GFB theory, it is often assumed that there is no sub-ambient pressure zone, since a force equilibrium on the top foil shows it is compensated by an extra lubricant supply from the bearing edges. Even tough this behavior reduces the axial flow, it is still useful in order to account for axial thermal gradient (from the shaft surface for instance).

Influence of temperature variations on the model

Dimensionless parameters

In this section, we sum up the different characteristic dimensions and the related dimensionless parameters used in the energy equation :

Parameter name

Symbol Dimensionless form Heat capacity

c p cp = c p • 1/c p0 Temperature T T = T • 1/T 0 Thermal expansion coefficient α α = α • 1/α 0 Thermal conductivity k k = k • 1/k 0
Pe ρ cp ũ ∂ T ∂ x + ṽ∂ T ∂ ỹ + R a L 2 w ∂ T ∂z = α 0 T 0 N d Λ α T ũ∂ p ∂ x + R a L 2 w∂ p ∂z + ∂ ∂ ỹ k * ∂ T ∂ ỹ + N d μ * ∂ ũ ∂ ỹ 2 + ∂ w ∂ ỹ 2 (3.1)
With :

P e = ρ 0 c p0 R a ΩC 2 k 0 L (3.2) N d = µ 0 (RΩ) 2 T 0 k 0 (3.3) α = - 1 ρ ∂ρ ∂T p (3.4)
Where P e is the Peclet number and N d the dissipation number.

Hypothesis

In order to solve this equation, a similar procedure as for the Reynolds equation has been done. Taking into account the thin-film hypothesis, the negligible terms have been removed from the original equation. On the left hand side of the equation are the convective terms, while on the right hand side of the equation, the first term corresponds to the fluid dilatation, the second term corresponds to the turbulent conductive term and the last term correspond to the turbulent dissipative term.

In the fluid compressibility term, the volume expansivity at constant pressure α must be determined in order to solve the equation. The lubricant is supposed to be compressible and its density is a non-linear function of pressure and temperature. However, we assume here that it can be calculated using the ideal gas law as a first approximation, in order to be able to solve the equation numerically. Further study maybe required to improve this part of the model.

In that case, the density is proportional to the inverse of the temperature :

ρ = M p R 1 T (3.5)

ThermoHydroDynamic model for GFBs

And the thermal expansion is :

α = M p ρRT 2 (3.6)
Finally, with ρ 0 calculated as the density in the surrounding environment conditions, under the ideal gas law assumption, one find that :

ρ 0 = MP 0 RT 0 (3.7)
It leads to the following expression in the vapor phase (the thermal expansion coefficient is considered to be zero in the liquid phase) :

α = p ρ T 2 (3.8) α 0 = 1 T 0 (3.9)
The ideal gas assumption is valid for several reasons and in that particular case, it gives very good approximations far from the vapor pressure value. Besides, when close to the transition and in the transition zone, the compressibility terms become less important in magnitude, since it tends to a liquid behavior.

Boundary conditions for the lubricant film

Circumferential direction

In the circumferential direction, for the numerical treatment of the boundary condition, we must give the last node number m the same parameter values as node number 1 ( for which parameters values are calculated) for them to correspond to a single geometrical location.

As a consequence, the node after node number m in the rotational direction is node number 2, then 3, etc. and the node before node number 1 is node number m -1, then number m -2, etc. This condition is very important since the energy equation is a parabolic equation and with flow-oriented calculations.

If there is a groove in the bearing, then our model is able to describe the mix between the recirculation flow and the inlet flow. However, no specific predictive model concerning the groove could be integrated since the effect really depends on the groove geometry. Therefore we choose to use a model in which the fluid temperature at the groove is a function of the percentage of inlet flow, set by the user. When this percentage is zero, it corresponds to a bearing without a groove.

Radial direction

The radial (cross-film) direction is fundamental for the understanding of the heating process inside the bearing. Shear stress along this direction generates the major part of dissipation.

At the top foil, an adiabatic condition is considered. As we saw in the literature review, this assumption gives good results for a model which is dedicated to describe the lubricant behavior. It also allows us to test the model in the worst conditions in terms of running conditions.

Because of the shaft high rotational speed, we assume that at one axial location, the shaft surface temperature is constant along the whole circumference. We also assume that the shaft temperature is constant along the radial direction. As a result the heat transfer in the shaft is modeled as a 1D heat equation along the axial direction. In order not to be to restrictive, we set the temperature at the shaft edges equal to the surrounding temperature as boundary conditions. (It can easily be replaced by heat flux boundary conditions for a specific case study)

The no-slip condition and the steady-state ensure the continuity of the temperature at the lubricant/shaft interface. Therefore, the steady-state 1D heat equation in the shaft can be simplified as follow [LEE 10] :

k sha f t πR 2 a ∂T S ∂z 2 -k S 2πR a 0 ∂T S ∂y dx = 0 (3.10)
Where k sha f t is the shaft thermal conductivity, not to be mistaken for k S which is the lubricant thermal conductivity at the shaft surface.

Axial direction

Pressure gradient along the axial direction also has an effect on the temperature field. When the pressure gradient is negative, a part of the lubricant enters the bearing from its edges and helps the cooling process. However, the axial temperature gradient is detrimental since it generates stress inside the bearing structure which can be hazardous, especially for the thin metal foils in GFBs.

ThermoHydroDynamic model for GFBs

The temperature at the bearing edge is set according to the axial flow orientation. If the fluid enter the bearing, the temperature is the ambient temperature. Otherwise, the temperature is calculated thanks to a linear approximation from the two previous nodes (backward scheme).

Turbulence effect on heat generation and heat transfers

Concerning the theory of turbulence that gives a new dissipative term and a new conductive term, a simplified explanation is to say that there is an equivalent viscosity and an equivalent thermal conductivity in the equation which account for the effect of turbulence on the fluid local properties. (Refer to the chapter 4 for a detailed explanation) This theory is really convenient for our model. Once the equivalent parameters have been calculated thanks to the turbulent theory, the energy equation is solved exactly as it would be in laminar conditions. In fact, it is common that both laminar and turbulent conditions exist in the bearing at the same time, depending on the location inside the lubricant film. thanks to this method, the THD problem will be solved in both area without any distinction.

Viscosity variations

Temperature variations have a direct impact on lubricant viscosity. For gas, we used an explicit formulation which gives the molecular viscosity as a function of the temperature :

[SUT 93] μ = T 3/2 1 + Su T 0 T + Su T 0 (3.11)
The Sutherland number Su is a constant which depends on the lubricant composition, which value can be find in the literature. This formulation is valid from 0 to 555K and gives a good approximation at low pressure (under 30 bars the error is inferior to 10%).

Conclusion on the ThermoHydroDynamic model for Gas Foil Bearings

Conclusion on the ThermoHydroDynamic model for Gas Foil Bearings

In this chapter, we implemented an important part of our THD model. Solving the 3D thin-film energy equation gives a temperature field which leads to the temperature-related viscosity variations. This essential part of a THD model takes into account cross-film temperature variations as well as the influence of inlet lubricant. This will affect global parameters as well as the bearing structure itself.

The turbulent model is directly related to the other parts of the THD model, as it influences the pressure and density fields, as well as the thermal transfers. In the next chapter, we will describe the numerical methods that are able to solve the 3D turbulent THD problem, in which several parameters are coupled.

Chapitre 4

High-speed GFBs and turbulent flow

GFBs have the ability to run at very high speeds and therefore a turbulent regime in the lubricant is very likely to occur. Turbulence in lubrication can have a strong influence on global parameters such as load capacity or power loss. Turbulence is a very complex phenomena to study and it is essential to understand the main principles that describe properly the phenomena. Turbulent models are very numerous and the most difficult job is to find a realistic model easy to integrate in a GFB THD model. One approach to describe a turbulent flow is direct numerical simulation. The major asset of this technique is that no approximation is needed since all the energetic scales are taken into account. However, a higher Reynolds number means more scales to be accounted for. Therefore, the computational capacity is the limit to high Reynolds number direct numerical simulation. We will not use this technique in this work.

Prandtl mixing length

The Prandtl mixing length is a parameter that corresponds to the characteristic length of a fluid particle. It can also be seen as the characteristic distance that a fluid particle travels before mixing with other particles due to turbulent motion. Prandtl assumed that the eddy viscosity could be seen as a function of the average velocity gradient and the mixing length. The mixing length is not a fluid property : it depends on turbulent fluctuation, geometry, etc. Therefore, empirical values are needed. Some of first studies cited above as well as the most recent theories use (directly or indirectly) the mixing-length theory [PRA 26]. The ones which use directly the mixing length are based on a scalar parameter and therefore, they are called "0 equation" models. The several models using mixing-length theory in turbulent lubrication were trying to approximate the mixing-length value in order to give accurate viscous coefficients in the Reynolds equation.

Concerning turbulent flow in journal bearings, numerous authors have chosen some variations on this theory. It allows an accurate description of the influence of the turbulence on both local and global bearing parameters, without any extra equation to solve. In that case, the eddy viscosity only depends on the local fluid parameters and therefore, it is called a 0-equation model.

Statistical models for hydrodynamic studies

Another way to study turbulence is to use statistical models. In this kind of models, the turbulent flow is described by the (statistical) average Navier-Stokes equation. Under specific conditions (See turbulent transition in section 4.5), turbulent flow is likely to occur in GFBs. The early development of turbulent lubrication theory [CON 65, NG 65, TAY 69] allowed a clear understanding of the direct influence of turbulence in lubrication.

High-speed GFBs and turbulent flow

Generally speaking, it is well known that turbulence locally modifies fluid characteristics, whose variations can be modeled by adding an extra viscosity to the molecular viscosity, called "eddy viscosity". In the more specific case of journal bearings, the result is that the pressure field is modified when turbulence occurs in GFBs, as well as the thermal behavior [CON 72, SAF 73, SZE 80].

Statistical models and THD studies

Other authors applied the model in order to improve the thermal prediction in a journal bearing, solving for example both the GRE and the energy equation for turbulent flow for a finite-length journal bearings, thurst bearings or seals [NAG 94, ABE 95, BOU 96, BRU 01a, TAN 04, SHY 12], or used variable density and specific heat for the lubricant [CHU 04].

These studies are also influencing the bearing structure conception since it is directly linked to the flow behavior [ADO 10, NIC 10]. San Andres [SAN 95] also studied turbulence in journal bearing for cryogenic fluids.

In order to be consistent with a 3D THD model, a model accounting for turbulence in the three main directions is required. Some of the 0-equation models cited above do not account for a local variation of the turbulence in the radial direction, whereas the main speed gradient is across the film. A few studies have been carried out considering a turbulent model in which the radial dimension is directly involved thanks to improved 0-equation models or 2-equation models, such as the kε model [KAT 83, MAN 09, YU 10, XU 11, WAN 11a].

The next section gives a more extensive explanation of theses models as well as theoretical background in order to understand why we singled out one model among the others.

0-equation model for HD study

Several models have been proposed in order to create a bearing model under turbulent conditions. Most of these models are based on the mixing-length theory and empirical considerations. In fact, the objective is to link the shear stress to the pressure in the bearing for viscous thin films.

In laminar regime, under the assumptions listed below, the relations between the pres-sure and the shear stress are : Constantinescu used the mixing-length theory in order describe the total shear stress [CON 65]. Getting the velocity functions in the film, we can use the mass conservation equation (2.22) to derive an equation similar to the classic Reynolds equation (under no slip assumption) :

∂p ∂x =
∂ ∂ x ρh 3 μ 1 gc x ∂ p ∂ x + R a L 2 ∂ ∂z ρh 3 μ 1 gc z ∂ p ∂z = Λ 2 ∂ ρh ∂ x (4.7)
Where gc x and gc z are the Constantinescu's turbulent functions respectively along the x and z direction, which depend on the local Reynolds number.

A different approximation of the total shear stress is preferred in the Ng and Pan model [NG 65], based on the eddy viscosity and the Boussinesq's hypothesis. They use a law of wall to calculate the speed functions and they integrate them to find a turbulent Reynolds equation. This involves the same type of equation than those from Constantinescu, but 4. High-speed GFBs and turbulent flow with coefficients g x and g z which depend on several parameters : film thickness, local Reynolds number, axial and circumferential pressure gradients.

∂ ∂ x ρh 3 μ g x ∂ p ∂ x + R a L 2 ∂ ∂z ρh 3 μ g z ∂ p ∂z = Λ 2 ∂ ρh ∂ x (4.8)
The most current formulation for g x and g z is a combined result of the initial formulation from Constantinescu, the work of Ng and Pan on the coefficient variations as a function of the local Reynolds number, and finally some hypothesis on the relation between velocity and pressure gradients in the boundary layer :

g x =
1 12. + 0.0136ℜ 0.9 l (4.9)

g z = 1 12. + 0.0043ℜ 0.96 l (4.10)
With the local Couette Reynolds number based on the shaft diameter :

ℜ l = ρR a Ωh µ (4.11)
The main drawback with this method is that the turbulent coefficients are already the results of an integration through the film thickness, which is perfectly fine as long as one only needs to solve the Reynolds equation. However, it cannot be implemented in a THD model using a 3D energy equation because one needs to introduce the local influence of turbulence in the energy equation, in the three directions. That is the reason why we choose another method for the THD model.

Eddy viscosity theory 4.3.1 The Navier-Stokes equations

These are the main equations that describe the fluid behavior for incompressible flow. We made the assumption of incompressible flow for the turbulent model as a first approximation in order to simplify the whole turbulent THD model. Even if this hypothesis loses information regarding the fluid behavior, the validation of our results as well as the literature on turbulent flow in GFBs confirms that the hypothesis is possible without being too detrimental.

ρ ∂v i ∂t + v j ∂v i ∂x j = - ∂p ∂x i + ∂ ∂x j µ ∂v i ∂x j + ∂v j ∂x i (4.12)
Eddy viscosity theory

The average Navier-Stokes equations

In this chapter, < • > is for statistical average and • is for statistical fluctuation. In order to describe the turbulent flow using a statistical model, one has to isolate the perturbation to turbulence which modifies the flow in average. That is the reason why we apply a statistical average process to the Navier-Stokes equations :

ρ ∂ < v i > ∂t + < v j > ∂ < v i > ∂x j = - ∂ < p > ∂x i + ∂ ∂x j µ ∂ < v i > ∂x j + ∂ < v j > ∂x i -ρ < v i v j > (4.13)
As we can see, this operation allows us to isolate the influence of turbulence to a single fluctuating term -ρ < v i v j >.

Stress tensor

In turbulent flow, the average stress tensor τ is written as follow :

τ i j = µ ∂ < u i > ∂x j + ∂ < u j > ∂x i (4.14)
In turbulent flow, the average Navier-Stokes equations (Equation 4.13) lead to a turbulent contribution to stress which can be written as the Reynolds tensor τ t :

τ t i j = -ρ < u i u j > (4.15)
As result, the sum of the laminar and turbulent contributions gives the total stress :

τ * = τ + τ t (4.16)
However, additional equations are needed in order to find the additional unknowns and solve the turbulent problem. In order to refine the model, additional equations are needed to describe these turbulent terms and to "close" the problem. The aim here is to find an equilibrium between realism and complexity.

High-speed GFBs and turbulent flow

Numerous models exist and their complexity and ability to accurately describe a turbulent problem vary from one to another. Some of these models are suitable to our THD model since the equations related to turbulence can be solved simultaneously with the Reynolds equation, the non-linear EoS and the energy equation without being too time-consuming.

Eddy viscosity

The eddy viscosity µ t is based on the Boussinesq's hypothesis which implies that there is a linear relation between the Reynolds stress tensor and the average speed gradients, by analogy with a newtonian fluid :

τ t i j = -ρ < u i u j >= µ t ∂ < u i > ∂x j + ∂ < u j > ∂x i (4.17)
As we can see in this model, the fluctuation is not calculated. It is now a function of the average velocity gradient (obtained with the Reynolds equation) and the eddy viscosity. In this model, the eddy viscosity is unknown and one additional equation is needed in order to find it.

Even if the Boussinesq's hypothesis is in the most general case a strong approximation [SCH 07], it gives numerous turbulent model for bearings that show good results when compared to experimental data.

The molecular viscosity µ is an intrinsic property of the fluid, but the eddy viscosity µ t is not, since it is based on statistical fluctuations. Next, we will consider several means to accurately model the eddy viscosity, which is the only unknown left in our system of equation. Once it is modeled, the system is closed.

3D eddy viscosity model

We need a 0-equation turbulent model which gives the eddy viscosity as a 3D function of local fluid parameters, and can be used as so in a THD model. We choose a 3D eddy viscosity model in which the influence of the turbulence is calculated through a modified version of the Reichardt empirical law.

The empirical coefficients are based on shear stress and velocity measurements in a pipe flow [REI 51] and adapted to fit the experimental data [NG 64]. In that case, the eddy viscosity is given by the following formula :

|τ| = µ * ∂ < u > ∂y 2 + ∂ < w > ∂y 2 (4.22)
Where κ is the Von Karman constant, δ l is the thickness of the laminar sublayer and y + is the dimensionless distance from the wall.

0-equation model applied to 3D THD study

This model gives good results for bearing predictions and has the ability to give a 3D eddy viscosity. Thus, it accounts for turbulence in both the 3D energy equation and the GRE.

Turbulent Generalized Reynolds Equation

In turbulent conditions, the GRE demonstration must be adapted. Thanks to the Boussinesq approximation applied to the thin-film Navier equations, the relations between pressure gradients and shear stress gradients simplify to :

∂ p ∂ x = ∂ ∂ ỹ μt ∂ < u > ∂ ỹ + ∂ ∂ ỹ μt ∂ < u > ∂ ỹ (4.23) ∂ p ∂ ỹ ≈ 0 (4.24) ∂ p ∂z = ∂ ∂ ỹ μt ∂ < w > ∂ ỹ + ∂ ∂ ỹ μt ∂ < w > ∂ ỹ (4.25)
From this point, the GRE demonstration is the same, except that the molecular viscosity is replaced by the equivalent viscosity that accounts for turbulence. The dissipation term in the energy equation corresponds to the heat generation due to the shearing in the film. When in laminar regime, this heat generation is the result of viscous dissipation due to the mean velocity field. In turbulent regime, the turbulent kinetic energy is transfered from the large to the small scales where it is dissipated.

The turbulent dissipation also generates heat due to the fluctuating velocity field and can be modeled thanks to an equivalent eddy viscosity. Following an analogous process to the turbulent Reynolds equation, under the hypothesis of homogeneous turbulence and the Boussinesq hypothesis, we simplify the turbulent dissipation term in the energy equation [BRU 01b] :

µ ∂u i ∂x j ∂u i ∂x j + ∂u j ∂x i = µ t ∂ < u i > ∂x j ∂ < u i > ∂x j + ∂ < u j > ∂x i (4.26)
Therefore, the turbulent dissipation term can be written as a function of the mean velocity only (and not the fluctuation terms) and the equivalent viscosity µ * = µ+µ t as follow :

µ * ∂ < u i > ∂x j ∂ < u i > ∂x j + ∂ < u j > ∂x i (4.27)

Turbulent conduction

The conductive term in the energy equation corresponds to the heat transfer due to the conductivity of the fluid. When in laminar regime, the conductivity is a fluid intrinsic property and the conductive heat transfer is the result of the mean temperature gradient. In turbulent regime, the fluctuating velocity generates an extra heat transfer.

The fluctuating movement transfers kinetic energy as well as heat. An analogous process to the Boussinesq hypothesis leads to the introduction of an eddy diffusion coefficient. Then, the following approximation can be done :

-ρ < u i ∂T ∂x i >= µ t ∂ < T > ∂x i (4.28)
Therefore, the turbulent conductive term can be written as the sum of the laminar and turbulent conductive terms and again, only as a function of the mean temperature. If we Transition from laminar to turbulent flow set that k * = k + k t :

∂ ∂x i k * ∂ < T > ∂x i (4.29)
With :

k t = µ t k ρPr t α f (4.30)
Where α f is the heat transfer diffusivity :

α f = k ρc p (4.31)
The turbulent Prandtl number Pr t is the ratio between the momentum eddy diffusivity and the heat transfer eddy diffusivity. Following the Reynolds analogy it is often assumed in lubrication that Pr t = 1 [HIN 59].

Transition from laminar to turbulent flow 4.5.1 Characteristic numbers

In the case where the bearing is fully turbulent, the turbulent model detailed above will predict the eddy viscosity. However, the transition between the laminar and turbulent regime is a complex phenomena. Different regimes can exist simultaneously at different locations inside the bearing.

We saw that some of the formula used in turbulence theories use the Reynolds number (equation(4.11)) which is the ratio between the inertia forces and the viscosity forces. Basically, a bearing can be describe as a rotating cylinder (the shaft) inside a hollow cylinder (the housing).

The theory says that for two cylinders like this, Taylor vortices develop when the local Taylor number Ta l reaches the value Ta c , the limit between the laminar and Taylor vortices regimes. When Ta l reaches 2Ta c , the transition between the Taylor vortexes regime and the turbulent regime, the flow becomes turbulent.

The local Taylor number depends on the local Reynolds number and the film thickness over bearing radius ratio. It can be written as :

Ta l = ρR a Ωh µ h R a = ℜ l h R a (4.32)
For the critical Taylor number, we choose to use the empirical relation [FRE 90] :

Ta c = 63.3 • ε 2 r -38 • ε r + 41.2 (4.33)
The value 41.2 comes from the theory of the Taylor vortices and is the limit value when the eccentricity ratio ε r is zero. This theory supposes that the three regimes (laminar, Taylor vortices, turbulent) can co-exist at the same time at different location in the fluid-film. When the regime is laminar, the eddy viscosity is zero. When the regime is turbulent, we compute the eddy viscosity thanks to the formula (4.18). Between these two situations, we apply a sinusoidal transition for the eddy viscosity, depending on the local Taylor number :

μt = κ y + -δl tanh( y + δl ) • sin Ta l -Ta c Ta c • π 2 (4.34)
Recent works [DEN 08] have shown that the transition regime is still the subject of important research and that more complex models can be applied (modified Reynolds equation, dynamic viscosity) in order to refine the description of the transition. The sinusoidal transition we choose to apply here is a first step and further study may help to refine it.

Conclusion on high-speed GFBs and turbulent flow

In this chapter, we summed up the main theoretical turbulent models. Considering the turbulent 3D THD model requirements we selected an appropriate 0-equation model which calculates a 3D eddy viscosity for the GRE and the thin-film energy equation, as well as a turbulent thermal conductivity for the thin-film energy equation. We also presented the transition between the laminar and the turbulent regime, and implemented a numerical model describing the local turbulent behavior in the lubricant film.

The next chapter describe the numerical methods and process we use in order to solve the 3D turbulent THD problem.

Chapitre 5 Numerical methods

In the previous chapters, we described all the parameters relative to the 3D turbulent THD problem. Some of these parameters can be computed directly and other are given by numerical resolution. The fact that most of these parameters are coupled is a central issue. In this chapter, we describe the numerical methods we implemented in order to solve the whole problem.

Numerical methods

Literature review

Finite Element Method

The Finite Element Method (or FEM) is widely used in both solid and fluid study for hydrodynamic lubrication problems. It is particularly well suited for the housing behavior in TEHD problems, or for complex geometry which implies an unstructured mesh. [BOU 85, KUC 00, PIF 00, BRU 03, SHA 03, KUC 04]. However, implementing this method requires a relatively complex code compared to others.

Finite Volume Method

Nowadays, the Finite Volume Method (or FVM) has become very popular, not only in fluid mechanics but also in tribology [CER 98, WU 00, ALS 01, ARG 02]. This method relies on the Finite Difference Method (FDM) but it is an integral method that describes a domain made of elementary cells.

One of the major FVM assets is that it conserves properties across the cell boundaries. Some modified FVM schemes have been created in order to cope with the lack of efficiency of the "standard scheme" in particular cases (convergence issues when reverse flow occurs). Other techniques based on the FVM are also applied in bearing study like the "hybrid scheme" [HAT 02].

Finite Difference Method

FDM has been one of the first method used for solving hydrodynamic lubrication problems. The FDM theory is based on simple principles but it can solve rather complex THD problems in an efficient and accurate way. Therefore, a lot of recent studies use FDM to solve hydrodynamic or THD problems [CHA 02, BOU 03, COS 03, BRU 05, AHM 10].

Besides, there is no doubt that FDM is very convenient when working with simple geometries such as the plain journal bearing or GFB profiles. It is also able to treat complex problems like reverse flow, using appropriate schemes. Finally, FDM is relatively easy to implement.

Various possibilities

It is worth noticing that other methods exist and are employed in more specific cases. For example Computational Fluid Dynamics (or CFD) is very useful when one wants to study a bearing with groove(s), where the hypotheses upon which the Reynolds equation

Programming language

Some very efficient commercial numerical solver exist. However, there are two main issues regarding the use of such tools in our case. The first issue is that we need to cope with very specific problems regarding phase transition. The second issue is that this particular work is part of a company large research plan which goal is to integrate the computer routine to a predictive program for the whole turbo-machine. Therefore the commercial numerical solver is not appropriate in this particular case.

We did not use high-level programming environments such as Matlab either, even if they offer numerous helpful features, intrinsic functions, visual display, etc. First, this type of environment is too slow to match our computational requirements. Besides, the model integration would be less convenient compared to the solution we chose.

The program used for the study has been written in FORTRAN, a language widely used for its simplicity, rapidity and capacity to solve complex physical problems with numerical methods. Then, the integration of our bearing model to the company turbomachinery model, also written in FORTRAN, will be simplified.

Main calculation flow chart

Equation of State

The EoS is a 3 rd order polynomial form : Density is the unknown and the pressure is one of the parameters. We found out that an efficient and easy way to solve this equation is to find the roots of the polynomial thanks to the dichotomy method. This method is very convenient for this particular problem since the physics ensure there is only one solution for a range of pressure between 0 and the vapor pressure, which is the vapor density, and another solution between the vapor pressure and the critical pressure, which is the liquid density.

Therefore, for a given set of pressure and temperature, we can easily restrain the domain for density in order to find the corresponding root, either the vapor density or the liquid density. If f ( ρ) is the dimensionless polynomial, the EoS is considered solved in one point when its absolute value is inferior to a given criteria. Then, the corresponding value of ρ is the solution.

Eddy viscosity

One of the main advantage of the formula we use for the turbulence model is that it gives directly a value for the eddy viscosity. Therefore, no specific numerical methods are required for the turbulent model, saving time and memory for the other calculations. However, the non-linear THD problem is solved iteratively and the value of the viscosity must be calculated at each node and each iteration.

General Reynolds Equation

The GRE is solved using a Gauss-Seidel algorithm and a centered FDM discretization (adapted to this kind of elliptic equations) (Figure 5.1). Tests have shown that over-relaxation is suitable either to have an accelerated convergence in simple cases or to refine the solution in more complex situations. (for example, high rotation speed combined with high eccentricity ratio) The local relative error on pressure is computed at each point and each iteration as follow :

|p -p old | |p old | (5.1)
The GRE is considered solved when the maximum value is smaller than a given criteria.

Thin-film Energy Equation

The thin-film energy equation is solved with a Gauss-Seidel algorithm and a floworiented FDM discretization. In order to take into account the flow orientation inside the bearing, one has to use upward or backward scheme to improve the algorithm when using temperature derivatives, depending on whether the node will be influenced by the previous node or the following one. In order to cope with the small lack of accuracy of the algorithm at some points (mostly when a reverse flow occurs), we used a convergence criteria which combines local and global error.

The local relative error on temperature is computed at each point and each iteration as follow : The global relative error on temperature is computed at each iteration as follow :

∑ |T -T old | ∑ |T old | (5.3)
This algorithm is strong enough to give the solution of the thin-film energy equation even at very high-speed (100,000 R.P.M.) and high eccentricity ratio (0.9) which, combined together, are very severe conditions compared to the current industrial applications.

Conclusion on the Numerical methods

As we will see in the next chapter, the simultaneous resolution of all the equations using FDM with an appropriate scheme gives very good results in terms of strengh(High Peclet number, high Reynolds number, two-phase flow), precision and calculation time.

Chapitre 6

Results and discussion

In the previous chapter, we detailed different theoretical and numerical issues related to our model. In this chapter, we are going to confront the theory and the hypothesis with the results in order to analyse the reliability of the model. In this section, we investigate the accuracy of the numerical model regarding the relation between pressure and density, as well as the influence of the global temperature. In order to be able to compare our numerical model to Refprop, we run the model setting a constant global temperature in the film. These are not realistic running conditions but are perfectly adapted to the EoS validation.

We use a GFB which characteristics and running conditions are described in Table 6 This GFB is a standard model and we choose a high eccentricity ratio to generate a high pressure area. The lubricant is also a classic halogenated refrigerant, which is likely to be under two-phase flow conditions in this study.

Comparison with the Refprop database

We plot the density versus the pressure at bearing mid-length (because this is where we can get the larger range of pressure and density) for various temperatures, and compare the results from our simulation with REFPROP (Figure 6.1). In this graph, the pure liquid phase density from REFPROP is not displayed for readability reasons (its density is two order of magnitude higher than the liquid phase). Then continuous line

Results and discussion

shows the values given by REFPROP only for vapor phase whereas the doted line shows the ones given by our model in the vapor phase and in the mixture (not simulated by REFPROP).

Looking at the REFPROP data, there is a discontinuity at the vapor pressure value, where the density goes from the vapor to liquid phase value. The comparison on this figure confirms the good prediction in vapor phase, as well as the value of the vapor liquid transition. It also shows how the transition model allows a continuous transition between vapor and liquid.

Concerning the transition in our model, we observe the mixture density which starts above the vapor pressure and increases with a steep slope compared to the one in the vapor phase.

Ideal gas hypothesis validity

The results obtained with the non-linear EoS for isothermal hypothesis and far from the liquid/vapor transition are close to those obtained with the ideal gas law. These results confirm that the asymptotic trend of a compressible gas for relatively low pressures (i.e. far from the vapor pressure) can accurately be described with the ideal gas law.

Vapor phase lubricant under high pressure

What we consider to be "high pressure" here is a pressure which is close to the vapor pressure at a given temperature. When the bearing is running under these conditions, the lubricant can no longer be considered as an ideal gas and therefore, the ideal gas law does not match the accuracy requirements anymore. That's the reason why a non-linear thermodynamic model is useful : It describes the liquid/vapor transition and it improves the model accuracy in the vapor phase as well.

6.2 Lubricant behavior under severe running conditions

General hypothesis

The reader should note that we only discuss fluid-related results in this study. We do not claim to give an exhaustive explanation of the whole structural behavior. Studies regarding coupling between this fluid model and the foil bearing structure will be the subject of further works.

We investigate GFB performance with the characteristics and operating conditions described in Table 6.1.

The main purpose of this study is to describe the local behavior of fluid parameters and how they are coupled together. For this reason we make relatively simple but realistic assumptions concerning the lubricant film. This is an unavoidable step for a clear understanding of the bearing behavior. Among the major assumptions are the following :

-A rigid GFB. This is a possibly restrictive assumption, but we believe that a thorough understanding of the fluid behavior by itself is a necessary condition to understanding the overall GFB performance. The study of the structure (deformation, and heat transfer) will be presented in future publication(s).

-An imposed shaft temperature. In order to isolate the influence of specific parameters, this boundary condition is the most convenient one. Moreover, recent studies concerning GFBs show that shaft temperature control is one of the keys for proper thermal management. [SAN 10] -An adiabatic boundary condition at the top foil. Computing accurate boundary conditions for top foil requires a complex thermal model for the structure. However,

Results and discussion

this boundary condition is realistic and allows the validation of experimentally observed phenomena.

We limit this study to the influence of biphasic flow on local fluid parameters, since our goal is to analyze the coupling between pressure, density, equivalent viscosity and temperature. In order to get a good understanding of most of the interactions, the eccentricity ratio is fixed to 0.9 (a rather severe operating condition). The speed varies from 20000 R.P.M. up to 100000 R.P.M., which again is severe considering most industrial applications. The following range of test conditions is a good representation of fluid parameter overall behavior. In order to have a better understanding of how the formation of liquid in the lubricant film modifies the bearing behavior, we first localize the mixture zone. The size of the mixture area (Figure 6.2) increases with shaft rotational speed. Its location is around the minimum film thickness, where the pressure reaches its highest values.

Fraction of liquid

The fact that the mixture zone develops in the high pressure zone confirms that its study is important since this is where the hydrodynamic lift is generated. Therefore, variations of the lubricant properties in this region can have a strong influence on the bearing performances.

Finally, we observe that only a small fraction of the lubricant is in the liquid state (order of magnitude of a percent). It is a justification of some assumptions we made regarding the use of the energy equation in two-phase flow. At small speed (20000 R.P.M.) standard profile is found, with a peak of pressure around 180 • and no effect due to the two-phase flow. However, for higher rotational speeds, the combined effects of turbulence and two-phase flow lead to an unusual pressure profile. The effect of turbulence is expressed in the Reynolds equation by the equivalent viscosity which tends to increase the pressure for a given set of operating conditions when compared to laminar flow.

Pressure

Results and discussion

The effect of two-phase flow is more subtle. Since the vapor pressure is determined by the local temperature, the transition is clear. In the mixture phase, the peak of pressure is smaller than it would be for an hypothetical equivalent lubricant with no phase transition (using the ideal gas law for instance). In the mixture, the pressure tends to be close to the vapor pressure while the density strongly rises. The turbulence has no direct impact on density (Figure 6.4) as there is no viscosity term in the EoS. However, since the density is a function of pressure and temperature, turbulence and density are also strongly coupled. The direct coupling with pressure is evident on the density profile. In a zone around the minimum film thickness which can be as large as one sixth of the circumference, pressure reaches the vapor pressure and a mixture zone appears.

Density

An equivalent density takes into account the role of the liquid and therefore creates a large peak at minimum film thickness location. Density is also function of temperature since the vapor pressure is strongly influenced by temperature. That is why a 3D THD model is highly desirable in two-phase flow study. In two-phase flow, strong density gradients occur and we see that the mesh takes the transition phenomena into account correctly, since no obvious discontinuities are detected in the results. From a numerical point of view, a more refined mesh is necessary in order to cope with the strong variations in this zone. As we can see, the three regimes (laminar, transition, turbulent) can exist simultaneously inside the film. The general trend is that at high film thickness, the Reynolds number is high and the flow is turbulent whereas at small film thickness (where most of the hydrodynamic lift is generated), the Reynolds number is lower and the flow is laminar.

Influence of turbulence

The variation of molecular viscosity due to temperature is of the same order of magnitude (or less) than its initial value (set in the bearing environment). Compared to FIGURE 6.6: Taylor number at mid-length and mid-thickness (dimensionless) the molecular viscosity, the eddy viscosity increases up to one order of magnitude higher.

We conclude that, in turbulent conditions, eddy viscosity variations influence the flow more than the molecular viscosity in some of the zone in the lubricant. However, it is important to accurately compute the molecular viscosity since its most important variations occur at the most critical location (the mixture zone around the minimum film thickness location).

Temperature

Temperature field (Figure 6.8) is three-dimensional but most of the effects are visible on the temperature profile at mid-length. The following results are those for 100000 R.P.M. The same trend is observed for every rotational speed. From 0 • to about 60 • a smooth temperature rise is observed until a discontinuity in the temperature slope occurs, due to reverse flow generated by the pressure gradient. (Close to the top foil, Poiseuille flow dominates relative to Couette flow).

The temperature rise from this location to the beginning of the mixture area is quite standard. In the mixture area (from about 140 • to about 210 • ) there is strong variations FIGURE 6.7: Equivalent viscosity at mid-length and mid-thickness (dimensionless) in temperature which occurs because of the proportion of liquid it contains, and also because the fluid is compressed and then extended. From about 210 • to 270 • an axial pressure gradient allows lubricant to enter the bearing, which generates the drop in temperature. From 270 • to 360 • the heating process due to shear stress continues as the axial pressure gradient decreases.

3D THD analysis 6.3.1 Necessity of a finite-length bearing hypothesis

We want to confirm that the use of a 3D model is useful compared to one where a dimension is neglected. We cannot neglect the circumferential direction for obvious reasons.

The axial direction has often been neglected but as we saw in the literature review, the THD study requires much more focus on this direction. Several conclusions regarding the importance of the axial direction in our model are listed below :

-The pressure field (Figure 6.9) is clearly depending on the axial location and the high pressure area is centered on the mid-length location. In the high pressure area in almost one forth of the zone at the bearing edges, the pressure increase is at least 50% smaller than at mid-length.

-A similar phenomena is observed for density in two-phase flow (Figure 6.10). The mixture zone is only located at the center of the bearing but again not at its edges where the pressure increase does not reach the vapor pressure.

-On Figure 6.11, we also see thanks to the streamlines that the bearing edges are a zone of lubricant exchange with the environment (both inlet and outlet flow). This exchange has a tremendous influence on the temperature field since the outlet flow drains off some of the heat generated in the lubricant whereas the inlet flow helps cooling it down.

-The last reason is related to the fact that in some turbo-machineries, some GFBs are likely to be mounted in such conditions where the environment pressure varies from one edge to another.

Lubricant behavior in the cross-film direction

The goal of the THD problem is to accurately describe the phenomena related to energy dissipation and heat generation. We see on figure 6.12 that the temperature distribution at mid-length is very complex and that strong radial temperature gradients appear. This is the result of shear stress and lubricant compressibility.

Close to the shaft surface, the Couette flow prevails. Therefore the flow is mostly unidirectional and along the circumferential direction. Because of the high shaft rotational speed, the temperature is constant along this direction. When going from the shaft surface to the top foil, the Couette flow influence decreases and a higher percentage of the fluid velocity is due to the Poiseuille flow. Some reverse flow occurs when the pressure gradient is strong enough. This reverse flow generates computational problems that were discussed in the numerical method chapter. From a physical point of view, we can observe that the recirculation is likely to create a high-temperature zone since the cooling becomes less effective in this situation. flow, the eddy viscosity can be one order of magnitude higher than the molecular viscosity and that is the reason why turbulence cannot be ignored. The film thickness plays a key role in the turbulent flow. As the eddy viscosity increases with the film thickness, most of the time the influence of turbulence is small in the high pressure area where the film thickness is small, but its overall effect plays a key role in energy transfers.

We can correlate the variation of the dimensionless distance from the wall y + (Figure 6.13) with the eddy viscosity. The high turbulent energy area is located far from the rotor and top foil surfaces, and we visualize the importance of the near-wall effect in a FIGURE 6.12: Temperature (K) and streamlines at mid-length at 100000 R.P.M. large part of the film (particularly in the minimum film thickness area where the flow is laminar). This is a particularity of the thin-film turbulent flow, where a large part of the flow is strongly affected by the vicinity of the walls.

We can underline that when y + is too small even at mid-film thickness, the criterion used to determine whether or not the flow is turbulent (with the Taylor number) correctly indicates that the flow is laminar.

Global analysis 6.4.1 Validation for the ideal gas isothermal model

For validation, we run the THD model setting a constant temperature in the film and using a dimensionless density proportional to pressure instead of the non-linear EoS in order to compare it to the current program for the same running conditions. This comparison was made for a large range of various running conditions. It appears that the local difference between the values given by the developed model and the former one (an internal program developed under constant viscosity and temperature assumptions) This process allowed us to validate the program on comparable running conditions. Therefore the predictive results can be given much more credit.

A predictive 3D THD model

If we compare the lift force and the power loss given by the former isothermal single-flow model and the THD two-phase flow model (Figure 6.15 and Figure 6.16), we see that at low speeds, we get really close prediction. This is part of the previous validation process since the pressure increase is not great enough to reach the vapor pressure.

On the other hand, great discrepancy in the results occurs as soon as the speed is high enough so that the vapor pressure is reached at some point. Then, we clearly see the influence of two-phase flow on global parameters. It is interesting to see that the appearance of a mixture zone causes a loss in lift force and less power loss when the speed is higher than 40000 R.P.M. This is due to the fact that the pressure in this zone is close to the vapor pressure.

Conclusion and prospects

In order to extend the field of applications for GFBs as well as there reliability, we studied their behavior in refrigerating gas. This is part of the ECS technology enhancement issue and we needed to consider complex phenomena in order to keep the main benefits of GFBs in these particular conditions. The fact that the refrigerating gas is working as a lubricant is part of the willingness to avoid external lubrication system.

The very high rotation speeds are also a major aspect in this field since it is strongly coupled with thermal management issues.

In this study, we showed the importance of an accurate description of the film parameters, which variations largely influence the bearing behavior. Among the principal theories, there are :

-Compressible lubricant, with an appropriate non-linear behavior when close to the vapor/liquid transition.

-Vapor/liquid transition and calculation of the mixture equivalent parameters.

-Turbulent flow for high-speed GFBs with a 3D eddy viscosity model.

-A 3D behavior for viscosity, particularly the cross-film variations. (temperature dependent)

-A 3D behavior for temperature, particularly in cross-film direction in order to be consistent with viscosity, but also in the axial direction in order to account for potential temperature gradient which considerably modifies the bearing 3D temperature profile.

The extend of the study includes a numerical model for the lubricant film behavior, but the consequences and implications go far beyond. The THD behavior of the film is very likely to have a strong influence on the bearing thermal structural behavior. The next step is to couple the lubricant THD model to a structural thermal model which is able to enhance thermal boundary conditions as well as to predict a map of the temperature at critical points in the structure.

An other interesting development is the integration of the flexible structure behavior for GFB types. This aspect of the study is complementary to an accurate prediction of the pressure field in the film. Once again, the precision of the lubricant model is mandatory in order to predict high stress zones in the structure which might damage the foils and be the cause of the bearing failure in the worst case scenario.

The first step will be a Von Mises stress calculation and on the long term, will be able to integrate an even more elaborated structural model in order to have a turbulent TEHD model. The structural dilatation aspect might be considered as well.

In parallel to further improvements, we will also adapt our model to thrust bearings. The fundamental difference will come from the geometry which will probably reveal new theoretical issues such as inertia effects.

We are also expecting to carry out an experimental campaign on a test rig which will consist in GFBs running in a close environment containing refrigerating gas. Therefore we are expecting a direct experimental two-phase flow validation.

Generally speaking, all these model enhancements will be integrated to the ECS predictive model in order to improve its overall efficiency.

[BRO 00] BROWN R., DRUMMOND G., ADDISON P.

Chaotic response of a short journal bearing.

FIGURE 2 :

 2 FIGURE 2: Gas foil bearings

FIGURE 3 :

 3 FIGURE 3: Air Conditioning Pack

FIGURE 4 :

 4 FIGURE 4: Air Conditioning Pack principle

FIGURE 5 :

 5 FIGURE 5: Air Cycle Machine

FIGURE 6 :

 6 FIGURE 6: Air Cycle Machine Principle

FIGURE 7 :

 7 FIGURE 7: Gas foil bearings

Sommaire 1 . 1

 11 GFBs characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 The early development of hydrodynamic lubrication . . . . . . . . . . . 12 1.3 Effect of the temperature on GFBs . . . . . . . . . . . . . . . . . . . . . 13 1.4 Predictive GFB models . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

FIGURE 2 . 1 :

 21 FIGURE 2.1: Comparison between the EoS (continuous lines) and REFPROP (doted lines)

FIGURE 3. 1 :

 1 FIGURE 3.1: Numerical discretization of the lubricant film

0

  -equation model applied to 3D THD study μt = κ y + -δl tanh(

4. 4 .

 4 2 Turbulent 3D thin-film energy equation 4.4.2.1 Turbulent dissipation

Sommaire 5 . 1

 51 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . 5.1.3 Finite Difference Method . . . . . . . . . . . . . . . . . . . . . . . 5.1.4 Various possibilities . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Programming language . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Main calculation flow chart . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Eddy viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.3 General Reynolds Equation . . . . . . . . . . . . . . . . . . . . . 5.3.4 Thin-film Energy Equation . . . . . . . . . . . . . . . . . . . . . . 5.4 Conclusion on the Numerical methods . . . . . . . . . . . . . . . . . . .

  Programming language is based become irrelevant [NAS 02, KOS 04, ISH 06]. Schumack [SCH 96] uses a pseudospectral method to solve a THD problem. Fatu et al. [FAT 06] also developed a new heat flux conservation algorithm to solve a complex big-end connecting-rod bearing TEHD problem.

  Conclusion on the Numerical methods|T -T old ||T old | (5.2)

FIGURE 6 . 1 :

 61 FIGURE 6.1: Validation of the EoS using Refprop

FIGURE 6 . 2 :

 62 FIGURE 6.2: Fraction of liquid at mid-length (dimensionless)

FIGURE 6 . 3 :

 63 FIGURE 6.3: Pressure at mid-length (dimensionless)

FIGURE 6 . 4 :

 64 FIGURE 6.4: Density at mid-length and top foil (dimensionless)

FIGURE 6 . 5 :

 65 FIGURE 6.5: Reynolds number at mid-length and mid-thickness (dimensionless)

FIGURE 6 . 8 :

 68 FIGURE 6.8: Temperature at mid-length and top-foil (dimensionless)

FIGURE 6 . 9 :

 69 FIGURE 6.9: Pressure field at 100000 R.P.M. (bar)

FIGURE 6 .

 6 FIGURE 6.10: Density at 100000 R.P.M. (dimensionless)

Figure 6 .

 6 Figure 6.14 displays the equivalent viscosity distribution. On this type of turbulent

FIGURE 6 .

 6 FIGURE 6.11: Temperature (K) and streamlines at mid-film thickness at 100000 R.P.M.

FIGURE 6 .

 6 FIGURE 6.13: Distance to the wall (dimensionless) at mid-length at 100000 R.P.M.

FIGURE 6 .

 6 FIGURE 6.14: Equivalent viscosity (dimensionless) at mid-length at 100000 R.P.M.

FIGURE 6 .

 6 FIGURE 6.15: Predicted lift force in two-phase flow

FIGURE 6 .

 6 FIGURE 6.16: Predicted power loss in two-phase flow

  GFBs characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 The early development of hydrodynamic lubrication . . . . . . . . . . . . 12 1.3 Effect of the temperature on GFBs . . . . . . . . . . . . . . . . . . . . . 13 1.4 Predictive GFB models . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 The Reynolds equation . . . . . . . . . . . . . . . . . . . . . . . 2.4 From the isothermal laminar model to the THD turbulent model . . . . . 2.4.1 A Generalized Reynolds Equation . . . . . . . . . . . . . . . . . 2.4.2 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . 2.4.3 Model improvements . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Conclusion on the hydrodynamic study and two-phase flow . . . . . . . . Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Lubricant film thermal study . . . . . . . . . . . . . . . . . . . . 3.1.2 Shaft thermal study . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Housing thermal study . . . . . . . . . . . . . . . . . . . . . . . Viscosity variations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion on the ThermoHydroDynamic model for Gas Foil Bearings . Statistical models and THD studies . . . . . . . . . . . . . . . . 4.2 0-equation model for HD study . . . . . . . . . . . . . . . . . . . . . . . 4.3 Eddy viscosity theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 The Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . 4.3.2 The average Navier-Stokes equations . . . . . . . . . . . . . . . 4.3.3 Stress tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Eddy viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.5 3D eddy viscosity model . . . . . . . . . . . . . . . . . . . . . . Turbulent 3D thin-film energy equation . . . . . . . . . . . . . . 58 4.5 Transition from laminar to turbulent flow . . . . . . . . . . . . . . . . . . 59 4.5.1 Characteristic numbers . . . . . . . . . . . . . . . . . . . . . . . 59 4.6 Conclusion on high-speed GFBs and turbulent flow . . . . . . . . . . . . 60 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1.1 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . 62 5.1.2 Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . 62 5.1.3 Finite Difference Method . . . . . . . . . . . . . . . . . . . . . . 62 5.1.4 Various possibilities . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2 Programming language . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3 Main calculation flow chart . . . . . . . . . . . . . . . . . . . . . . . . . 63 Simplification in the model . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Numerical discretization of the lubricant film . . . . . . . . . . . . . . .

	Table des matières
	Table des matières Table des figures Liste des tableaux Introduction Limiting assumptions 1 General considerations about gas foil bearings 1.1 2 Hydrodynamic study and two-phase flow 2.1 Literature review . 20 xi xv xvii 3 9 11 15 2.2.3 A modified Peng-Robinson equation of state . . . . . . . . . . . 21 2.2.4 Theoretical model validation . . . . . . . . . . . . . . . . . . . . 22 2.2.5 The Clapeyron formula . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.6 The vapor/liquid transition issue . . . . . . . . . . . . . . . . . . 23 2.2.7 A vapor/liquid transition model . . . . . . . . . . . . . . . . . . 24 2.3 ThermoHydroDynamic model for GFBs ρ Fluid mass density (kg.m -3 ) τ Stress tensor (kg.m -1 .s -2 ) a ρ EoS temperature-dependent coefficient Pr Prandtl number R Ideal gas constant (J.mol -1 .K -1 ) R a Bearing radius (m) 4.4.2 5 Numerical methods 61 Table des figures Liste des tableaux a min Minimal speed of sound in the mixture (m.s -1 ) T Local temperature (K) Nomenclature a V Minimal speed of sound in the vapor (m.s -1 ) T c Critical temperature (K) 5.1 69 6.1.3 Ideal gas hypothesis validity . . . . . . . . . . . . . . . . . . . . 70 6.1.4 Vapor phase lubricant under high pressure . . . . . . . . . . . . . 71 6.2 Lubricant behavior under severe running conditions . . . . . . . . . . . . 71 6.2.1 General hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.2.2 Fraction of liquid . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.2.3 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2.4 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.2.5 Influence of turbulence . . . . . . . . . . . . . . . . . . . . . . . 75 6.2.6 Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.3 3D THD analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.3.1 Necessity of a finite-length bearing hypothesis . . . . . . . . . . 77 6.3.2 Lubricant behavior in the cross-film direction . . . . . . . . . . . 79 6.4 Global analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.4.1 Validation for the ideal gas isothermal model . . . . . . . . . . . 83 6.4.2 A predictive 3D THD model . . . . . . . . . . . . . . . . . . . . 85 6.5 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 86 Conclusion and prospects 87 Bibliographie 89 1 Air Conditioning Pack . . . . . . . . . . . . . . . . . . . . . . . . . . . i 2 Gas foil bearings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii 3 Air Conditioning Pack . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Air Conditioning Pack principle . . . . . . . . . . . . . . . . . . . . . . 5 Air Cycle Machine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Air Cycle Machine Principle . . . . . . . . . . . . . . . . . . . . . . . . 7 Gas foil bearings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Comparison between the EoS (continuous lines) and REFPROP (doted lines) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 5.1 Pressure flow chart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Validation of the EoS using Refprop . . . . . . . . . . . . . . . . . . . . 6.2 Fraction of liquid at mid-length (dimensionless) . . . . . . . . . . . . . . 6.3 Pressure at mid-length (dimensionless) . . . . . . . . . . . . . . . . . . . 6.4 Density at mid-length and top foil (dimensionless) . . . . . . . . . . . . . 6.5 Reynolds number at mid-length and mid-thickness (dimensionless) . . . . 6.6 Taylor number at mid-length and mid-thickness (dimensionless) . . . . . 6.7 Equivalent viscosity at mid-length and mid-thickness (dimensionless) . . 6.8 Temperature at mid-length and top-foil (dimensionless) . . . . . . . . . . 6.9 Pressure field at 100000 R.P.M. (bar) . . . . . . . . . . . . . . . . . . . . 6.10 Density at 100000 R.P.M. (dimensionless) . . . . . . . . . . . . . . . . . 6.11 Temperature (K) and streamlines at mid-film thickness at 100000 R.P.M. . 6.12 Temperature (K) and streamlines at mid-length at 100000 R.P.M. . . . . . 6.13 Distance to the wall (dimensionless) at mid-length at 100000 R.P.M. . . . 6.14 Equivalent viscosity (dimensionless) at mid-length at 100000 R.P.M. . . . 6.15 Predicted lift force in two-phase flow . . . . . . . . . . . . . . . . . . . . 6.16 Predicted power loss in two-phase flow . . . . . . . . . . . . . . . . . . . 2.1 Dimensionless fluid parameters for the Reynolds equation . . . . . . . . 25 2.2 Dimensionless fluid parameters for the GRE . . . . . . . . . . . . . . . . 28 3.1 Dimensionless fluid parameters for the energy equation . . . . . . . . . . 42 6.1 Test GFB characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . 69 < • > Subscript for statistical average α Volume expansivity at constant pressure (K -1 ) α f Heat transfer diffusivity • Subscript for statistical fluctuation • * Exponent for the sum of the laminar and turbulent parameter values • t Exponent for the turbulent regime • 0 Subscript for the reference value • F Subscript for the top foil • L Subscript for the liquid phase • S Subscript for the shaft • V Subscript for the vapor phase • x Subscript for the circumferential coordinate • y Subscript for the radial coordinate • z Subscript for the axial coordinate ∆. Difference between vapor and liquid values of a given variable δ l Thickness of the laminar sublayer (m) ε r Eccentricity ratio γ Adiabatic index κ Von Karman constant Λ Bearing number µ Dynamic viscosity (Pa.s) ω ρ EoS acentric factor ℜ l Local Reynolds number A w Van der Waals coefficient for molecular interaction B EoS constant coefficient (for a given fluid) B w Van der Waals coefficient for volume exclusion C Bearing clearance (m) C . Phenomenological constants in Clapeyron's formula c p Heat capacity (J.kg -1 .K -1 ) D Bearing diameter (m) F 0 Viscosity integral coefficient F 0 (m 2 .s.kg -1 ) F 2 Viscosity integral coefficient F 2 (m.s) F 3 Viscosity integral coefficient F 3 (s) g Ng-Pan's turbulent functions G 1 Viscosity integral coefficient G 0 (m.s) G 2 Viscosity integral coefficient G 2 (s) gc Constantinescu's turbulent functions h Film thickness (m) k Thermal conductivity (W.m -1 .K -1 ) L Bearing axial length (m) M Molar mass (kg.mol -1 ) m ρ EoS constant coefficient (for a given fluid) N d Dissipation number p Pressure (Pa) P c Critical pressure (Pa) p sat Vapor pressure (bar) Pe Peclet number Ta c Critical Taylor number Ta l Local Taylor number u,v,w Circumferential, radial and axial velocity components (m.s -1 ) x,y,z Circumferential, radial and axial coordinates (m) y + Dimensionless distance from the wall 3.1 xii xiii xv xvii xix xx xxi

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 Modern lubrication theory basis . . . . . . . . . . . . . . . . . . 17 2.1.2 The short bearing approximation . . . . . . . . . . . . . . . . . . 17 2.1.3 The need for numerical methods . . . . . . . . . . . . . . . . . . 18 2.1.4 The Reynolds model . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.5 A Generalized Reynolds Equation (GRE) . . . . . . . . . . . . . 19 2.1.6 Thermodynamic equation of state (EoS) . . . . . . . . . . . . . . 19 2.2 Thermodynamic model for two-phase flow -A non-linear equation of state 20 2.2.1 Symbol description . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.2 Choosing the equation of state . . . . . . . . . . . . . . . . . . . 3 Behavior of refrigerant in laminar isothermal biphasic running conditions 25 2.3.1 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . 25 xi Table des matières 2.3.2 3.1.4 3-D energy equation models . . . . . . . . . . . . . . . . . . . . 3.1.5 Relations between bearing parameters and thermal behavior . . . 3.1.6 Thermal management . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The need for a 3D ThermoHydroDynamic model . . . . . . . . . . . . . 3.2.1 General considerations . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Specific considerations related to two-phase flow . . . . . . . . . 3.2.3 A 3D THD model . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Influence of temperature variations on the model . . . . . . . . . . . . . 3.3.1 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . 3.3.2 3D turbulent thin-film energy equation . . . . . . . . . . . . . . . 3.3.3 Hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Boundary conditions for the lubricant film . . . . . . . . . . . . . 3.3.5 Turbulence effect on heat generation and heat transfers . . . . . . 3.4 4 High-speed GFBs and turbulent flow 4.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Direct numerical simulation . . . . . . . . . . . . . . . . . . . . 4.1.2 Prandtl mixing length . . . . . . . . . . . . . . . . . . . . . . . . 4.1.3 Statistical models for hydrodynamic studies . . . . . . . . . . . . 4.1.4 4.4 0-equation model applied to 3D THD study . . . . . . . . . . . . . . . . 4.4.1 Turbulent Generalized Reynolds Equation . . . . . . . . . . . . . 5.3.1 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3.2 Eddy viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.3.3 General Reynolds Equation . . . . . . . . . . . . . . . . . . . . 64 5.3.4 Thin-film Energy Equation . . . . . . . . . . . . . . . . . . . . . 64 5.4 Conclusion on the Numerical methods . . . . . . . . . . . . . . . . . . . 65 6 Results and discussion 67 6.1 Thermodynamic behavior validation . . . . . . . . . . . . . . . . . . . . 69 6.1.1 Validation conditions . . . . . . . . . . . . . . . . . . . . . . . . 69 6.1.2 Comparison with the Refprop database . . . . . . . . . . . . . .

  Modern lubrication theory basis . . . . . . . . . . . . . . . . . . . 2.1.2 The short bearing approximation . . . . . . . . . . . . . . . . . . . 2.1.3 The need for numerical methods . . . . . . . . . . . . . . . . . . . 2.1.4 The Reynolds model . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.5 A Generalized Reynolds Equation (GRE) . . . . . . . . . . . . . . 2.1.6 Thermodynamic equation of state (EoS) . . . . . . . . . . . . . . . 2.2 Thermodynamic model for two-phase flow -A non-linear equation of state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Symbol description . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Choosing the equation of state . . . . . . . . . . . . . . . . . . . . 2.2.3 A modified Peng-Robinson equation of state . . . . . . . . . . . . Conclusion on the hydrodynamic study and two-phase flow . . . . . . .

	Literature review
	2.1 Literature review
	2.1.1 Modern lubrication theory basis

Sommaire 2.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 2.2.4 Theoretical model validation . . . . . . . . . . . . . . . . . . . . . 2.2.5 The Clapeyron formula . . . . . . . . . . . . . . . . . . . . . . . . 2.2.6 The vapor/liquid transition issue . . . . . . . . . . . . . . . . . . . 2.2.7 A vapor/liquid transition model . . . . . . . . . . . . . . . . . . . 2.3 Behavior of refrigerant in laminar isothermal biphasic running conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . . 2.3.2 The Reynolds equation . . . . . . . . . . . . . . . . . . . . . . . . 2.4 From the isothermal laminar model to the THD turbulent model . . . . 2.4.1 A Generalized Reynolds Equation . . . . . . . . . . . . . . . . . . 2.4.2 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Model improvements . . . . . . . . . . . . . . . . . . . . . . . . . 2.5

  The short bearing approximation can be found in recent theoretical works carried out by Reddy et al. [SUD 97], Brown et al. [BRO 00], Yang et al. [YAN 01], Keogh and Khonsari [KEO 01], Souchet et al. [SOU 04], Taylor [TAY 04], Song et al. [SON 05] and Ighil et al. [IGH 08] to name but a few. We can also cite Reason and Narang [REA 82] who combine long and short bearing approximations in order to describe the whole range of bearing width. The short bearing approximation is also useful when trying to create a simplified predictive model, such as Naïmi's low Sommerfeld number model [NA 10].

TABLE 2

 2

.1: Dimensionless fluid parameters for the Reynolds equation

  Shaft thermal study . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.1.3 Housing thermal study . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1.4 3-D energy equation models . . . . . . . . . . . . . . . . . . . . . 38 3.1.5 Relations between bearing parameters and thermal behavior . . . . 39 3.1.6 Thermal management . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 The need for a 3D ThermoHydroDynamic model . . . . . . . . . . . . . 40 3.2.1 General considerations . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.2 Specific considerations related to two-phase flow . . . . . . . . . . 41 3.2.3 A 3D THD model . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3 Influence of temperature variations on the model . . . . . . . . . . . . . 42 3.3.1 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . . 42 3.3.2 3D turbulent thin-film energy equation . . . . . . . . . . . . . . . . 3.3.3 Hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Boundary conditions for the lubricant film . . . . . . . . . . . . . . 3.3.5 Turbulence effect on heat generation and heat transfers . . . . . . . 3.4 Viscosity variations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion on the ThermoHydroDynamic model for Gas Foil Bearings

Sommaire 3.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1.1 Lubricant film thermal study . . . . . . . . . . . . . . . . . . . . . 33 3.1.2

  3-D energy equation models in order not to neglect the axial thermal gradient [BUK 06, CAI 06, PEN 06, AHM 07, SIM 08, SIN 08, SAN 09]. In 2006, Hannnon [HAN 06] proposes a Generalized Universal Reynolds Equation for variable fluid properties and variable geometry for self-acting bearings. However we do not need to use such a general approach since our study focus on common geometries. Dobrica and Fillon [DOB 06] confirm that 3-D algorithm is mandatory if one wants to predict axial thermal behavior accurately. Michaud et al. [MIC 07] find that minimum film thickness is overestimated in a 2-D model. Furthermore, thermal elastic deformation analyses require temperature variations in the three directions [FIL 08, MAO 08], especially for more complex geometry such as pads [GLA 06].

TABLE 3

 3

.1: Dimensionless fluid parameters for the energy equation 3.3.2 3D turbulent thin-film energy equation

  Even though an analogous analytical relation can not be proposed in turbulent regime, empirical coefficients can link the total shear stress to the pressure gradient.

				∂τ xy ∂y	(4.1)
		∂p ∂y	≈ 0		(4.2)
		∂p ∂z	=	∂τ zy ∂y	(4.3)
	The hypothesis needed for these relations are :
	-Isothermal flow
	-Thin lubricant film
	-Continuous environment
	-External mass forces negligible compared to viscous forces
	-Negligible inertia forces
	∂p ∂x	=	∂τ xy ∂y	+	∂τ t xy ∂y	(4.4)
	∂p ∂y	≈ 0				(4.5)
	∂p ∂z	=	∂τ zy ∂y	+	∂τ t zy ∂y	(4.6)

  Vapor phase lubricant under high pressure . . . . . . . . . . . . . . 6.2 Lubricant behavior under severe running conditions . . . . . . . . . . . 6.2.1 General hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Fraction of liquid . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.3 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.4 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.5 Influence of turbulence . . . . . . . . . . . . . . . . . . . . . . . . 6.2.6 Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 3D THD analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	6.1 Thermodynamic behavior validation
	6.1.1 Validation conditions

Sommaire 6.1 Thermodynamic behavior validation . . . . . . . . . . . . . . . . . . . . 6.1.1 Validation conditions . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.2 Comparison with the Refprop database . . . . . . . . . . . . . . . 6.1.3 Ideal gas hypothesis validity . . . . . . . . . . . . . . . . . . . . . 6.1.4 6.3.1 Necessity of a finite-length bearing hypothesis . . . . . . . . . . . 6.3.2 Lubricant behavior in the cross-film direction . . . . . . . . . . . .

  .1.

	Characteristic	Value Dimension
	Bearing
	Length	27	mm
	Shaft diameter	28	mm
	Clearance	90	µm
	Eccentricity ratio	0.9	-
	Shaft speed	100000	R.P.M.
	Lubricant
	Pressure	1	Bar
	Temperature	20	• C
	Name	R245fa	-
	Viscosity	12.3	µPa.s
	Molar mass	134.05	g.mol -1
	Heat capacity	976.9 J.kg -1 .K -1
	Thermal conductivity 0.012 W.m -1 .K -1

TABLE 6 .

 6

1: Test GFB characteristics

  Proceedings of the Institution of Mechanical Engineers, Part J : Journal of Engineering Tribology, vol. 214, n o 4, 2000, p. 387-400, Prof Eng Publishing. [BRU 01a] BRUNETIERE N. Etude théorique et expérimentale du comportement thermohydrodynamique des garnitures d'étanchéité. Thse de doctorat, Thèse de Doctorat de lUniversité de Poitiers, France, 2001. [BRU 01b] BRUNETIERE N., TOURNERIE B., FR ÊNE J. Etudes TEHD des garnitures mcaniques. XVme congrs franais de mcanique, 2001. [BRU 03] BRUNETIERE N., TOURNERIE B., FRENE J. TEHD Lubrication of Mechanical Face Seals in Stable Tracking Mode : Part 1 -Numerical Model and Experiments. Journal of Tribology, vol. 125, 2003, page 608. [BRU 05] BRUCKNER R., DELLACORTE C., PRAHL J. Analytic Modeling of the Hydrodynamic, Thermal, and Structural Behavior of Foil Thrust Bearings. rapport, 2005, NASA/TM2005-213811. [BUK 06] BUKOVNIK S., D ÖRR N., CAIKA V., BARTZ W. J., LOIBNEGGER B. Analysis of diverse simulation models for combustion engine journal bearings and the influence of oil condition. Tribology International, vol. 39, n o 8, 2006, p. 820 -826. 4th AIMETA International Tribology Conference. [CAI 06] CAIKA V., BUKOVNIK S., OFFNER G., BARTZ W. Elasto-hydrondynamic journal bearing model with pressure, temperature and shear rate dependent viscosity. AITC-AIT 2006 International Conference on Tribology, 2006. [CER 98] CERMAK J. A control volume-based discretization of the Reynolds equation for the numerical solution of elastohydrodynamic lubrication problems. International Journal for Numerical Methods in Fluids, vol. 26, n o 8, 1998, p. 977-986, Wiley Online Library. [CHA 02] CHANG Q., YANG P., MENG Y., WEN S. Thermoelastohydrodynamic analysis of the static performance of tilting-pad journal bearings with the Newton-Raphson method. Tribology International, vol. 35, n o 4, 2002, p. 225 -234. [CHR 41] CHRISTOPHERSON D. A new mathematical method for the solution of film lubrication problems. AR-CHIVE : Proceedings of the Institution of Mechanical Engineers 1847-1982 (vols 1-196), vol. 146, n o 1941, 1941, p. 126-135, Prof Eng Publishing. [CHU 04] CHUN S. Thermohydrodynamic lubrication analysis of high-speed journal bearing considering variable density and variable specific heat. Tribology International, vol. 37, n o 5, 2004, p. 405-413, Elsevier. [CON 65] CONSTANTINESCU V. Theory of turbulent lubrication(Turbulent lubrication theory based on mixing length

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2012ISAL0133/these.pdf © [M. Garcia], [2012], INSA de Lyon, tous droits réservés

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2012ISAL0133/these.pdf © [M. Garcia], [2012], INSA de Lyon, tous droits réservés 2. Hydrodynamic study and two-phase flow

Remerciements

appreciated the help and support that were given to me by my supervisors and colleagues