Numerical model of face seals for down-hole tools

(1)

Numerical model of face seals

for down-hole tools

Jeremy Wyart

Master thesis for the degree of Engineering Mechanics

Supervised by: Dr. Andrew Parry, Schlumberger Examiner: Pr. Henrik Alfredsson, KTH

(2)

Abstract

** Intentionally left blank **

Acknowledgements

In the present work, I present the work done for the master thesis I have done for my dou-ble degree at KTH together with the Ecole polytechnique. I realized it at Schlumberger in France.

I would like to thank my manager Andrew Parry for the work we have done on this project together. It has been a pleasure to do my thesis with him. I also woud like to thank Jeremy Cochain for the couple of months we worked together. I would like to thank the M&S team I was part of, Roel Van Os, Viet Tung NGuyen and Clement Texier for their warm welcoming and to have helped me out when needed. Finally, I would like to thank Pr Henrik Alfredsson whoaccepted to follow this thesis as an examiner.

(4)

1 Introduction

Schlumberger is the world’s leading supplier of technologies in the oil & gas industry. Present in 85 countries with 120 000 employees on the five continents. Schlumberger owns two business-orientated divisions :

- SLB is a supplier of products and services for exploration, drillings and completions of wells

- WesternGeco is the world’s largest seismic company and provides advance acqui-sition and data processing services.

Created in 1926 by the brothers Marcel and Conrad Schlumberger in Paris, the com-pany is based on their invention of resistivity measurement method developed around 1920. From 1929, resistivity measurements were done in every part of the world, for example Argentina, Ecuador, India, Japan, the Soviet Union, Venezuela or the USA. This technique historically helped a lot the growth of the oil business.

1.1 History

s

Figure 1: Photo of the spindletop well in 1901

The first oil well has been drilled by Colonel Edwin Drake in Pennsylvania in

1859. At this time, the well produced

4m3/day. The major discoveries occured

in the beginning of the 20th century,

with the first major discovery in 1899 in Spindletop (Texas) in the USA. The field

produced 16 000 m3/day. This period saw

the emergence of major oil Companies. They improved their development by in-tegrating the production, the refining and the transport of oil until commercializa-tion. In 1870, the Standard Oil is founded by Rockefeller, becoming a giant in the in-dustry that led to its dissolution in 1911 by the US supreme court. At some point, the company had the monopoly of oil in the USA. The company split in 33 smaller companies and most of them are today still in activity: Amoco, Exxon, Mobil, Chevron.

(5)

on core samples or cuttings. The first log showed that electrical measurements in a drill hole could help identify geological formations traversed by the drill. The Schlumberger brothers, Marcel and Conrad, stood behind this invention that they developed in France and tested at their family domain in Alsace, France. Later, they went abroad to sell and develop their method. At first, the big oil companies were not very interested. But after the success in Pechelbronn, France, where the Schlumberger brothers showed the efficiency of their method, companies were willing to try. Since then, the electrical log has been further developed, but the basic principle is still the same as in the original design of 1927.

t

Figure 2: Portraits of Conrad and Marcel Schlumberger

1.2 Oilfield activities and life of a well

The exploration and production of oil is now done using science and modern technologies. The first step is the oil prospection. When hydrocarbon traps are found, drilling and production can begin.

The construction of a well is executed in several steps. First the well is drilled

and logged. Then, casing is installed and cementing carried out. Since drilling, logging, casing, and cementing are performed step-by-step, they are repeated several times. After that, production tubing and packers are installed. Finally, perforations are done in the casing to permit the oil to flow into the tubing.

1.2.1 Oil prospection

To discover potential oil traps, large areas are first photographed by airplanes or satellites

and experts study the formation of rocks to evaluate the potentiality. This survey

(6)

Figure 3: The principle of an offshore seis-mic survey. The waves are reflected back by the different rock layers.

To continue the search and to be able to view oil and gas reservoirs that are buried under thousands of meters of sea or rock, seismic surveys are executed. They can be performed on land or at sea but the principles are the same: sound waves penetrate the many layers of rock. Each boundary reflects a part of the sound back to the surface. The rest continues down-ward. On the surface, special devices – geophones – pick up the reflected sounds. Depending on how long the reflection time is, the type of geological formation can be inferred.

Figure 4: The principle of an onshore seismic survey. The waves are reflected back by the different rock layers.

(7)

1.2.2 Drilling and Logging

Figure 5: Schematic principle of logging

while drilling. The drilling is called direc-tional drilling as it can be adapted to results from the logging.

The drilling is nowadays done with rotary drilling tools. The rotary drilling rig uses a rotary bit with rows of teeth that pen-etrate the rock The bit is attached to a drill pipe that consists of several joints of pipe. Joints are added to the drill pipe as the hole becomes deeper. Afterwards, the cuttings of rock must be moved out of the way. Otherwise, the drill bit would be hindered. To permit this, fluid circu-lates in the well. This fluid, called drilling mud, transports the cuttings to the sur-face where they are sorted out so that the drilling mud can be recycled in the

bore-hole. The fluid enters the well through the drill pipe and goes out through the drill bit. A huge pump on the surface moves the mud circulation system.

Figure 6: Schematic principle of casing. The deeper the hole is the smallest it becomes.

Logging is a set of techniques, which consists of col-lecting and registering in real time geological informa-tion about the inside of the earth. The logging per-mits, eventually, to confirm the surveys that have been done earlier. Sometimes they are right and there are hydrocarbons, sometimes there is none. Only one of seven exploration wells is developed into a production well. Much has occurred since the Schlumberger broth-ers succeeded with the first electric log in 1927. To-day there exists many different variants, each with its own specialty and field of application. The main types are: wireline logging and logging while drilling (LWD). Wireline logging is used in open and cased hole. LWD is performed, as the name indicates, while drilling.

1.2.3 Cementing and casing

To prevent the well bore from collapsing, it has to be cased. That means putting down a strong steel pipe and then cementing it.

(8)

deepens.

The casing provides a complete isolation of different zones, supports the casing, and protects the casing string. It is very important that fluids cannot migrate from one formation to another. This is to prevent, for example, oil leakage in nearby fresh-water reservoirs. Cement contains silica, alumina and iron oxide. Cement for wells sets through a chemical process that does not require air. This process consists of very complex chemistry.

1.2.4 Perforating

