A p-version finite element gas bearing program for pivoted calculations based on false transient response

FINITE ELENENTS IN ANALYSIS AND DESIGN Finite Elements in Analysis and Design 16 (1994) 1-14

ELSEVIER

A p-version finite element gas bearing program for pivoted calculations based on false transient response Steve H. N g u y e n * The University of Kansas, Mechanical Engineering Department, Lawrence. KS 66045, USA Received April 1992; revised July 1993

Abstract

An efficient finite element program for pivoted gas bearing design calculations, called p-GAS, is presented. The program is constructed based on the modified, steady-state, compressible lubrication Reynolds equation and static equilibrium equations, p-GAS provides pivoted solutions by tracking false transient responses. The present gas bearing finite element program yields accurate numerical solutions. It conveniently provides higher-order solutions since this program utilizes a hierarchical finite element formulation, and it requires less computing time than does even the most efficient known-to-date factored implicit finite difference scheme. The success of p-GAS rests on three factors, namely, the hierarchical p-version finite element method, which reduces the total degrees of freedom needed for modeling the bearing systems; the artificial transient approach which rapidly brings the bearing attitude to the correct attitude; and the Gaussian frontal solver which accelerates the matrix solving time by reducing the in-core storage requirement. Experiments with reduced integrations, which save the computing effort, yield acceptable results. Directions on further reduction of the computations are briefly discussed.

1. Introduction

The pressure generated due to fluid flow between movable bearing surfaces either pushes the bearing surfaces apart (positive pressure) or pulls them closer together (negative pressure). To maintain proper clearance, forces are applied to the bearing system. Depending on how the force parameters are handled, one of two types of static calculations is performed. The first type, called the fixed-attitude calculation, computes the required force and its location on the bearing for achieving a given clearance or spacing. However, it is of more interest to solve the inverse of this problem or the second type, called the pivoted calculation, which computes the attitude once the

*Currently at the University of Texas Health Science Center, School of Dentistry, Houston, TX 77225, USA. 1994 Elsevier Science B.V. SSDI 0 1 6 8 - 8 7 4 X ( 9 3 ) E 0 0 7 3 - A

2

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1-14

specified applied force and its location are given. Pivoted calculations are more elaborate because for a given spacing, the corresponding gas bearing lift and moments generated in the bearing channel often do not balance the loading condition on the bearing• As a result, the spacing must be continually adjusted and pressure calculations are performed based on the new spacing until the balancing condition is achieved. The final spacing can be largely different from the starting one. Pivoted calculations are used very frequently in the analysis and design processes since constraints often introduce "what-if" questions. In computer industries, these two types of calculations are heavily exploited because compressible lubrication plays a very important role in the design of any magnetic storage interface• This interface consists of a runner (rotating disk or moving tape) serving as a storage medium, and a slider acting as a floating vehicle on the magnetic medium for writing information onto and reading information from the magnetic surface. Fig. 1 illustrates a slider used in today's magnetic disk drives and an applied force on the slider by means of a suspension. This is the well-known Winchester slider, having two bearing rails for supporting the gas lift. On each rail a tapered region near the leading edge acts as a pressure stabilizer. The gas pressure acts on the two rails and makes the slider "fly" above the magnetic disk a distance He, a pitch angle ~, and roll angle /3, as shown in Fig. 2. These three quantities are often referred to as the slider attitude. Although this bearing configuration looks rather simple, there is no analytic solution describing the pressure distribution on the rails. Bearing designers must resort to numerical methods. In this computational arena, the finite difference (FD) method has played a dominant role [-1-12] in the analysis and design of modern magnetic slider bearings even though many FD schemes are limited to only rectangular bearings within which load-supporting regions must be parallel to the main bearing boundaries, while the finite element (FE) method has been very much under utilized. It is clear in the lubrication literature [,13-18] that the FE method is very general, and that it poses no difficulty with irregular bearing geometries, complicated clearance, and mixed boundary conditions. The preference for the FD method can be explained by three reasons. First, the FE

Suspension

Yf

0

Moving Disk Surface

(XL, YL)

) x

XT

U u

u = I U I cos (0) v=-IU

Fig. 1. Physical slider-disk configuration and plan views.

I sin (0)

S.H. Nguyen / Finite Elements in Analysis and Design 16 (1994) 1-14

