A new hydrodynamic lubrication theory for bilinear rheological fluids

I4lzL

ELSEVIER

I. Non-Ncwton~an

Fluid

Mech..

56 (19951 25.1 266

A new hydrodynamic lubrication theory for bilinear rheological fluids C . W . Wu

Received

‘.*,

H.X.

25 June 1993; in revisrd

Sun

form

b

14 July 1994

Abstract Hydrodynamic lubrication of a bilinear rheological fluid plays an important role in the so-called electrorheological (ER) “smart’ bearings. There may not. however, be a reliable and efficient algorithm for the non-Newtonian lubrication problem of bilinear rheological fluids. In this paper, a finite element trethod (FEM), together with a mathematical programming solution technique, are proposed for solving such a non-Newtonian lubrication

problem. lubricated Tichy’s should

This method was used to study the hydrodynamic lubrication by a bilinear rheological lubricant. The computed results

of a journal

were compared work for the Bingham rhcologieal model. It was shown that the Bingham be regarded as a limiting case of the bilinear rhcological model.

Ke~or&

Bilinear

rheological

model;

Bingham

model:

bearing with model

Lubrication

1. Inrroduction

The bilinear rhenlogical model, which was also called the biviscous model [I]. is often used to study the electrorheological (ER) “smart” hearings. Stangroom [2] and Dimarogonas and Kollias [ 31 described the interesting performances of the ER “smart” bearings. Tichy [4] proposed a differential method with an iterative

l Corresponding author. Present address: Advanced Materials Technology Materials Science and Engineering. North Carolina State University. Raleigh,

0377-0257/95/W9.S0

SSDI

:r> 1995 - Elsevier

0377-0257(94)01277-6

Science

B.V.

All

rights

reserved

Laboratory, Dept. of NC 27695-7517. USA.

254

C. W. Wu. HA’. Sun / J. Non-Newtonian

Fig. I. Bilinear

Flurd Mech.

rheological

56 (1995) 253. 266

model

algorithm to study Bingham plastic flow in a journal bearing. but he studied only the Sommerfeld type boundary condition which is most easily dealt with in lubrication mechanics. Thi:, paper will give a reliable and generalized finite element method (FEM) for the lubrication mechanics of bilinear rheological fluids. Consider the bilinear rheological model shown in Fig. I. The rheological relationship can lx written as ’ = { S’tign,,) + q[j’ 0 when the viscosity ratio r cal model for small r,, and model becomes Newtonian toward infinity.

2. FEM and mathematical 2.1.

Cottstitutive

control

(1) - sign(3)r 0/q er] II1: 1”:” = qc/q >> I, the model approaches the Bingham rheologiapproaches the limiting shear model [ 51 for large T,,. The with viscosity 1 if q, = 0 or with viscosity qr if q, tends

programming equations

for

bilinear

rheological

model

For completeitess, some assumptions often used in lubrication mechanics are repeated as follows. (I) The pressure is constant across the film thickness. (2) The viscosity and density of the lubricant are constant. (3) The lubricant is in laminar flow. (4) There is no slip at the interfaces of the lubricant and lubricated surfaces. (5) The viscous heating is negligible. In the case of twodimensional flow, the bilinear rheological relationship can be written (see Fig. 2) as T = tJ.(f - A(” + A’y,

(2)

C. W. WI. HA’. SW / J. Non-Newonion

Fig. 2. Geometry

FM

Mech.

of the control

Fig. A”’

= AC

2, it can

be shown

variables

which

satisfy

the

that

-AD

= (T - T”)/tl -(t -(T

255

vanablc~

where L(I) and I”’ are called the constitutive control following relations: (a) I(‘)> 0 if T > T,, and I”’ = 0 if T < To, (b) i,(” > 0 if 7 < -TV and A”] = 0 if T L - 5”. (c) AlI’ il2, = 0, 20) 2 0, I’?’ 2 0. From

56 (1995) 253. 266

- f”)/k

-f")y.

(3)

i’2’ = -(r + *,,)$$

(4)

Similarly

where 5 = v,/r]. Combining Eqs.

(2) -(4)

gives

the

following

yield

functims:

2%

C. W. Wu. HA’. Sun I J. Non-Newtonian

Fig. 3. Lubrication

The yield relationships: jlil

functions

~0’ = 0,

and

control

variables

f = au/t$

Substituting gives

p

54 (1995) 253-266

film geometry.

< 0 if j, “‘~0 and/“)=0 that the lubricated surfaces slide in the x direction can be obtained [6]:

Supposing fluid velocity of Cameron

Using

.I”’

Fluid Mech.

i”’

(i = 1,2)

have

the

following

ifi”‘> in the x direction (see Fig. 3), the from Eq. (2) referring to the work

yields

‘79. (7) into

(5a)

and

(5b).

and

I v ap u,-u, -_-__ 1-5 A""+z.---+__ cqa.r h

in:roducing

slack

variables

P”’ and

c(‘).

(1) _ j,(z)) dy - 2 + 0”) = 0, @a)