Upon completion of all the previous steps the well is almost ready for production. The casing and the cement are sealing the well from the oil and gas reservoirs. To let the oil enter the well, communication has to be established between the reservoir and the well. A special gun is lowered in the well and at the oil-producing zone where it shoots rather small holes in the casing. These holes are called perforations and do not only penetrate the casing but also the cement and a short distance into the rock.

As can be seen in fig. 7, many improvements have been brought to the perforating

during the 20th century. Today, it is possible to perform directional perforating i.e.

explosions that are not symmetrically distributed around the well.

Figure 7: The development of perforating.

1.3 Face seals and down-hole tools

Face seals are mechanical components that can be found on down-hole tools. They ensure the sealing between two environments: the inner fluid, usually motor oil and the outer fluid which is usually drilling mud as explained in 1.2.2.

(9)

It is important to notice that very different operating conditions are studied. Some tools have a life span of years when the face seal is running at constant speed. This is the case for example for ESP, electrical submersible pumps which stay at the bottom of wells to pump oil toward the surface. Leakage is here a critical parameter as the oil reserve is limited and must be sufficient for years. On the other hand, drilling tools experience numerous starts and stops as well as pressure pulses. In addition, in every cases, high pressures and high temperature are present.

Face seals are composed of a rotor and a stator (see figure 8) separated by a micro-scopic film fluid. This film ensures the cooling of the faces in contact and reduces the friction. The two faces are maintained close to each other by a spring. The spring can be attached on the rotor (as on figure 8) or on the stator. The face seal studied later on has a spring loaded rotor.

(10)

2 Presentation of the model

A face seal is modeled as two annular faces separated with a fluid of micro-metric

size. Our model is based on a specific seal from a specific tool but everything

ap-ply to all kind of seals present in down-hole tools. This work is based on a previous work done by R.F Salant in 2002. (see [Salant and Cao, 2005], [Cao and Salant, 2002], [Cao and Salant, 2003] for more details).

A table of all the parameters name used is recalled in appendix A

The model has been developed in Octave, an open-source software used for numerical computation. We used also several model in ANSYS mechanical and ANSYS fluent to model specific physical behaviour or to validate with a finite element analysis parts of our model.

2.1 Geometry

This section describes the face seal geometry given for this study (see figure 9). The hexagon shape has been modified to a circle for axisymmetric assumptions. The seal is inner pressured by a lubricating fluid. The outer diameter is in contact with the drilling mud. The rotor is made of silicon carbide and the stator is made of steel.

(11)

2.2 Multiphysics study

The complexity of the study of face seals is the coupling of several multi-scale physics interacting deeply with each other (see 10). Analysis of different phenomena are car-ried out and coupled with each other. Through numerical iterations, an equilibrium is reached. In addition, we have several scales involved, as the film fluid has a thickness of about 1 micron and the radius of the face seal is a few millimetres. This variation of scales from one parameter to another amplifies the effects on the study.

The analysis includes:

• A film fluid analysis with Reynolds equation. This analysis computes the pressure within the fluid, influencing the opening of the face seal as well as the leakage rate. • A contact mechanics analysis modelling the asperity-asperity interactions between the two faces. The repulsive force derived out of it is used to compute the opening force.

• A body dynamics analysis computing the faces position equilibrium.

• A thermal analysis computes the heat flux due to the viscous strain in the film fluid as well as the friction done by asperity contact.

• A deformation analysis derives the mechanical deformations of the faces due to pressure forces as well as the thermal deformations due to the heat flux at the interface.

(12)

(13)

3 2D Axisymmetric model

The model used is shown in figure 11. The fluid film has a local thickness h. h is

the sum of hm which is the minimal film thickness (the point where the two faces are

at the closest) and the deformation dh. Originally, the faces are taken as planar and deformations are computed taking the planar state as reference. The deformations are

referred to as dh. pi is the inner pressure at the inner radius ri. po is the outer pressure

at the outer radius ro.

Going from a steady model to a transient model has been done through the Duhamel’s method given in [Salant and Cao, 2005]. At each time step, we consider that all the forces are in equilibrium but for the heat transfer. This enables us to use the influence coefficient method as described in section 3.6 and accounting the thermal history on the heat flux computed at time step t.

Figure 11: 2D axisymmetric geometry used for the analysis

3.1 Force Balance

The floating seal is in equilibrium in the axial direction. The rotor is fixed and the stator experiences axial translation. Consequently, an axial force balance must be written on

the stator. We define two forces, the opening force Fopen and the closing force Fclose

that must be equal:

(14)

The opening force is simply written as the sum of the fluid pressure and the contact pressure due to asperities.

Fopen= 2π

Z ro

ri

prdr + Fcontact (2)

The expression of the closing force contains the hydraulic force of the face seal and has been studied in [?]. A face seal is said to be compensated when the inner pressure has an effect on the closing force. Figure 12 shows that a part of the closing pressure is due

to the inner pressure pi (red surfaces in figure 12). If we define Fs as the spring force,

FP I the closing force due to the inner pressure and FP O the closing force due to the

outer pressure, we can write

Fclose = Fs+ FP I + FP O = Fs+ pi π 4(d 2 B− d2I) + po π 4(d 2 O− d2B)

dB is the area-average bellows diameter defined such that d2B− d2I= d2O− d2B. Given

these expressions, we can write the closing force as

Fclose = Fs+ K · π 4 · (pi− po)(d 2 O− d2I) (3) with K = pi(d 2 B− d2I) + po(d2O− d2B) (pi− po) · (d2_O− d2_I) (4)

The validation of this procedure is given in section 5.4.

3.2 Fluid Mechanics Analysis

Using the hypothesis of lubrication (Reynolds equation) and the axisymmetric steady state, we can derive the pressure distribution from the following simplified Reynolds equation: ∂ ∂r(rh 3∂p ∂r) = 12µr ∂h ∂t (5)

Then the fluid is discretized using the control volume method, described in [Patankar, 1980] and reminded page 16. Following this, the Reynolds equation becomes

ap pi+ aw pi−1+ ae pi+1+ b = 0 (6)

The coefficient at time t + dt are calculated as follows

(15)

Figure 12: Force balance on the rotor. The opening force is composed of the pressure at the interface and the closing force is composed of the spring force and a force depending of the inner and outer pressures applied at the outer surfaces of the rotor. The inner