Fig. 2. Free-bodydiagrams. method for compressible lubrication studies did not emerge as a general tool until 1974 [17]. By this time FD methods had been intensively tested. Second, the convenience in describing the discretized bearing mesh by the FD methods often wins the analyst's heart. And finally the short computing time needed by FD schemes further encourages this economical direction. With the fastest FD scheme known-to-date being the factored-implicit FD scheme (FIFD) developed in 1980 [4], modern bearings such as the 3370 and 3380 Winchester sliders can be solved very rapidly since the FIFD scheme operates in a splitting direction manner. As a result, FE methods were occasionally used for analyzing bearings that are difficult to model by the FD scheme [14,16,19]. With respect to the first reason, there is nothing one can do to undo history, but there are numerous improvements that can be done to bring the FE approach to a competitive level. Today, software that provides automatic meshing procedures and geometry definitions has removed the mesh generation burden (second reason). The last reason dealing with the computational speed has received much attention. In 1984, Garcia-Suarez et al. [19] implemented the upwind formulation which helps speeding up the calculations. In 1991, Crone et al. [20] used parallel processing to solve a low-order FE formulation. In an effort to reduce further the required computations to make the FE method more attractive, Nguyen recently proposed p-version hierarchical FE formulations [21-24]. Nguyen [24] was able to model a low flying Winchestertyped slider having side flows by fewer degrees of freedom without affecting the accuracy and stability of the numerical solutions for fixed-attitude analysis, thus demonstrating large reductions in the computational effort. In the present investigation, a gas bearing program, which utilizes the hierarchical FE formulation proposed in [24], is constructed for modeling modern magnetic slider designs. Emphases are placed on computing times and accuracy for pivoted analysis. The program, called p-GAS, computes steady-state pivoted solutions by tracking false transient responses. Comparison is then made between results obtained from p-GAS and the FIFD scheme. The p-FE results are also compared against experimental and other numerical results when possible. Essential components of p-GAS and the underlying theory are presented.

2. The new gas bearing program

The gas bearing program developed here provides pivoted steady-state solutions, making use of the compressible hierarchical FE formulation presented in [24]. This section first introduces the governing equations. Then, it describes basic subroutines used in p-GAS together with the

4

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1-14

numerical methods. Finally, a numerical procedure, which outlines execution steps performed in the program, is given. The approach used is called the false transient method because no timedependent equations are considered at all. 2.1. Governing equations.

The pressure profile P generated due to fluid flow between the bearing surfaces separated by a clearance /4 is described by the modified Reynolds lubrication equation (see [25] for more details). Here the steady-state, isothermal, compressible lubrication equation that includes the molecular mean free path slip-flow effect was used: V . ( H Z ( P H + 62aPa)VP - @ U P H ) =

O,

on (2,

(1)

on C1,

(2a)

on C2,

(2b)

with the boundary conditions P = Pa, ( H 2 ( P H + 62~Pa)VP- @ U P H ) . n

= q,

where p is the gas viscosity, 2a the gas mean-free path, Pa the ambient pressure, U the runner velocity vector, C1 and C2 constituting the bearing boundary, ~ the bearing domain (bearing area) and n the outward normal unit vector. Here the bearing surfaces are the magnetic slider/disk interface. The positive pressure present in the slider/disk bearing pushing the slider away from the disk surface is balanced by a force F applied on the slider at x v and Yv to prevent the slider from "flying away". When summing forces in the vertical or spacing direction H, and summing moments about the corner point B (see Fig. 2), we obtain W-

(3a)

( F + rag) = Fo,

W ( x L -- xe) -- F(xL -- XF) -- mg(xL -- Xc) = F I ,

W y p - F y v - mgYc = F 2 ,

(3b)

where m denotes the slider mass, 9 the gravitational acceleration, xL the slider length: Fo, F~, and F 2 the unbalanced force and moments, respectively, Xc and Yc the slider's center of mass, Xp and yp the center of pressure, and W, Wx, Wy the pressure lift and moments, respectively. The pressure lift, moments and center of pressure are computed as tw, wx, wy) =

Ii,x,y) (P - Pa)dQ,

x, = wx/w,

= wy/w.

t3c)

The probability of selecting a spacing distribution, H, in (1) so that Fi (i = 0, 1, 2) are negligible, i.e. the slider is in equilibrium, in one attempt is extremely minute. Very often (1)-(3) are solved repeatedly until the unbalanced terms vanish, and in this process H is continually adjusted. 2.2. Basic subroutines