1 . '2'---Y ---,.

ap-----+--L+ u, eqax h

1-c

hap 2:?j ax

+-I h(i(I) _ i(2))dJ _ Ir! + p’21 h sr) 3l =

0,

(W A’0

“‘I)

= 0 ,

rl”‘20.

u”’

2 0

(i = 1.2).

( 8~)

C. W. Wu. H.X. Sun / J. Non-Newtonian

Fluid Mech. 56 (1995) 253.-?ti

251

The Eqs. (8a) -(8c) are caLledthe constitutive control equations of the lubrication problem of bilmear rheologicdi fluids. 2.2.

Reynokdr

r.vpe equation

awd

its oariation

ftfrcncriot:al

Substituting Eq. (6) into the following volume flow continuity equation: u(x,,:v)dy + ; = 0,

(9)

yields a Reynolds type equation.

(js")

-

/I h if2’ dy dy - -2 (>.“I - I.“‘) 4,s sIi

1

Its variational functional [7-91 is

( IO)

where the control variables A”’ and I.(” will be determined by the constitutive control equations (8a) -( 8~). Because fluids can not resist large tensile stresses,it is generally assumed, in lubrication mechanics, that a negative film pressure is impossible. In other words, the variational functional ( 1I) should have p ~0 as an additional constraint besides the constraint equations (8a) -(8c). Consequently, the problem mentioned above can be summarized as the following optimization problem [8-IO] with complementary constraints. Find film pressurep min

J(P),

(124

s.t. p tp, ap/a.r, 5, i+. X, y) + ~(1)= 0, j.lG

p 20.

~(1)= 0.

2.(')2 0,

di) -> 07

(W (i = 1.2).

(12c) (i2d)

Eqs. (12a)-( 12d) can be used to deal conveniently with the Reynolds pressure boundary condition. The Reynolds boundary condition will be automatically satisfied as long as the selectedanalysiszone is large enough to include the positive pressure zone. In other words, the size of the analyzed zone will not affect the computed results of the pressuredistribution as long as it is large enough to cover the whole positive pressurezone. Too large an analyzed zone will, however, lead to much computing work. The mathematical proof for a similar problem was given in Ref. 1) I].

258

C. W. Wu, H.X.

Sun 1 J. Non-Newlonim

2.3. Finite element discretization

Fluid Mech.

and mathematical

56 (1995) 253-266

programming

Since the tilm pressure is assumed to be constant in the y direction, the pressure elements &e linear. However, L”‘(.w,y) is still a function of both x and ), so plane elements have to be used to determine the yielding state. In order to distinguish between them, the latter are called sub-elements or constitutive control elements. It is assumed that a pressure element contains r sub-elements and each of the sub-elementshas only one yielding state, i.e., A(” is constant in each sub-element. If the lubricated zone is divided into m pressure elements, there are m x r sub-elements in the system. In the present analysis,the following dimensionlessparameters are used: x = x/B,

y=

Y/kin

9

H = h/h,,,

T = tU,/B,

I\“’ = ;.“I/, mm/u,, P =phL/(~BVrk v: = u,/u,, yCi) = #)hmi,/V,.

Q = a/( CL). u: = UJU,,

vr = VI Ix=” + U*lx=fh

u=vf+v:,

fi = aH/aT,

After finite element discretization and the introduction of a Lagrangian multiplier vector, L, Eqs. (12a)-( 1213)can be rewritten in the following discretization form: min J(P) = f PKP - P’ [BU - WH - Q + DA] - L*P,

(134

s.t. CP + MA + F + V = 0,

(W

Fig. 4. Schematic

of a journal

bearing.

C. W. Wu. H.X. Sun 1 J. Non-New~onion

Fluid Mech.

56 (1995) 2%

266

259

9U e

I ^^

=0.

ILP”

a

IH

I *10

E OH

\

6

0 0

n/2

2n/3

W/6

n

0 Fig. 5. Pressure profiles for diRerent values of 5 under the half-Sommerfeld rz = IO;-.---, r: = 0: o. Tichy’s solu;ion [4].

boundary

condition.

-- -.

PA=O,

v 2 0,

A 20,

(13c)

LTP=O,

L >c,

P r 0,

(134

where the related matrices and vectors are defined in the Appendix. Letting 6J(P)/6P = 0 in Eqs. ( I3a) -( 13d) yields the following linear complementary problem: (14a) VA=O,

YZO,

A 2 0,

(14’4

LTP=O, L 20. P20. (144 where Q’ = BU - Q - WI& Eq. ( 14) can be solved using many methods. Here, the Le;ake method [IO] was used.

3. Numerical

analysis

for

a journal

bearing

A journal bearing rotating with velocity w is shown in Fig. 4. The dimensionless film thickness [4] H =( I +E cos O)/(1 -E), where 0 =x/R, E = e/c and

c=R,,-R,<
260

C. W. Wu. H.X.

40

Sun / 1. Non-Newronian

Fluid Mech.

56 (1995) 253-266

e =O. 8

30

p”

Half

somerfeld

20

10

0 2n/3

W/6

R

7n/6

49 Fig. 6. Pressure profiles

for dilTerent

values of 5 under differen*

boundary

conditions

s 9 Fig. 7. Core formation and yielding zones under diRerent boundary conditions (t a 100. r; = 10. E aO.8). B Sub-elements with core formation (no yielding) under the half-Sommerfeld boundary condition; H sub-elements with core formation (no yielding) under the Reynolds boundary condition; -, boundaris of yielding and no yielding zows under the Sommetfeki boundary condition [4].

C. W. Wu, H.X.

Sun I J. Non-Newonion

Fluid bitch.

56 (1995) 253-266

261

eccentricity ratio E = 0.8. The lubricated zone 0 = 0 -+ 71is divided into 30 pressure elements with 31 nodes. One pressure element contains 14 sub-elements: therefore, the total number of sub-elementsis 14 x 30 = 420. In order to compare the present work with Tichy’s work [4] for the Bingh::m model, the computed results are given using Tichy’s dimensionless parameters: rg = r,c/(~oR), P* =pc2/(qoR2) = PAL/ (I - E)~.Here. the length of the bearing is B = RR. From Fig. 5 it can be shown that when the viscosity ratio, 5, increases, the film pressure decreases. When 5 = 100, the film pressure distribution of the bilinear rheological model is so close to that of the Bingham rheological model [4] that a diffcrcnce can hardly be found. When rt = 0, the bilinear rheological model becomes a Newtonian model with viscosity q. Secondly, in the outlet zone (n < 0 < 2x), the Reynolds pressure boundary condition is used to analyze the cavitation lubrication. In mathematical language, p = 0 at the inlet 0 = 0. but p = 0 and ap/at? = 0 at some point 0 = OunknDWn in the outlet zone, with Ountnawn lying between A and 2n not known a priori. This is a free boundary problem resulting in some difficulties in numerical analysis. A detailed discussionabout the Reynolds boundary condition can be found in many standard texts on lubrication theory [6]. The pressureprofiles under the Reynolds boundary condition are shown in Fig. 6 for different values of 5 using the present method. The analyzed zone is 0 = O-, 38n/30 consisting of 38 pressure elements and 14 x 38 = 532 sub-elements(see Fig. 7). That is to say. in the numerical analysis,we can try to take tue pressure boundary condition p = 0 at both 0 = 0 and 0 = 38n/ 30. From Fig. 6 it can be seen that the point at which the positive pressurevanishes

Fig.8(a)

C. W. WU, H.X.

Sun I J. Non-Nenronian

Fluid Mrch.

56 (W95)

253-266

(b Fig. 8. (a) Velocity profiles under the half-Sommerfeld boundary condition CC = 100. r: = IO, E = 0.8). (b) Velocity profiles under the Reynolds boundary condition (I = 100, rg = 10, a = 0.8).

e Fig. 9. E&t

of yielding

shear stress on the pressure distribution

(C = 100. t = 0.X).

C. W. Wu. H.X.

Sun I J. Non-Nerronion

Fluid Mech. 56 (1995) 253 -266

263

0 0

V6

+a

%lS

aq3

6116

%

WI?

8 Fig. 10. EtTect of eccentricity

rtio

on the pressure

distribution

(< = 100. rl = IO).

is at 35x/30 (i.e., fIunknom= 35x/30) for 5 = 10 and at 34x/30 (i.e., fIunLnown = 34x/ 30) for 5 = 100. Theoretically speaking, the terminated analysis point can, therefore, be taken at any value from 34x/30 to 2x when 5 = 10, or from 35n/30 to 2n when t = 100. Fig. 7 shows the yielding states when t: = 100 and E = 0.8 for the half-Sommerfeld and Reynolds boundary conditions, respectively. It is found that the so called “core” zone (no yielding) for the bilinear rheological model under the half&mmerfeld condition is almost the same as the right half of that for the Bingham rheological model under the Sommerfeld condition. In the inlet zone, however, the yielding zone (ITI> q,) size for the Reynolds boundary condition is a little smaller than that for the half-Sommerfeld boundary condition. In the outlet zone the situation is opposite to that in the inlet zone. Figs. 8(a) and 8(b) give the velocity profiles for the half-Sommetiftild and Reynolds boundary conditions, respectively.At the cross-sectionsof 8 = 0 and x/2 the velocity profiles for the two boundq conditions are near% identical, bui in the outlet zone, some deviation exists. Fig. 9 shows the etfect of the yielding shear stresson the pressuredistribution for e = 100 and E = 0.8. Fig. 10 shows the effect of the ec&nttLity ratio E on the pressure distribution for e = 100 and ~0’= IO. Both Figs. 9 and IO were obtained using the same elements as used in Fig. 6 under the Reynolds boundary condition. It can be seenthat the effect of the eccentricity ratio E is more important than that of the yielding shear stress.

264

C.W.

Wu. H.X.

Sun i 1. Non-Neu~fonim

Fluid Mech.

56 (1995) 253-266

The hydrodynamic lubrication problem of a bilinear rheological fluid is studied using the FEM and a mathematical programming method. The comparison of Tichy’s work on the Bingham model and the present work shows that the numerical method described in this paper is a reliable and general one for bilinear rheological fluid hydrodynamic lubrication. The theory and numerical method given in this paper can be used to study the general non-Newtonian flow of bilinear rheological fluids in a journal bearing The Bingham plastic fluid is only a limiting case of bilinear rheological fluids. A unified Reynolds tylx equation was obtained. It can be easily used to study both squeezing flow and shear flow of bilinear rheological fluids. The cavitation hydrodynamic lubrication problem of the bilinear rheological model was also first worked out using a mathematical programming method. Taking the Reynolds pressure boundary condition as an example. the hydrodynamic lubrication of a journal bearing is studied. It is shown that there is a considerable difference between the pressure profiles using the Reynolds and half-Sommerfeld pressure boundary conditions, but there is little difTerence between the velocity profiles. In lubrication mechanics, the Reynolds pressure boundary condition is widely applied to the analysis of journal bearings. The half-Sommerfeld pressure condition is rarely used because the flow continuum equation can not be satisfied at the point of the minimum film (here ap/ax is not continuous, see Fig. 5). The Sommerfeld boundary condition is almost never used in engineering applications because no fluid support force will be obtained. Wats et al. [12] and Tichy [4]. respectively, used the half-Sommerfeld and Sommetfeld pressure boundary conditions to study Bingham plastic flows. Using their methods to deal with the Reynolds boundary condition would, however. be difficult. In our method the Reynolds boundary condition is automatically satisfied in the solution process. No additional work is necessary. If t is larger than 100. the bilinear rheological model seems to approach the Bingham rheological model. There is reasonable agreement between the present results (6 = 100) and Tichy’s analyses [4] on a Bingham plastic fluid; therefore, the present method can also be used to analyze a Bingham fluid lubrication problem.

Achwiedgment This work Foundations,

was jointly China.

supported

by the National

&Ieadix: MaMces aad vectora ia Eqs. (Ml)-(14c)

and

Liaoning

Natural

Science

C. W. Wu. ti.X.

Fig. 41. Schematic

Sun / J Non-Newtonian

of the eth pressure

Fluid Mech.

clcmcnt

56 (1995) ?5j

and the r&red

-766

265

suh-elemencs

where L.,,.is the dimensionless length of the eth pressure element. When uniform nodes are used in the film thickness (y direction). [D] = i

[ a!!

c n I JI.,.

[D,., - D,.] dX.

[D,,] = [D,.,,D,+ . . .D,.,l.

al

[Cl = [Cc, G,.

.

.C:.,IT.

where f4,, is the Y coordinate value of the shape center of thejth sub-element rn the eth pr,s;;re element and fi,, is the dimensionless film thickness corresponding to the shape center, as shown in Fig. Al.

265

C.W.

We,

[~ 34 S u n

/ d. Non-Pl(,~tottlan

Fluid

I~l,.ck. 5 0 ( 1 9 9 5 )

253

260

w h e r e E ~ i s a v ~< r s q u a r e m a t r i x o n l y w i t h t t n i { e l e m e n t s i n i t s d i a g o n a l a r × v square matrix in whiGh all the cletnents a~e unit elements,

u~

_

Y * - - U*'

{F:'}~(--~2I'~. {A}

{A;'? .....

~~r/ "

-- u,*

~)

AW, A;¥ .....

--Yon'*'*"

[1, 1. . . . AW .....

[1

I

and

!o i s

. .

,I]'. A,?J . . . . .

A ; £ ). . . . .

AL=,,) . . . . .

AL~2} - .

Refereoecs [ l [ S.12).I~. W i l s o n , S q u e e z i n g I l o w o f a B ] n g h a m m a t e r i a l , J. N o n - N e ~ l o n i a n F l u i d M e c h . , 4 7 (1983) 211--219. [2] $ . E . S ~ a ~ g r o o m , Tl~e B i n g h a m p l a s t i c mode~ Q~ E R f l u i d s ~ n d i t s i m p l i e a l i o n , P r o c e e d i n g s o f t h e S e c o n d I n l e r n a t l o n a l C o n f e r e n c e o n E R F l u i d s , T e d m i ¢ , L a n e 0 s l e r . P A . 1990, p p . 41--52. I3] A . D i m a l ' o g o n a s a n d A . K o l l i a s , E t ~ t r o r h e o l o g i c a l f l u i d - c o n t r o l l e d ' s m a r t " j o u r n n l b e a r i n g . T r i b o l o g y T r t m s . . 35 (1992) 611--618. [4] J . A . T i e h y , l t y d r o d y n a m i c l u b r i c a t i o n tl~t~ory f o r t h e B i n g h a m p l a s t i a f l o w m o d e l , $. R h ¢ o l o g y , 35 (1991) 477--496. [ ~ B.O. 3 a e o b s o n t~nd • J . F l a m r o e k , N o n * N e w t o n [ a n f l u i d m o d e ] i n w r p o r a t M i n t o ¢ ] a s t o h y d m d y n a r a i ¢ ] u b r i ~ t i o n o f t ' ~ t a n g a l ~ r ¢ o n l a e t s , A S M E J. Tri[~ology, 106 ( I 9 8 4 ) 2 7 5 ~ 2 8 2 . 16] A . C a m e r o n , B a s i c L u b r i c a t i o n T h e o r y , 3 r d e d n . , E l l i s H o r w o o d L t d . , Clliehe~t~r, 1981. 17) K.~I. H n e b n e r , T h e F i n i l e El~:txl~nt M ~ l h o d f o r Etllginc~.'rs , J. W i l ~ y , N e w ~(ork, 1995 [8] W . X . Z h o n g a n d C . W . W n , E l a s t i c p l a s t l ¢ e o o t ~ t s u s i n g p a r a m e l t i c p r o g r a m m i n g , i n M . ] ] . A l i a b ~ d i a n d C . A . B ~ b b l a ( E d s . ) , C o m p u t a t i o n a l m e t h o d ~ i n C o n t a c t ptx~blelns; C o m p u t a t i o n a l M e c h a n i c s P u b l i o a t i o n s L t d . , ]993, C h a p t ~ 9, p p . 305 352. 19] C . W . W u e t aL, P a r a m e t r i z v a r i a t i o n a l p r i n c i p l e o f viscu3plastta l u b r i ~ t i o n m o d e l . A S M E J. T r i b o l o g y , 114 (1902) "/31--935. I10] ~ D . C h a n g . B a a i c T h e o r y o f E u g i n e e r i n g S t m t n t u r z Optirnlv~ttiott D e s i g n , ~Vate r C o n s e r v a n c y a n d I ] l ~ t r l e i t y p t l b l i s h ~ r , B¢ijln~, C h i n a . 1984, [ i n C h l n e s ¢ ) . [ 11] C . W . C r y e r , T i l e n l c t h o d o f C h H g t o p h e r s o n f o r s o l v i n g f r ~ b o u n d a tW p r o b l e m s f o r i n f i n i t ej o u r n a l h e a r i n g s b y m e a n s o r f i n i t e d i f f a r ~ n ~ . M a t h . C o ~ l p u t . , 2S (1971) 435 443. 112] S. W a d a e t al., B e h a v i o r o f a B i n g l x a n x~ l i d i n h y d r o d y n a m i c l u b r i c a t i o n ; p a r t 1, G e n e r a l tkeor~d, B u l l . J S M E , 16 ( 9 2 ) ( ] 9 7 3 ) 422 431.