pressure applies on all the red surfaces (summed as FP I) and the outer pressure applies

on all the green surfaces (summed as FP O)

kW = ri−1h3i−1 ∆r kP = rih3i ∆r kE = ri+1h3i+1 ∆r (8) (9) To solve this equation, the only parameter needed is the film thickness distribution. h must be specified and coupled to the other analyses through an iterative scheme.

(16)

A note on the control volume method

The control volume method presented by Patankar considers a one dimensional diffusion equation d dx k ·dT dx + S = 0 (10)

Every grid point P has two neighbors W and E (see figure 13). To derive the dis-cretization equation, we employ grid points and values in figure 13. The integration of (10) integrated over the control volume gives

kdT dx e − kdT dx w + Z e w Sdx = 0 (11)

With linear interpolation , we evaluate the derivatives at points e and w. We obtain

ke(TE− TP)

δxe

−kw(TP − TW)

δxw

+ S∆x = 0 (12)

that gives equation (6) by rearranging the terms.

Finally, considering a constant conductivity for a given control volume with an equidistant grid, we write the heat flux at the interface, we obtain the interface conductivity as ke= 0.5 kE +0.5 kP −1

that gives expression (8).

(17)

3.3 Contact analysis

Once the film thickness distribution is determined, the contact pressure can be found. The model chosen for the contact pressure is a plastic asperity deformation model developed in [Lebeck, 1991]. The model deals with the contact of two rough surfaces by transforming it in a contact between a flat plane and an equivalent rough surface. The asperities that are higher than the distance between the two surfaces are truncated and considered in contact. The aspirities are taken as isotropic in the two spatial direc-tion and follows a Gaussian distribudirec-tion with a specified standard deviadirec-tion. The local contact pressure is given by:

Pc= H Z ∞ h(r) 1 σs √ 2πe (−z2 2σ2s)dz (13)

The parameters are the hardness of the material H, the equivalent roughness of the

surface σs and the local film thickness h.

The equivalent roughness is derived from the roughness of the two surfaces:

σs=

q σ2

1+ σ22

The contact force in the axial direction is consequently derived by integration of the contact pressure over the interface.

Fcontact= 2πH Z ro ri " Z ∞ h(r) 1 σs √ 2πe (−z2 2σ2_s)_dz # rdr (14)

3.4 Heat generation analysis

The performance of a mechanical face seal is affected by the fluid pressure and the contact pressure. But one of the dominant phenomena is the thermal deformation of the faces. In addition, the thermal history of the face seal must be taken into account as the diffusion of heat is not instantaneous. The thermal history is taken into account by the Duhamel’s method and explained below.

The frictional instantaneous heat generation is derived with the expression: q_{f ric}inst = 2πfsω

Z ro

ri

pcr2dr (15)

where fs is the friction coefficient and ω is the shaft speed.

The viscous heat generation follows the expression:

qinst_visc = 2πµω2

Z ro

ri

r3

(18)

where µ is the fluid viscosity dependent of the temperature.

Sallant in [Salant and Cao, 2005] expressed the temperature field as

T − T0 =

Z τ =t

τ =0

φ(r, z, τ, t − τ )dτ (17)

where φ(r, z, τ, t) is Duhamel’s auxiliary function, corresponding to the time response of the system to a heat source. φ is found from a numerical experiment using a finite element analysis and applying a constant heat flux at t = 0+. The resulting temperature at node i is found to obey the following empirical law

δTi− δToi

Tsi− Toi

= P φ(, t)

jTijqj(Qk() − Qk(0))

= f (t) (18)

where f(t) is the Duhamel’s function. The Duhamel’s function that Sallant used was

f (t) = (1 + b_t)−1. A better fitting of the response has been developed, taking the shape

f (t) =1 + C1 exp(− t τ1 ) + C2 exp(− t τ2 ) (19) with            C1 = 0.5735 C2 = 0.26547 τ1 = 0.03242 τ2 = 0.13707

The temperature influence coefficient are derived with an ANSYS finite element method analysis described in section 3.6.

3.5 Deformation Analysis

The surface deformation is separated in two causes: mechanical and thermal deforma-tion. The mechanical deformation are caused by the fluid pressure and the contact pressure while the thermal deformation is due to the heat generation. Both types of deformation are computed through the Influence Coefficient Analysis (see section 3.6). The influence coefficient method assumes that all the deformation are linear dependent of the load. The advantage is to reduce the computation cost and the accuracy remains very good as proven in section 5.3.

The mechanical deformation is computed from

dhm(i) =

X

j

IFij(p(j) + pc(j)) + IQijqj (20)

(19)

dhth(i) − dht=0th (i) = Z t τ =0 X j IQ_ij[qj(τ ) − qj(0)] ∂ ∂t[f (t)] (21)

Once the deformation is calculated, the film thickness becomes

h = hm+ dhm+ dhth (22)

3.6 Influence Coefficients

The influence coefficient method is based on the principle that the deformation and the temperature vary linearly with the load. Two kind of influence coefficients are calculated: the mechanical and the thermal.

For the mechanical influence coefficients, a 2D model of the rotor and stator is built where the interface is divided in n elements. For each of these elements a standard input pressure is made and the resulting deformation at nodes is printed. These elements

builds the n · (n + 1) Influence Coefficient Matrix IFij.

(20)

The two set of thermal coefficients, both for temperature computation and the ther-mal deformation computation, are computed with an ANSYS mechanical model. The first thing to include is the heat transfer coefficient at the surface of the elements of the seal. This coefficients are calculated using a global model and a local model in ANSYS Fluent. The global model enables us to determine the field of temperature as seen in figure 14.

Figure 15: Local model and Temperature Field of the face seal in ANSYS Fluent giving the heat transfer coefficients at the boundaries

The local model ( see figure 15), simplified to only keep the silicon carbide parts corresponding to the face seal itself. It is used to compute a set of heat transfer coef-ficients. The input is the temperature field taken from the global model and the heat transfer coefficients are written in a text file HeatCoeff global.prof. This text file is then transformed into an ANSYS mechanical readable file with a Fortran function.

3.7 Output parameters

The leakage rate is written as

Q = πρrh

3 6µ

dp

(21)

3.8 Computational procedure