In this section, essential subroutines are introduced, and then a numerical procedure which integrates these subroutines is presented, p-GAS consists of three main subroutines: INPUT: Organization of all input parameters resides in this subroutine. Data are read in, including the slider geometry such as length, xL, rail width, b, bearing width, YL, and the taper

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1 14

5

angle, 7, located at the leading end (x = 0); the applied load, F, and its location, xr and Yv; the gas parameters :/~, 2a, and Pa; initial guesses on the spacing HB at the trailing end (x = XL), pitch angle c~,and roll angle/3; the bearing discretized information; and the runner velocity vector U consisting of two components, u in the slider length direction x and v in the slider width direction y as shown in Fig. 1. If the slider is positioned at an angle 0 with respect to U in such a way that the slider's trailing edge is closer to the spindle than the leading edge is, this is a negative skew angle (0 < 0) by convention. If the leading edge is closer to the spindle, 0 > 0. Fig. 1 demonstrates a negative skew. SPACING: To establish a systematic way of describing the slider/disk spacing distribution, a quantity called the bearing contour, Ho, is introduced. Ho provides information on the slider features. First, consider the slider which consists of a flat slab resting on the z = 0 plane, i.e. H0 = 0. The taper-flat slider configuration in Fig. 1 calls for a cutout portion in the taper region. Since there is no need for specifying the nonload supporting recess portion between the slider rails, Ho is simply expressed as Ho=(xT--x)tan7

ifx~
and

H0=0

ifX>XT

for all y values. For other contours, H0 will be modified further. This initial spacing H for this taper-flat slider is then computed as H = Ho + HB + (XL -- X)a + yfl.

The spacing distribution or clearance varies continually as the slider approaches the equilibrium state. Spacing corrections are done by changing the spacing at B, pitch and roll angles as follows: (4)

H = Ho + (Hn + AH) + (XL-- X)(a + A~) + y(fl + Aft).

In this form, (4) is more convenient to program. It rapidly provides H(x, y) once the spacing at B, pitch and roll angles are known. Note that Ho is precomputed once and is reused in every iteration. REYNOLDS: This subroutine computes the pressure information based on the spacing distribution (4). Due to the nonlinear nature of the compressible Eq. (I), FE solutions must be obtained iteratively as described in [17-19, 21-24]. The formulation used in p-GAS is taken directly from [24], originally derived for the fixed attitude analysis, in which the solution P is repeatedly corrected until convergence in P is reached for a fixed spacing. In this process, the pressure correction term 7/in (5) is solved based on the most recent pressure solution P, and the pressure solution is updated, i.e. Pt,ew~ P(old) + ~. Then the computational process repeats the same cycle again until convergence criteria are reached. The fixed attitude calculation process is summarized in Fig. 3. For the pivoted analysis presented in this paper, the pressure correction term is computed =

~

Spacin,_g[I R

e ~ o l d s ~

" I (Frontal~ - ~ , . o , , , ~ s ~ , , ~ , - ~ k

Fig. 3. A fixed-attitude calculation program flowchart.

o,t,vj

6

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1 14

once for each spacing, i.e. no convergence in P for any given spacing H. This is the false transient portion of the calculation. H is continually adjusted based on P, and for each pressure distribution P, one correction seems adequate: [K]e{qJ}e

=

{F}e,

(S)

where

D=

H 2

(HP + 6,~aPa).

P(~,q) = L N,(~,q)4', = [N] {4'}e,

~u(~,~/)= [N] {qJ}e,

(Sc)

(x,y,u,v) = EL] {(x,y,u,v)}e.

(5d)

i=1 l

H(~,q) = ~ L~(~,q)H, = EL] {H}e, i=1