Given all the models described in section 3.1 to section 3.7, a computational proce-dure is created to couple them. The computational sequence is described in figure 18. The figures 16 and 17 describes the way all the different routines and parameters work practically. Fluent model model3_nitrogen.dat.gz model3_nitrogen.caz.gz heat_transfer_coeff_abs_nitrogen.prof Fortran 90 h_proc.f90 heat_transfer_coeff_ansys_nitrogen.txt Ansys Mechanical thermal.txt Ansys Mechanical mechanical.txt Stator.dat Rotor.dat Stator.dat Rotor.dat Octave Main.m

(Include 14 routines : initiate_parameters.m, initial1.m, closeforce.m, adim_parameters.m, PressureBis.m, ContactForceTiers.m, HeatBis.m, DeformationBis.m, OpenForceBis,

tempera-ture.m, Performance.m, output.m, DeformOutput.m, PerformOutput.m)

Input..dat vt.dat performance4..dat hi.dat pi.dat ti.dat dthermal.dat dmechanical.dat set number of elements

set size of elements

set N nb of elt at interface set geometry set material parameters-set load

set N nb of elt at inter-face

set geometry set material parameters set load

Main.m set imax set input file names set hrelax

(22)

(23)

Initiate parameters : initiate_parameters.m Guess an initial h = hm +dh

Inititate influence coefficients : initial1.m

Adimensionalize the parameters : adim_parameters.m Calculate the closing force : closeforce.m Calculate the fluid pressure p : PressureBis.m Calculate the contact pressure pc : ContactForceTiers.m

Calculate the heat generation q : HeatBis.m Calculate the deformation dh : DeformationBis.m

Calculate Opening force : OpenForceBis.m

Calculate newh = hm+dh

If |new h - old h| < error

If |Cforce - Oforce|<error

Calculate new temperature new t : temperature.m

If |old t—new t| < error

t = λ*t+(1-λ)*old t

Set the new viscosity coefficient amu

Output the results in the result files.

h = hm + λ newh + (1—λ) oldh

Change the minimal film thicknes hm

No

h = new hm + dh

No

(24)

4 3D tilted model

Considering the 3D model, we created a 3D model with a tilt angle between the two faces. This tilt angle is called β. The azimuthal variation of film thickness implies some changes in equations, especially for computation of the pressure.

Figure 19: 3D model geometry. The rotor and the stator have a relative tilt angle that induces an azimuthal change of the fluid film thickness

4.1 3D pressure routine

The azimuthal variation implies some change in the pressure computation. First, the asperity have an impact on the fluid flow as described in [Patir and Cheng, 1979]. An average flow model for the asperity effect on the flow is taken. The Reynolds equation is now: 1 r ∂ ∂r ϕrh3 12µ ∂p ∂r +1 r ∂ ∂θ ϕh3 12µr ∂p ∂θ = ω 2 ∂h ∂θ (24)

where ϕ is defined as the pressure flow factor, taken as ϕ = 1 − Ce−ηhmeanh _{.The tables}

for C and η are given in [Patir and Cheng, 1979].

(25)

plausible. A work on cavitation understanding has consequently been done. Three dif-ferent examples from literature and finite element analysis have been compared to our routine forecasting cavitation.

This has been described by Payvar in [Payvar and Salant, 1992]. The Reynolds equa-tion given in equaequa-tion (24) is valid for non cavitating zones. In case of cavitaequa-tion, we have ∂ ∂θ ρ ρc h = 0 (25)

For this, a new form of the Reynolds equation is derived, using the cavitation index F and a single variable Φ for both regions with and without cavitation. The value of the cavitation index accounts for the presence or absence of cavitation at each node as described below.

In cavitating regions, the cavitation index F (r, θ) = 0 and we have the following equations:      p = pcav F = 0 ρ ρc = 1 + (1 − F )Φ (26) In non-cavitating zones, ( F Φ = p−pcav pin−pcav F = 1 (27)

Using both equations (26) and (27), we derive the Reynolds equation for cavitation: 1 r ∂ ∂r ϕrh3 12µ ∂(F Φ) ∂r +1 r ∂ ∂θ ϕh3 12µr ∂(F Φ) ∂θ = ω 2(pin− pcav) ∂ [(1 + (1 − F )Φ)H] ∂θ (28)

This equation can be non-dimensionalized if needed and gives: 1 ¯ r ∂ ∂ ¯r ϕ¯r¯h3∂(F Φ) ∂ ¯r +1 ¯ r ∂ ∂θ ϕ¯h3 ¯ r ∂(F Φ) ∂θ = γ ∂ [(1 + (1 − F )Φ)H] ∂θ (29) with γ = 6 µ ω r 2 i (pin− pcav)h2mean (30) The values of γ is useful to study qualitatively the probability of cavitation. Cavita-tion occurs for high γ.

(26)

Figure 20: 3D grid-points geometry for discretization of Reynolds equation apΦP + anΦN + asΦS+ aeΦE + awΦW = b (31) where                      an= rk_∆rn∆θFN as = rk_∆rs∆θFS ae = k_r∆θe∆rFE aw = k_r∆θw∆rFW ap = − rk_∆rn∆θFN +rk_∆rs∆θFS+k_r∆θe∆rFE+k_r∆θw∆rFW + γ (1 − FP)her∆r b = γ(he− hw)r∆r − (1 − FW)hwr∆rΦ (32) 4.2 Influence coefficients

The influence coefficients should be taken as symmetric. At first, the deformation given by ANSYS for an input load is not symmetric. In addition, very small variations are exponentially increased through numerous iteration, especially in case of small defor-mations (near axisymmetry state for example). Consequently, we forced the symmetry of the influence coefficients by computing only half of the influence coefficients with ANSYS. The complete set of influence coefficients is then created by symmetry of this halved set.

4.3 Deformation routine

(27)

computation of the stator deformation. An axisymmetric deformation is computed for the rotor when the stator deformation contains an azimuthal variation.

(28)

5 Validation procedure for 2D model

Some of the routines have to be verified against analytical or numerical solutions. As most of the local equations from the 2D model are still valid, some other physical phe-nomena have to be taken into account.

5.1 Pressure Validation

5.1.1 2D steady validation

For several input film thickness h, the analytical and numerical pressure has been plotted.

We took care to ensure several ratios of hmax

hmin to obtain different pressure profiles.

The integration of Reynolds equation (equation (5)) gives the following expression :

p = ro Z ri C1 rh3dr + C2 (33) h = constant

For a constant h, we obtain the following expression

p = C1 h3 0 ln(r ri ) + p(ri) (34) with C1 = p(ro) − p(ri) ln(ro/ri) h3₀ For a linearly varying h = a · r + b

p = p(ri) + C1· δ0ln(r) + δ1ln(ar + b) − δ2 a(ar + b)− δ3 2a(ar + b)2 (35) with                  C1 = (p(ro) − p(ri)) h

δ0ln(ro) + δ1ln(aro+ b) −_a(arδ_o2_+b)−_2a(arδ_o3_+b)2

i−1

δ0 = _b13

δ1 = −_ba3

δ2 = −_ba2

δ3 = −a_b

(29)

(a) Correlation between numerical and analyti-cal computation of pressure for a constant h

(b) Correlation between numerical and analyt-ical computation of pressure for a linearly in-creasing h

(c) Correlation between numerical and analyt-ical computation of pressure for a linearly in-creasing h

(30)

5.1.2 2D transient validation

In order to ensure that the sub-model is valid, a case for which an analytical method exists has been tested. The film thickness distribution is chosen as

h = 1 − t

c

The pressure routine solves the non-dimensionalised equation

∂ ∂ ¯r(¯r¯h 3∂ ¯p ∂ ¯r) = c¯r ∂¯h ∂t (36) where      ¯ r = _rr i ¯ p =_pp amb ¯ h = _hmeanh (37)

The boundary conditions are arbitrary chosen as (

¯

p = 7 at r = 1 ¯

p = 1 at r = ro/ri

For this test case, the analytical solution is

¯ p =C2+ C1 ¯ h3 ln(r) − c¯¯r2 2τ ¯h3 (38) with ( C1 = ¯ h ln(ro/ri) h po− pi+_{4τ ¯}c_h3( ro ri − 1) i C2 = pin+_{4τ ¯}c_h3

The solution given in figure 22 shows a perfect fitting between analytical results and numerical results.

(31)

Figure 22: Analytical verification of transient pressure calculation

5.2 Validation of the model with 2D steady axisymmetric solution

A sample run with constant pressure difference and constant shaft speed has been carried out. The parameters are the same as in the 2D axisymmetric run. We can see that the 2D transient model is consistent with the results given by the 2D steady model.

5.3 Validation of the influence coefficient method

(32)

(a) Stator mechanical deformation (b) Stator thermal deformation

(c) Rotor mechanical deformation (d) Rotor thermal deformation

(e) Temperature at the interface

(33)

(e) Temperature at the interface _{(f) Heat flux and pressure at the interface}

(34)

(e) Temperature at the interface _{(f) Heat flux and pressure at the interface}

(35)

5.4 Validation of the Closing Force

The closing force formula in equations (3) and (4) has also been verified with a finite

element analysis. This is done by carrying out an ANSYS finite element model to

determine the reaction force on the spring due to the application of pi = 400 kP a and

po = 100 kP a. This case is without spring pre-tension so the spring force does not

come into play. Using the formula derived for closing force we get 89.6N . This value corresponds to the value found in ANSYS below (figure 26).

Note that the closing force equation does not take into account the change in dB due

to deformation of the bellows in the radial direction caused by the pressure.

The spring force is a function of axial pre-setting of the spring as shown in the graph below in figure 27. The axial pre-setting is the initial displacement of the spring compared to a reference.

(36)

(37)

6 Validation procedure for 3D model

The complexity of such a code requires some validation as small changes in some pa-rameters can change drastically the resulting behavior.

6.1 Validation of the fluid pressure routine

The complexity of cavitation happening in a face seal imposed us to carefully validate our pressure computation. We developped a solver that we validated with other solvers as well as with litterature results. Three cases have been studied: the work done by Payvar who created a solver, Lebeck’s results and a ANSYS Fluent finite element analysis. Payvar’s work on cavitation

Two articles dealing with cavitation have been studied. A work on reproducibility has been first done usign the code given by Salant and described in [Cao and Salant, 2003]. This code is based on the work of ([Payvar and Salant, 1992]) expanded and applied for a different configuration by [Xiong and Wang, 2012].

Reproducibility We first reproducedthe results in the case given in [Payvar and Salant, 1992].

The results presented in the publication and the ones we reproduced are given in figure 28. The parameters are given in table 1.

Non-dimensionalized film thickness _h h

mean = 1 + (B − 1)

η−1

η0−1+ A cos

2₍N θ 2 ).

We were unable to reproduce the results in [Payvar and Salant, 1992]. However, as shown later on, we have consistent solution between Payvar’s code and the chosen pressure code. (see figure 28)

Parameters Value

Inner radius ri 0.110 m

Outer radius ro 0.123 m

Sealed pressure ps 5.2697 MPa

Ambient presure pamb 0.4459 MPa

Cavitation pressure pcav 0.2047 MPa

B 0.5

A 0.625

η0= r_ro_i 1.116

Source parameter γ variable

Table 1: Parameters for the test case given in [Payvar and Salant, 1992]

(38)

(a) γ = 51 i i i i i i i i n i i n i u i i i i i i i j m i i i i i i u i i i M i i L i i i i i i i H M i n U ! - . ! ! l l l i : i ) ! U l \ l l ] l ! l l U l U l U l l l l l l l l l l U ! l l l i m i l l l l n i m i i i ' . m i i m i i m m i i i m i i m i u m i i m i m i m i i i n ; i u i LI : i ; s 11111111: u 11 LiiA 11 n 111 M 111 l i n 1111111 u 11

7 = 8 . 5

7 = 5 1

7 = 8 5

! l U ! J ! ! U U ! ! i ! ! H J l X l ! i n i U J . U U l j m J U m i J U l l U l l l l i i u i i i i m i i i j u m u i i i i i i u J j a i J t n u i i i u i i m i M i i m ; m ! ! L ! i ! ! i m i u u u u J i i L i i m u u r L j i i i u i u i i m i m i

Fig. 7 Development of the cavitation zone with increasing values of the parameter y. A = 0.625, B = 0.5, ij0 = 1.116.

the parameter A was allowed to vary to study the effect of the

wave amplitude. Thus the two main parameters studied were

7 and A. A 30 X 30 grid was used for all the calculations.