The subscript e implies the element level, and the comma denotes partial differentiation with respect to the variable that follows. {¢}e and {qJ}¢ are the elemental generalized pressure and incremental pressure vectors. They are termed "generalized" because they contain not only nodal quantities but also higher derivative values. Here r denotes the total degrees of freedom per element, ~ and r/the natural coordinate variables varying between - 1 and + 1, and l the total nodes configuring the element geometry. While the Lagrangian interpolation functions L(~, q) in (5d) are well-known and are available in any standard FE text, the hierarchic basis functions N(~, q) are not as popular because of their late arrival. N(~, r/) for quadrilateral elements are obtained by taking tensor products of the one-dimensional hierarchic functions N(~) and N(q) (see [21] for derivations). Components of N(~), N(r/), [N(~, r/)] and {¢}e are given in Table 1, where p~ denotes the polynomial level in the element length direction, and p, the polynomial level in the element width. [Nx] and IN.y] in (5a) and (5b) are related to [ N j and [N~] by the Jacobian mapping matrix derived directly from the element coordinates, x and y, given in (5d). The formulation (5) yields a non-symmetric FE system of equations. Based on this formulation, first {4'} is assumed to be the ambient pressure. Eq. (5) at the global level provides the global generalized incremental pressure vector { 7j} when solved. Because the global fluidity matrix [K] can be very large, it is not feasible to start the matrix solving procedure after all elemental matrices are assembled. The elemental matrix and vector are reduced on element-by-element basis by the Gaussian frontal algorithm. The details of this solver are documented in [26]. Briefly, instead of solving for the global vector {7j} after all element matrices are assembled, the frontal solver starts reducing the element matrices on the element-by-element basis in REYNOLDS. Equations that cannot be reduced immediately remain in-core, and those fully summed are reduced and stored off-core on disk. This procedure minimizes the in-core storage requirement, thus speeding up the solving time.

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1-14 Table 1 Hierarchic C ° shapefunctions and nodal variables for quadrilateral elements Location

N(~, t/)

q~ P

-- l, -- 1 0, -- 1

No(t/)No(~)

No(rl)N~(O

(s = 2,3 . . . . . pc)

l, -- 1 1,0 1,1 0,1 - 1, 1

No(q) Nl(~) Nt(tl)N,(O

(t = 2,3 . . . . . p~)

N~(tl)Ns(¢) Nt(t/)No(~)

(s = 2,3 . . . . . Pc)

-- 1,0

X,(tl)No(0

(t = 2,3 . . . . . p~)

Ss(tl)S,(O

(s = 2,3 . . . . . p~), (t = 2,3 . . . . . p~)

P P

Xt (t/) X,(~)

0,0

P

N 0 ( o ) = (1 - to)/2, Nl(to ) = (1 + 09)/2, N o ( o ) = (090 - ae(to))/O! 1, ao(to) =

if 0 = even,

to,

dts+0p P ~',I' -

if 0 = odd,

~(s dr/'

N E W T O N : In this subroutine, new flying attitude is determined, and convergence is checked. Once the pressure correction vector { ~} is available from R E Y N O L D S , the pressure is corrected; P(new) = P(old) + t/J

=,,

{~b}(new)

=

(6)

{~D}(old)"-F {~[/}

and the gas bearing lift, moments and center of pressure are c o m p u t e d according to (3c). A new slider attitude is estimated based on error ratios derived from equilibrium Eqs. (3a) and (3b), respectively, as rn=

r~ =

(w (

F-TFmg

)

1 ,

(7a)

)

W(XL -- X~) -- 1 F ( x L -- x v ) + m g ( x L -- XC) '

r~ =

)

1 , F y F + mg Yc

(7b)

where (7a) is obtained by dividing both sides of (3a) by (F + mg), and (7b) is derived similarly from (3b). These ratios indicate the "distance" of the current solution from the equilibrium position. The increments in H, ~ and fl are estimated as Ai = iasgn(ri) ,

i = H,~,fl,

for Iril >1 0.40,

(8a)

Ai = 2.5iasgn(ri),

i = H,~,fl,

for Iril < 0.40,

(8b)

where the subscript a denotes allowable values, and the symbol sgn indicates the sign of the argument inside the parenthesis. The computational process stops, and the slider attitude (Hn, ~ and fl) is accepted as the correct solution if and only if the following criteria are satisfied: Iril <%e,

i = n,o~,fl.

(9)

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) I 14

1(Fr°ntal) I

I

[k=k+l[

No

[

Fig. 4. The current pivoted-calculation program flowchart.

Otherwise, information from (8) is then used for adjusting the clearance H via (4) before another pressure calculation loop begins.

2.3. Numerical procedure Having described the mathematics and algorithms used in each subroutine, we are now at the point of integrating all subroutines together to yield the desired gas bearing program. Operations performed in p-GAS are outlined graphically in Fig. 4. Clearly, in every pressure calculation loop k, the clearance H (4) is established in SPACING, the pressure correction { ~} and the pressure {P} are computed based on (5) in REYNOLDS and the frontal solver, and N E W T O N computes new flying parameters (8) and checks for convergence (9). Within each loop, {P} is corrected once. As soon as (9) is satisfied, the computation process stops.

3. Numerical experiment In this study, uniform p levels were selected (p = p~ = p,) throughout. [K]e was computed exactly by Gauss-Legendre quadrature with (p + 1) sampling points in each element direction. Since the size of the element matrix [K]e in (5) increases as the order of approximation goes up, it is clear that the computational effort devoted to the calculations of [K]e could become tremendous if very high p levels were chosen. To reduce the computing time, first, more elements were chosen to compensate for the selection of low p levels; second, FE meshes were numbered so as to minimize the matrix front width; and third, the converged solution the slider attitude and pressure profile for a given p level was used as the starting solution for the next higher p level calculation. To shorten the computing time further, we used reduced integrations which require one less sampling point in each element direction for a given p level. Note that reduced integrations are often used for other purposes [27]. Since p-GAS conveniently provides higher-order solutions based on the hierarchical basis functions given in Table 1, p-FE models do not require as many elements as do other models. The FE mesh used consists of 60 elements (15 unequal spacings in the length or x direction, and four equal spacings in the width or y direction) as shown in Fig. 5(a). Near high pressure gradient areas such as the trailing edge and the taper end, grids are more compact to prevent instability. Computed solution are then compared with published simulations and experimental results, if available. To obtain the computing time required by the FIFD model, a computer program using

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1 14 IlU[ I I lml I

ILIIIIIIIIIIIIIIIIIIIllll ::

®

IIIII I I lUl I

9

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

:

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

:

®

Fig. 5. (a) p-FE mesh, (b) FD mesh. the F I F D scheme similar to the one presented in [8,9] was written for this project for comparison purposes. The F D mesh used consists of 66 points and 49 points in the x and y directions, respectively, shown in Fig. 5(b). This is a typical F D mesh, which is unequally spaced and dense to account for the rapid change in pressure for high bearing number flow conditions. Throughout this study, the compressible gas was selected to be air at 25°C, SI units were used, and the following data were applied to all calculations, unless indicated otherwise: Pa = 1 0 1 3 2 0 ,

P = 1.84 x 10-5,

= 40 x 10 - 6 r a d , fla =

10x 10-6rad,

fl = 0,

H B = -- 350x 10 -9,

'~a = 65 X 10 - 9 ,

e = 0.005,

H a = 40 x 10 - 9 ,

~a = 25 x 10 - 6 rad,

9 = 9.806.

p-GAS was written in double precision fortran format. The program ran on an engineering workstation, Apollo series 400, and p-GAS obtained the required C P U time by calling CTIME, and IMSL subroutine. Example 1: A 3370-type Winchester slider

A taper-flat Winchester 3370 slider was selected. The following slider dimensions and force parameters were used: a = 0 . 6 3 5 x 1 0 -3, x r = 0.380 x 10 - 3 ,

x L = 4 . 0 6 0 x l 0 -3,

y L = 3 . 0 5 0 x l 0 -3,

7 = 17.453 x 10 . 3 rad (or 1°),

Y c = Y v = 1.525x10 -3,

F=0.146403,

m=7.0xl0

Xc = x r = 2.030 x 10 - 3 ,

-5

To bearing modelers, this is called the 15 g load case because the applied force F and the slider weight total 0.14709 N or 15 g. Case 1: N o skew condition (0 = 0 °)

The slider center of mass (Xc, Yc) is positioned at 0.0381 from the center of the spindle rotating at 60 revolutions per second, which implies u -- 14.36336 and v = 0 at the slider's center of mass. The disk velocity components vary depending on the radial distance from the center of rotation ("small disk" effect). The computed spacings H at point B (x = xL, y = 0), pitch and roll angles, and the computing time needed are compiled in Table 2(a). To demonstrate the convergence of the results, solutions of four p levels are tabulated. At the p level one (or linear approximation) within each element, the slider attitude was not accurate as expected from such a coarse mesh. Improved results

S.H. Nguyen/ Finite Elements in Analysis and Design 16 (1994) 1 14

10

Table 2(a) 3370 slider calculations for a zero skew position (0 = 0 °) Type

p

HB X 10 9

~ X 10 6

fl X

10 6

FE

1 2 3 4

245.7 309.6 309.9 309.6

30.8 46.9 47.4 47.4

5.9 7.0 6.9 6.9

(5.7) (6.9) (6.8) (6.8)

p-model

(256.0) (310.1) (309.6) (309.5)

311.9 a 311.0 b

F D models

(53.8) (48.3) (47.5) (47.4)

50.3 48.0

Time ( C P U s) 1.05 3.82 1.74 1.24

(0.86) (3.14) (1.11) (1.05)

6.7 6.8

11.61

Time (CPU s)

aRefs. [8,9]. bFD mesh shown in Fig. 5.

Table 2(b) 3370 slider calculations for a negative skew position (0 = - 12 :') Type

p

H B x 10 9

~xx 10 6

flx 10 6

FE

2 3

232.0 (233.0) 228.7 (228.4)

62.5 (63.5) 63.9 (63.8)

12.8 (12.6) 13.2 (13.2)

231.7" 234.2 b

61.2 63.3

17.2 12.6

p-model F D models

5.21 (4.39) 5.17 (5.08) 12.98

"Refs. [8,9]. bFD mesh shown in Fig. 5.

are obtained with the p level two (quadratic), p level three (cubic), and p level four (quartic). Changes in the slider attitude between the quadratic, cubic, and quartic solutions are very small, indicating convergence of the results. From Table 2(a), the quadratic approximation seems adequate; however, the cubic results are often needed for assurance. This is equivalent to refining the mesh and rerunning the model often seen with FD and the conventional h-version FE methods. Pressure profiles beneath the bearing pad (not shown) do confirm the convergence of the pressure profiles as the p levels go up. For the choice of mesh selected (60 elements), there is clearly no need to run very high-order solutions. The pressure profile for the quadratic approximation (p = 2) appears adequate. When compared against other results, the p-FE results agree well with FD results. With regard to the computing time, FD calculations utilizing the fastest known-to-date FIFD scheme consume more time than p-GAS. To further shorten the required computing time, reduced integrations were selected, and results are reported in parentheses. Case 2: A negative skew conditions (0 = - 12 °) For this case, u = 14.04949 and v = 2.98631. The computed results are given in Table 2(b). Changes between the quadratic and cubic solutions are very small, indicating convergence of the results. Agreement between p-FE and FD results is clear, except a slight deviation in the roll value

S.H. Nguyen/ Finite Elements in Analysis and Design 16 (1994) 1-14

11

reported in [8, 9] is observed. This is probably due to the coarser mesh (56 x 42) used in [8, 9]. Also, for a large skew angle, side flows become significant, and we model the slider by unequally spaced F D grids in both x and in y directions to account for the rapid change in the pressure profile, while in I-8,9-1 grids in the y direction were uniform. It is remarkable that p-GAS again takes less computing time. Results from reduced integrations are reported in parentheses. Example 2: A 3380-type Winchester slider

In the second example, a taper-flat Winchester 3380 slider is modeled. The following slider dimensions and force parameters were used: a = 0.400 x 10-3, 7=

XL = 4.000 X 10-3,

15.0xl0-3rad,

m=5.0x10-5,

YL = 3.000 x 10-3,

Xc = XF = 2 . 1 O O x lO - 3 ,

X T = 0.40 X 10-3,

y c = YF = 1 . 5 0 0 X 1 0 - 3 ,

0 = 0 °.

Case 1: A normally loaded slider (F = 0.14660) The disk velocity components are fixed, u = 40.0, v = 0.0, everywhere in the bearing channel. The c o m p u t e d results, including the slider attitude and C P U time, are presented in Table 3(a).

Table 3(a) 3380 slider calculations for a normally loaded and zero skew case Type

p

H B x 109

~ × 106

FE p-model

2 3

285.3 (283.9) 285.2 (284.7)

104.7 (104.7) 105.4 (105.5)

295.0 ~ 284.8 b 290.0 c

106.8

Other results

fl × 106 -0.1 (-0.1) 0.0 (0.1) 0.0

Time (CPU s) 4.46 (3.73) 1.52 (2.54) 12.85

"FE result 1-20]. bFD mesh shown in Fig. 5. CExperimental result [10]. Table 3(b) 3380 slider calculations for a heavily loaded and zero skew case Type

p

HB × 109

~ × 106

fl × 106

FE p-model

2 3

127.9 (129.8) 129.0 (129.3)

11.5 (13.8) 13.6 (13.6)

-0.1 (-0.1) 0.0 ( - 0 . 1 )

129.0 a 130.4 b

14.5

Other results

aExperimental result [20]. bFD mesh shown in Fig. 5.

0.0

Time (CPU s) 6.18 (5.13) 2.55 (1.82) 13.15

12

S.H. Nguyen/Finite Elements m Analysis and Design 16 (1994) 1 -14

Convergence of results is evident from Table 3(a). The computed results agree with published FE results [20], experimental results [10], and with FD calculations. Slight discrepancy between the current result and the published FE results reported in [20] is observed. This is expected since [20] modeled the slider with low-order triangular elements. With regard to the computing time, p-GAS continues to take the lead in economizing the required CPU time. Results from reduced integrations reported in parentheses deviate slightly. Case 2: A heavily loaded slider (F = 0.9806) The disk velocity components are constant everywhere, u = 44.5, v = 0. Note that a very large force is applied on the slider (to gas bearing designers, this is called the 100 g load case), pushing the slider much closer to the disk surface. To prevent large excursions in the false transient response due to this low flying condition, maximum allowable increments (H a, ~a,/3a) in the slider attitude are halved. For FD calculations, this corresponds to the halving of the time step. Computed results are given in Table 3(b), showing agreement between the current solutions, FD and experimental results. Less computing time is used by p-GAS as evident from Table 3(b). Reduced integrations continue to provide acceptable results.

4. Summary Modern bearing design solutions obtained by a gas bearing finite element program, called p-GAS, are presented and compared against FD and other published results. The program is constructed based on the modified steady-state lubrication equation and slider's static equilibrium equations, p-GAS yields steady-state pivoted solutions by tracking the false transient response until variations in the bearing lift and moments are below some prescribed tolerance. For the 3370 and 3380 Winchester-type sliders being tested here, the present gas bearing finite element program clearly shows that the program provides accurate results and that it requires less computing time than does even the most efficient finite difference scheme known to date. The success of p-GAS relies on three factors, namely, the hierarchical method which reduces the total degrees of freedom needed for modeling the bearing systems; the artificial transient approach which rapidly approaches the correct attitude; and the Gaussian frontal solver which accelerates the matrix solving process by taking advantage of the off-core storage and therefore decreasing the in-core storage requirement. These three factors substantially reduce the computational effort. Here we were able to keep the CPU time reasonable by using quadratic approximations and checking the results by cubic approximations. One can use the current program for obtaining very high approximation results, but this is not recommended since there is a point beyond which any further increase in the approximation level only raises the computing cost. Effort to build FE meshes is much less here since reasonably coarse meshes can be selected. Also p-GAS conveniently outputs higher-order solutions, which provide a clear picture of the convergence process of the slider attitude. Experiments with reduced integrations show acceptable results. The computing time would be further shortened if elemental matrices were condensed before sending to the Gauss frontal solver. The situation would be significantly improved if the program were written as an efficient adaptive system, which could identify regions needed to be

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1-14

13

subdivided from an initial coarse mesh and/or regions that need to be computed with high p levels only.

References [1] V. CASTELLIand J. PIRVlCS,"Review of numerical methods in gas bearing film analysis", Trans. ASME J. Lubr. Technol. 90, pp. 777 792, 1968. [2] J.W. WHITE,"A variable wave speed technique for the steady-state analysis of low-clearance gas bearings", Trans. A S M E J. Lubr. Technol. 99, pp. 339-345, 1977. [3] A SERENYand V. CASTELLI,"Numerical solution of Reynolds equation with slip boundary conditions for cases of large bearing number (A > 300)", Trans. A S M E J. Lubr. Technol. lfll, pp. 64 66, 1979. I-4] J.W. WHITEand A. NIGAM,"A factored implicit scheme for the numerical solution of the Reynolds equation at very low spacing", Trans. A S M E J. Lubr. Technol. 102, pp. 80 85, 1980. [5] C. BAGC!and A.P. SIN6H, "Hydrodynamic lubrication of finite slider bearing: effect of one dimensional film shape, and their computed aided optimum designs", Trans. A S M E J. Lubr. Technol. 105, pp. 48-66, 1983. [6] J.W. WHITE, "Flying characteristics of the zero load slider bearing", Trans. ASME J. Lubr. Technol. 105, pp. 484-490, 1983. [-7] K. KOGURE,S. FUKUI,Y. MITSUYAand R. KANEKO,"Design of negative pressure slider for magnetic recording disks", Trans. A S M E J. Lubr. Technol. 105, pp. 496-502, 1983. [8] J.W. WroTE, "Flying characteristics of the 3370-type slider on a 5 1/4 inch disk Part I: static analysis", in: Tribology and Mechanics of Magnetic Storage Systems, ASLE Publ., SP-16, pp 72-75, 1984. [9] J.W. WHITE, "Flying characteristics of the 3370-type slider on a 5 1/4 inch disk - part II: dynamic analysis", in: Tribology and Mechanics of Magnetic Storage Systems, ASLE Publ. SP-16, pp. 77-84, 1984. [10] K.L. DECKERT,S.A. BOLASNAand H.S. NISHIHIRA,"Statistics of slider bearing flying height", in: Tribology and Mechanics of Magnetic Storage Systems, ASLE Publ., SP-22, pp. 1-5, 1987. [11] S. FUKH and R. KANEKO,"Analysis of ultra-thin gas film lubrication based on linearized boltzmann equation: First report derivation of generalized lubrication equation including thermal creep flow", J. Tribol. 110, pp. 253 262, 1988. [12] P.E. RAAOand J.W. WroTE, "Entrance and stationary roughness effects on the load carrying capacity of a wide wedge gas bearing", J. Tribol. l l l , pp. 41-48, 1989. [13] M.M. REDDI,"Finite element solution of the incompressible lubrication problem", Trans. A S M E J. Lubr. Technol. 91, pp. 524 533, 1969. [14] M.M. REDDI and T.Y. CHu, "Finite element solution of the steady-state compressible lubrication problem", Trans. A S M E J. Lubr. Technol. 92, pp.495-503, 1970. [15] J.F. BOOKERand K.H. HUEBNER,"Application of finite element methods to lubrication: an engineering approach", Trans. A S M E J. Lubr. TechnoL 94, pp. 313 323, 1972. [16] S.M. ROHDE,"Finite element optimization of finite stepped slider bearing profiles", Trans. ASLE 17, pp. 105-110, 1974. [17] S.M. ROHDEand K.P. OH, "Higher order finite element methods for the solution of compressible porous bearing problems", Int. J. Numer. Methods Eng. 9, pp. 903-911, 1975. [-18] Y. MITSUYA,"Molecular mean free path effect in gas lubricated slider bearing (an application of the finite element solution)", Bull. Jpn. Soc. Mech. Eng. 22 (167), pp. 863-870, 1979. [19] C. GARCIA SUAREZ,D.B. BoGy and F.E. TALKE,"Use of an upwind finite element scheme for air-bearing calculation", in: Tribology and Mechanics of Magnetic Storage Systems, ASLE Publ., SP-16, pp. 90 96, 1984. [-20] M. CRONE,M.S. JHON, B. BHUSHANand T.E. KARIS,"Modeling the flying characteristics of a rough magnetic head over a rough rigid-disk surface", J. Tribol. 113, pp. 739-749, 1991. [21] S.H. NGUVEN,"p-version incompressible lubrication finite element analysis of large width bearings", J. Tribol. 113, pp. ll6 ll9, 1991. [22] S.H. NGUYEN, "p-version finite element formulation for analyzing infinitely wide gas bearing problems", Finite Elements in Analysis and Design 9, pp. 141-148, 1991.

14

S.H. Nguyen/Finite Elements in Analysis and Design 16 (1994) 1 14

[23] S.H. NGUYEN,"Solution of Reynolds equation with slip-flow effects by the hierarchical finite element formulation", Tribol. Trans. 34, pp. 565 568, 1991. [24] S.H. NGUYEN,"Use of a p-version finite element formulation for compressible lubrication calculations", Comm. Appl. Numer. Methods 9, pp. 139-145, 1993. [25] W.A. GRoss, Fluid Film Lubrication, Wiley, NewYork, 1962. [26] B.M. IRONS, "A frontal solution program for finite element analysis", Int. J. Numer. Methods Eng. 2, pp. 5--32, 1970. [27] O.C. ZIENKIEWlCZ,The Finite Element Method, McGrawHill, NewYork, 3rd edn., 1983.