Figure 3 shows the pressure distribution along the

circum-ferential line t\ = 1.04 as a function of 6 when A = 0.625 and

7 was allowed to vary from 8.5 to 171. For low values of 7,

i.e., 7 < 35, no cavitation zone was observed. The curves for

7 = 51 and 7 =171 show the increasing extent of the cavitation

zone. Figure 4 shows the variation of pressure along r\ = 1.04

as a function of wave amplitude/! at 7 = 171. For low values

of A, i.e., A < 0.0625, no cavitation was observed along the

line 17 = 1.04 but at A = .0625 cavitation is imminent. The

A =0.375

' , : . ' : ; ; H I i n i i n i t i i i i i i n i i i u i i n i i i i i m i i n i i M r M r r ! ? : : ! : : i ! i i i i i i H i i m i ! i i i i i i i i i i i i i i H i i i i m i i ! ! U i i m < ) 1 1 : ; : 111! i m n 11111: i i u 11111111111111 n i m 11 > i>! 111111 > : : ! ! : ' . i u i m ] ' , u i i ! i ! n i u m i i i ] i m i n i : m i i ! ! i n n ; H i : : i s u n n m i 11 • i i i m m i i i i i m m m i i u m m i m i ! ! : : m m m i i i i ! i i L i i i m i n i u i i u m i l 111111 m m m i : ; ! i : i ! i m i i 1111 i U L i : < ' . i i i i m i i 111 m i i i i i i i m i i i i i u i ! !

A =0.625

Fig. 8 Development of the cavitation zone with increasing wave am-plitude parameter A. y = 171, B = 0.5.

curves for A = . 125 and A = 0.625 show the rising maximum

pressure and the increasing extent of the cavitation zone.

Fig-ures 5 and 6 show similar plots of the quantity — along the

Pc

line 7) = 1.04. They indicate that the partial film content is

continuous at film rupture but has a discontinuity at film

reformation as indicated by the reformation condition in Eq.

(7). Figures 7 and 8 show a pictorial view of the development

and enlargement of the cavitation zone with increasing 7 and

A, respectively. These figures show two neighboring periodic

segments of the seal face with the left edge representing the

radial line of maximum film thickness. Figure 9 shows the

(b) γ = 51 (c) γ = 85 i i i i i i i i n i i n i u i i i i i i i j m i i i i i i u i i i M i i L i i i i i i i H M i n U ! - . ! ! l l l i : i ) ! U l \ l l ] l ! l l U l U l U l l l l l l l l l l U ! l l l i m i l l l l n i m i i i ' . m i i m i i m m i i i m i i m i u m i i m i m i m i i i n ; i u i LI : i ; s 11111111: u 11 LiiA 11 n 111 M 111 l i n 1111111 u 11

7 = 8 . 5

7 = 5 1

7 = 8 5

the parameter A was allowed to vary to study the effect of the

wave amplitude. Thus the two main parameters studied were

7 and A. A 30 X 30 grid was used for all the calculations.

Figure 3 shows the pressure distribution along the

circum-ferential line t\ = 1.04 as a function of 6 when A = 0.625 and

7 was allowed to vary from 8.5 to 171. For low values of 7,

i.e., 7 < 35, no cavitation zone was observed. The curves for

7 = 51 and 7 =171 show the increasing extent of the cavitation

zone. Figure 4 shows the variation of pressure along r\ = 1.04

as a function of wave amplitude/! at 7 = 171. For low values

of A, i.e., A < 0.0625, no cavitation was observed along the

line 17 = 1.04 but at A = .0625 cavitation is imminent. The

Journal of Tribology

A =0.375

A =0.625

curves for A = . 125 and A = 0.625 show the rising maximum

pressure and the increasing extent of the cavitation zone.

Fig-ures 5 and 6 show similar plots of the quantity — along the

Pc

line 7) = 1.04. They indicate that the partial film content is

continuous at film rupture but has a discontinuity at film

reformation as indicated by the reformation condition in Eq.

(7). Figures 7 and 8 show a pictorial view of the development

and enlargement of the cavitation zone with increasing 7 and

A, respectively. These figures show two neighboring periodic

segments of the seal face with the left edge representing the

radial line of maximum film thickness. Figure 9 shows the

JANUARY 1992, Vol. 114 / 203

(d) γ = 85 (e) γ = 428 i i i i i i i i n i i n i u i i i i i i i j m i i i i i i u i i i M i i L i i i i i i i H M i n U ! - . ! ! l l l i : i ) ! U l \ l l ] l ! l l U l U l U l l l l l l l l l l U ! l l l i m i l l l l n i m i i i ' . m i i m i i m m i i i m i i m i u m i i m i m i m i i i n ; i u i LI : i ; s 11111111: u 11 LiiA 11 n 111 M 111 l i n 1111111 u 11

7 = 8 . 5

7 = 5 1

7 = 8 5

the parameter A was allowed to vary to study the effect of the

wave amplitude. Thus the two main parameters studied were

7 and A. A 30 X 30 grid was used for all the calculations.

Figure 3 shows the pressure distribution along the

circum-ferential line t\ = 1.04 as a function of 6 when A = 0.625 and

7 was allowed to vary from 8.5 to 171. For low values of 7,

i.e., 7 < 35, no cavitation zone was observed. The curves for

7 = 51 and 7 =171 show the increasing extent of the cavitation

zone. Figure 4 shows the variation of pressure along r\ = 1.04

as a function of wave amplitude/! at 7 = 171. For low values

of A, i.e., A < 0.0625, no cavitation was observed along the

line 17 = 1.04 but at A = .0625 cavitation is imminent. The

Journal of Tribology

A =0.375

A =0.625

curves for A = . 125 and A = 0.625 show the rising maximum

pressure and the increasing extent of the cavitation zone.

Fig-ures 5 and 6 show similar plots of the quantity — along the

Pc

line 7) = 1.04. They indicate that the partial film content is

continuous at film rupture but has a discontinuity at film

reformation as indicated by the reformation condition in Eq.

(7). Figures 7 and 8 show a pictorial view of the development

and enlargement of the cavitation zone with increasing 7 and

A, respectively. These figures show two neighboring periodic

segments of the seal face with the left edge representing the

radial line of maximum film thickness. Figure 9 shows the

JANUARY 1992, Vol. 114 / 203

(f) γ = 428

Figure 28: Comparison between the reproduced results (left) and Payvar’s results (right) with parameters given in table 1. The results show no cavitation in case γ = 51 or γ = 85

(39)

(a) γ = 51 (b) γ = 51

(c) γ = 85 (d) γ = 85

(e) γ = 428 (f) γ = 428

(40)

We also verified that the pressure calculated without cavitation was consistent with a direct solver of the pressure field. The two fields of pressure were equivalent (including negative pressure).

Validation with Lebeck In [Lebeck, 1991], an example of a face seal pressure field

calculation is given. Our code reproduced the results with an good accuracy. The results are given in figure 30.

The film thickness is taken as h = h0+ ha∗ (1 − cos(nθ))

Parameters Value

Cavitation pressure pcav 0 MPa

Viscosity µ 680e-6 Pa/s

Asperity height hmean 0.1e-6 m

minimal film thicness h0 0.284e-6 m

film thickness variation ha 0.1e-6 m

rotation speed ω 188.5 rad/s

Table 2: Parameters used to define Lebeck’s case

(a) pressure field computed by our code (b) pressure field computed by Lebeck

(41)

Fluent comparison A Fluent model has been produced with the following parameters

in table 3. The geometry of the face seal is taken as: h = h0+ hacos(nθ/2)2. In Fluent,

no boundary conditions on pressure can be set up. Consequently, an inner radial velocity of 0.25m/s has been taken. The results are given in 31.

Parameters Value

Sealed pressure ps variable

Ambient presure pamb variable

Viscosity µ 0.02 Pa/s

Asperity height hmean 7e-8 m

minimal film thicness h0 5e-6 m

film thickness variation ha 5e-6 m

rotation speed ω 188.5 rad/s

normal velocity of the fluid ur 0.25 m/s

Number of peaks of h n 3

(42)

(a) pressure field computed by our code

(b) pressure field comparison at outer radius, inner radius and at the centered radius

Figure 31: Comparison between the pressure field computed by Fluent at ri, ro and

(ri+ ro)/2 and our code with parameters given in table 3. The results give the same

(43)

6.2 Validation with 2D axisymmetric code

As to validate the solution, a comparison between the 3D code in axisymmetric geometry and 2D solution from the steady 2D code has been done. The result is given in figure 32. The parameters used in this case are given in table 4

Parameters Value

Viscosity µ 0.003 Pa/s

Asperity height hmean 7e-8 m

rotation speed ω 377 rad/s

Hardness H 1e9 Pa

Friction coefficient fs 0.1

Spring Force Fspring 80 N

Table 4: parameters used to define the axisymmetric case

(a) film thickness (b) pressure

(44)

7 Results for a sample case for 2D transient conditions

Two kinds of results have been obtained for the same initial parameters and initial state. The first one is done varying linearly the shaft speed. The second one is done varying the pressure difference between the inner pressure and the outer pressure.

7.1 Results with regard to shaft speed

This run models a serie of start and stop of the shaft. The time conditions are shown in figure 33. The results versus time is shown in figure 34 and the plot of the parameters through the interface at several moments in one period are shown in figure 35.

One can notice that for a transient analysis, the film thickness and the leakage are smaller than in the steady analysis.

(45)

(a) Fluid film thickness (b) Leakage rate

(c) Heat Flux (d) Temperature

(e) Maximum Contact Pressure

(46)

(a) Fluid film thickness (b) Temperature

(c) Pressure (d) Contact Pressure

(47)

7.2 Results with regard to pressure pulses

This section represents the pressure pulse given in the bottom hole. This pressure pulses load the face seal and may be one of the cause of the failure of face seal. The response to such pressure pulse has been plotted in figures 36 for the load, figure 37 for the time response and in figure 38 for the plot of the interface at several time steps in the periodic cycle.

The results shows that the code is perfectly suitable to study the pressure pulse at the conditions to redefine correctly the time steps and the frequency of the pulses (normally 12Hz).

(48)

(a) Fluid film thickness (b) Leakage rate

(c) Heat Flux (d) Temperature

(e) Maximum Contact Pressure

(49)

(c) Pressure (d) Contact Pressure

(50)

8 Results for 3D model for several cases

Three test cases has been computed with different tilt angle β from β = 1 · 10−7 to

β = 1 · 10−5. The results are summarized in the table 5, and the figures for β = 1 · 10−7

are plotted in figure 39 and 40. The other are put in appendices.

β = 1 · 10−7 β = 1 · 10−6 β = 1 · 10−5

Ambient Pressure Pamb 1.0132 hPa 1.0132 hPa 1.0132 hPa

Sealed Pressure Ps 7.09582 hPa 7.09582 hPa 7.09582 hPa

Maximum Pressure P 88 bar 83 bar 53 bar

Minimum Pressure P 0.9 bar 0.9 bar 0.9 bar

Cavitation Pressure Pcav 0.9 bar 0.9 bar 0.9 bar

Max Film Thickness h 0.64 µm 0.68 µm 1.1 µm

Min Film Thickness h 0.36 µm 0.37 µm 0.52 µm

Max Deformation dh 0.28 µm 0.26 µm 0.19 µm

Min Deformation dh 0 µm 0 µm 0 µm

Max Contact Pressure Pc 2.8 bar 1.8 bar 0.0002 bar

Min Contact Pressure Pc 0 bar 0 bar 0 bar

Max Temperature T N/AoC N/AoC N/A oC

Min Temperature T N/AoC N/AoC N/A oC

Pressure opening force fvisc N/A N N/A N N/A N

Contact opening force fcont N/A N N/A N N/A N

Table 5: Results for the three chosen cases

(51)

(c) Fluid Pressure (d) Contact Pressure

(e) Heat flux (f) Radial Velocity

(52)

(a) Fluid film thickness (um) (b) Initial Geometry (um)

(c) Fluid Pressure (bars) (d) Contact Pressure in bars

(e) Heat flux (f) Radial Velocity in m/s

(53)

9 Conclusions

(54)

References

[Ayadi et al., 2014] Ayadi, K., Tournerie, B., and Brunetiere, N. (2014).

Optimisa-tion de la Modélisation des Régimes de Fonctionnement des Garnitures Mécaniques

d’Etanchéité : Analyse Théorique et Expérimentale.

[Bruneti`ere, 2010a] Bruneti`ere, N. (2010a). An analytical approach of the

thermoelas-tohydrodynamic behaviour of mechanical face seals operating in mixed lubrication. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineer-ing Tribology, 224:1221–1233.

[Brunetière, 2010b] Brunetière, N. (2010b). Les garnitures mécaniques Étude théorique

et exp´erimentale. Habilitation `a diriger des recherches, Universite de Poitiers.

[Bruneti`ere and Modolo, 2009] Bruneti`ere, N. and Modolo, B. (2009). Heat transfer in

a mechanical face seal. International Journal of Thermal Sciences, 48(4):781–794. [Cao and Salant, 2002] Cao, B. and Salant, R. F. (2002). Analytical model of an mwd

tool mechanical seal. (Schlumberger internal report).

[Cao and Salant, 2003] Cao, B. and Salant, R. F. (2003). Analytical model of an mwd tool mechanical seal. (Schlumberger internal report).

[Djama¨ı et al., 2010] Djama¨ı, A., Bruneti`ere, N., and Tournerie, B. (2010). Numerical

Modeling of Thermohydrodynamic Mechanical Face Seals. Tribology Transactions, 53(3):414–425.

[Incropera et al., 2007] Incropera, F. P., DeWitt, D. P., Bergman, T. L., and Lavine, A. S. (2007). Fundamentals of Heat and Mass Transfer, volume 6th of Dekker Me-chanical Engineering. John Wiley & Sons.

[Lebeck, 1991] Lebeck, A. O. (1991). Principles and Design of Mechanical Face Seals. Wiley-Interscience, 1st edition.

[Lienhard IV and Lienhard V, 2012] Lienhard IV, J. H. and Lienhard V, J. H. (2012). A Heat Transfer Textbook, volume 82. Phlogiston Press, fourth edition.

[Patankar, 1980] Patankar, S. (1980). Numerical heat transfer and fluid flow. McGraw-Hill Book Company.

[Patir and Cheng, 1979] Patir, N. and Cheng, H. S. (1979). Application of Average Flow Model to Lubrication Between Rough Sliding Surfaces.pdf. Transactions of the ASME, 101:220–230.

[Payvar and Salant, 1992] Payvar, P. and Salant, R. F. (1992). A Computational

(55)

[Salant and Cao, 2005] Salant, R. F. and Cao, B. (2005). Unsteady Analysis of a Me-chanical Seal Using Duhamel’s Method. Journal of Tribology, 127:623–631.

(56)

Appendices

A Table of variables 57

B Discretization of 3D Reynolds equation with cavitation 58

(57)

A

Table of variables

Octave Physical

Definition

name variable Unit

number of nodes imax ∅

shaft speed an an [rpm]

shaft speed omga ω [rad/s]

inner diameter ri ri [m]

outer diameter ro ro [m]

non-dimensionalized radius r

ambient pressure pamb pamb [Pa]

non-dimensionalized inner pressure pin pin

pamb

non-dimensionalized outer pressure po po

pamb

non-dimensionalized fluid pressure p _pp

amb

non-dimensionalized contact pressure pc pc

pamb

reference temperature tref Tref

non-dimensionalized temperature t _TT

ref

heat flux q q [W]

mean surface roughness hmean hmean [m]

non-dimensionalized film thickness h _h h

mean

non-dimensionalized minimum film thickness hm hm

hmean

non-dimensionalized spring force fspring Fspring

pamb·r2i

non-dimensionalized closing force Cforce Fclose

pamb·r2i

non-dimensionalized opening force Oforce Fopen

pamb·r2i

viscosity amu µ [kg m−1s−1]

leakage gleakm Q [m3/h]

(58)

B

Discretization of 3D Reynolds equation with cavitation

This describes in detail the discretization of the Reynolds equation including cavitation used in section 4.1. Z [29]rdrdθ (39) ⇔ ϕrh3∂(F Φ) ∂r dθ n − ϕrh3∂(F Φ) ∂r dθ s + ϕh 3 r ∂(F Φ) ∂θ dr e − ϕh 3 r ∂(F Φ) ∂θ dr w (40) = c ([(1 + (1 − F )ϕ)hrdr]_e− [(1 + (1 − F )ϕ)hrdr]_w) ⇔rkn∆θ ∆r FNΦN− rkn∆θ ∆r FPΦP − rks∆θ ∆r FPΦP + rks∆θ ∆r FSΦS + ke∆r r∆θ FEΦE − ke∆r r∆θ FPΦP − kw∆r r∆θ FPΦP + kw∆r r∆θ FWΦW = c (he− hw) r∆r + c ([(1 − F )Φh]e− [(1 − F )Φh]w) r∆r (41) with k = ϕh3 (42) ⇔rkn∆θ ∆r FNΦN+ rks∆θ ∆r FSΦS+ ke∆r r∆θ FEΦE+ kw∆r r∆θ FWΦW − rkn∆θ ∆r FN + rks∆θ ∆r FS+ ke∆r r∆θ FE + kw∆r r∆θ FW ΦP = c(he− hw)r∆r + c ((1 − FP)ΦPhe− (1 − FW)ΦWhw) r∆r (43) ⇔rkn∆θ ∆r FNΦN+ rks∆θ ∆r FSΦS+ ke∆r r∆θ FEΦE+ kw∆r r∆θ FWΦW − rkn∆θ ∆r FN + rks∆θ ∆r FS+ ke∆r r∆θ FE + kw∆r r∆θ FW ΦP = c(he− hw)r∆r + γ ((1 − FP)ΦPhe− (1 − FW)ΦWhw) r∆r (44)

The final discretization equation is :

rkn∆θ ∆r FNΦN+ rks∆θ ∆r FSΦS+ ke∆r r∆θ FEΦE+ kw∆r r∆θ FWΦW − rkn∆θ ∆r FN+ rks∆θ ∆r FS+ ke∆r r∆θ FE+ kw∆r r∆θ FW + γ (1 − FP)her∆r ΦP = c(he− hw)r∆r − γ(1 − FW)hwrΦW∆r (45)

(59)

(60)

C

Results for 3D case

(61)

(a) Fluid film thickness (b) Initial Geoemtry

(62)

(63)

(a) Fluid film thickness (b) Initial Geoemtry

Numerical model of face seals for down-hole tools