SciELO - Scientific Electronic Library Online

 
vol.38 número4Baro-diffusion and thermal-diffusion in a binary fluid mixture confined between two parallel discs in presence of a small axial magnetic fieldReduction of sulfur levels in kerosene by Pseudomonas sp strain in an airlift reactor índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Articulo

Indicadores

  • No hay articulos citadosCitado por SciELO

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google

Bookmark


Latin American applied research

versión impresa ISSN 0327-0793

Lat. Am. appl. res. v.38 n.4 Bahía Blanca oct. 2008

 

Numerical investigation of inertia and shear-thinning effects in axisymmetric flows of carreau fluids by a galerkin least-squares method

R. R. Martins1, F. S. Silveira1, M. L. Martins-Costa2 and S. Frey1

1 Laboratory of Computational and Applied Fluid Mechanics (LAMAC)
Mechanical Engineering Dept. - Universidade Federal do Rio Grande do Sul - Porto Alegre/RS - Brazil
frey@mecanica.ufrgs.br

2 Laboratory of Theoretical and Applied Mechanics (LMTA)
Mechanical Engineering Dept. - Universidade Federal Fluminense - Niterói/RJ, Brazil
laura@mec.uff.br

Abstract — This article presents a finite element simulation of Carreau flows through an abrupt contraction. The employed mechanical model consists in using the Carreau viscosity equation to characterize the shear-thinning fluid behavior, giving rise to a generalization of Navier-Stokes equation containing a non-linear diffusion term. A Galerkin-Least Squares methodology approximates the mechanical model circumventing the Babuška-Brezzi condition, which consists of adding to the classical Galerkin method mesh-dependent residuals, resulting from least squares of the Euler-Lagrange equations. Numerical results for both velocity and pressure fields accounting for shear-thinning and fluid inertia effects have been obtained for an axisymmetric 4:1 sudden contraction with Carreau number ranging from 0 to 100, power-law exponent from 0.2 to 1.0 and Reynolds number from 2 to 100. These results have shown good agreement with the literature.

Keywords — Non-Newtonian Fluids. Carreau Equation. Sudden Contraction Flow. Galerkin-Least Squares Method.

I. INTRODUCTION

Non-Newtonian fluid behavior is present in a wide class of fluids, encompassing suspensions of particles (e.g. blood, paint, ink, and food products), polymer melts and solutions. Many mathematical models describe the non-linear behavior of the material functions obtained experimentally, some of them built within the framework of an axiomatic approach of the non-linear theory of continuous media.

Numerical approximations for shear-thinning flows have deserved a special attention of a large number of researches in non-Newtonian fluid mechanics. Kim et al. (1983) have employed a classical finite element approach to study the roles of fluid inertia and shear-rate dependence of a Carreau viscosity field. The authors have concluded that the effect of increasing either shear-thinning or fluid inertia decreases the upstream vortex size at the sudden contraction flow. Laminar fully-developed flows in an eccentric annular geometry for power-law fluid have been considered by Fang et al. (1999). The employed finite difference scheme gives rise to numerical results for a range of annular radius (0.2≤r≤0.8), eccentricity (0≤ε≤0.8) and power-law coefficient (1≥n≥0.2), with these two latter parameters playing a strong influence on the fluid behavior. Roberts et al. (2001) introduce pseudoplastic empiric models based on Cross and Ellis viscosity functions, taking into account the material yield stresses, used to describe commercial fluids such as drilling mud. Reis Junior and Naccache (2003) simulated, via a finite volume methodology, non-Newtonian fluids through axisymmetric sudden expansion and contraction flows, investigating the influence of rheological parameters on these flows. Power-law fluid flows through a 3:1 sudden expansion were simulated by Manica and De Bortoli (2004), employing a finite difference technique coupled to a Runge-Kutta scheme. The performed simulations have shown bifurcation for values of power-law index between 0 and 2, which occurs either at the beginning of shear thinning or for strong shear thickening. A finite volume method, employing the SIMPLE pressure-correction strategy coupled to the QUICK difference scheme, has been used by Neofytou (2005) to study four non-Newtonian viscous models in a lid-driven cavity flow: power-law, Quemada and modified Bingham and Casson ones.

In the present work the roles of fluid inertia and shear-rate dependency on viscosity for an inelastic non-Newtonian fluid have been analyzed. A Galerkin/least-squares (GLS) methodology was employed to simulate pseudoplastic flows obeying Carreau viscosity function through an axisymmetric 4:1 sudden contraction. This GLS methodology overcomes classical Galerkin shortcomings at high advective flows by adding mesh-dependent terms (functions of the residuals of the Euler-Lagrange equations evaluated elementwise) to classical Galerkin formulation, enhancing its convergence without upsetting its consistency, since the residuals of the Euler-Lagrange equations are satisfied by their exact solutions. Stabilization is important not only to circumvent Babuška-Brezzi condition but also to preserve numerical stability in locally advective dominated flow regions due to the material behavior of shear-thinning liquids.

In all numerical simulations, a structured non-uniform finite element mesh has been employed, consisting of 7,154 equal-order bilinear quadrilateral elements (Q1/Q1), in order to approximate velocity and pressure fields. The Carreau number was investigated from 0 to 100, the shear-thinning coefficient from 0.2 to 1.0 and the Reynolds number from 2 to 100.

A. Some Notation

The problems considered herein are defined on a bounded domain Ω⊂ℜN=2 with a polygonal or polyhedral boundary Γ, formed by the union of Γg - the portion of Γ where Dirichlet conditions are imposed - and Γh - the portion subjected to Neumann boundary conditions. A partition Ωh of into elements is performed in the usual way: no overlapping is allowed between any two elements, the union of all element domains ΩK reproduces and a combination of triangles and quadrilaterals for the two-dimensional case can be accommodated. Quasiuniformity is not assumed (Ciarlet, 1978).

As usual, L2(Ω), L02(Ω),H1(Ω) and H01(Ω) stand for Hilbert and Sobolev functional spaces, respectively, as follows,

(1)

with C0(Ω) representing the space of continuous functions in Ω.

II. MECHANICAL MODELING

Considering an isothermal axisymmetric flow (ui = 0), it suffices to solve mass and linear momentum balance equations, since Cauchy stress tensor is symmetric – automatically satisfying the balance of angular momentum. Let Ω represent an arbitrary fluid volume with boundary Γ, the conservation of mass may be expressed as

(2)

where ur and uz are the non-zero components of the fluid velocity, ρ its mass density and nr and nz are the radial and axial components of the unit outward normal to the surface Γ.

A local form of the continuity equation may be obtained by applying Gauss divergence theorem to Eq. (2) and considering it valid for any arbitrary volume Ω.

(3)

Restricting further the discussion to incompressible fluids, the equation of continuity, Eq. (2) , reduces to

(4)

From the linear momentum conservation (Bird et al., 1987), the balance of linear momentum may be mathematically expressed by

(5)

in which, considering axisymmetric flows, Trr, Trz, Tzr, Tθθ and Tzz, are the non-zero components of Cauchy stress tensor and ρfr and ρfz represent the radial and axial components of the body force per unit volume, respectively.

Assuming the stress tensor T decomposition in a spherical and a deviatoric portion T= - p I+τ – with p representing a non thermodynamic pressure, I the unit tensor, and τ the deviatoric stress tensor – the local Eulerian description of the motion equation, Eq. (5), may be obtained by applying Gauss theorem and the continuity equation (Eq. (2)-(4)) and by considering it valid it for any arbitrary fluid volume Ω

(6)

In inelastic non-Newtonian fluid flows, the deviatoric stress tensor may be related to fluid kinematics by a generalization of Newton law for viscosity. A non-linear dependence of the deviatoric stress on the rate of strain tensor may be introduced by considering a generalized Newtonian constitutive law,

(7)

in which D is the symmetric portion of the velocity gradient – denoted by rate of strain tensor – and the shear rate viscosity. For an axisymmetric flow, the components of the deviatoric stress tensor may be expressed as

(8)

In the present work the viscosity function is given by Carreau constitutive equation (Bird et al., 1987)

(9)

with η0 and η being asymptotic values of fluid viscosity at zero and infinite shear rates, respectively, λ being a characteristic time equal to the reciprocal of the shear rate at which shear thinning begins and, finally, (n-1) represents the power-law slope of the logarithmic of viscosity function, . The shear rate scalar, , represents the Frobenius norm of tensor D, a mathematical measure for the shear rate when simple shear flow is assumed – namely: .

III. FINITE ELEMENT APPROXIMATION

Combining continuity and motion equations, Eq. (4) and (6), with Carreau law, Eq. (9), the following boundary value problem may be stated for steady-state laminar flow,

(10)

in which Ω represents the internal domain, Γg and Γh are the portions of the boundary Γ on which Dirichlet and Neumann boundary conditions were imposed, respectively, and are the non-zero components of the stress vector on .

A. Galerkin Least-Squares formulation

The finite element approximation of Eq. (10) was built in by employing the usual fluid dynamics subspaces for velocity (Vh) and pressure (Ph) (Ciarlet, 1978),

(11)
(12)
(13)

with Rm (m=k,l) denoting the polynomial spaces of degree m, Ωh a usual finite element partitioning (Ciarlet, 1978) and ΩK the domain of the K-finite element of Ωh.

Based on the above definitions of velocity and pressure subspaces, Eq. (11)-(13), a Galerkin least-squares formulation for the problem presented in Eq. may be stated as: Find the pair (uh, ph) ∈ VhgxPh, such that, for all (vh, qh) ∈ VhgxPh :

(14)

in which the stability parameter τ is defined as

(15)

with denoting the p-norm on ℜN and the parameter mk being determined from the error analysis introduced by Franca and Frey (1992).

Remark: When the parameter τ is equal to zero in the GLS formulation defined in Eq. (14)-(15), the classical Galerkin formulation is recovered.

B. Nonlinear scheme

When the shape functions for uh, ph, vh and qh are introduced in the GLS formulation, Eq. (14)-(15), an algebraic residual set of equation is yielded in the form:

R(U) = 0 , (16)

where U is the vector of degrees of freedom of uh and ph and R(U) is given by the set of matrices

(17)

In Eq. (17), matrices [K] and [G] are, respectively, originated by the diffusive and pressure terms of Eq. (17), and N(u), GT, F, respectively by the advective, incompressibility and body force ones. (The Nτ, Kτ, Gτ, and Fτ in Eq. (17) are generated by the least-squares terms of Eq. (14)-(15)) To solve system Eq. (14)-(15) a quasi-Newton method (see for instance, Dalquist and Bjorck, 1969) has been implemented where the jacobian matrix was updated only at each two or three iterations.

ALGORITHM:

I. Estimate vector U0 and set the number of iterations (m) to update jacobian matrix J(U).
II. Set k = 0, j = 0, ε = 10-7.
III. If k -int(k / m)*k =0, then j = k.
IV. Solve for incremental vector ak+1:

(18)

where R(U) is given by Eq. (18) and J(U) defined by

(19)

V. Compute vector Uk+1:

(20)

VI. Compute , then do k = k+1 and go to step III; otherwise, store solution Uk+1 and exit from the algorithm.

Remark: As initial solution estimates, null velocity and pressure fields were employed in the ALGORITHM. In order to improve algorithm convergence, a continuation method (Dalquist and Bjorck, 1969) was implemented on advective terms of Eq. (16)-(17).

IV. NUMERICAL RESULTS

In this section results employing the Galerkin/Least-Squares (GLS) finite element approximation for the problem defined by Eq. (10) are presented and discussed. The numerical computations consider the flow of a Carreau fluid – Eq. (9) – through an axisymmetric abrupt 4:1 contraction. The computational results have been obtained by employing a finite element code developed at Laboratory of Applied and Computational Fluid Mechanics (LAMAC) at Federal University of Rio Grande do Sul.

A. Code validation

The computational implementation of the GLS formulation has been validated by the benchmark of the lid-driven cavity problem for Newtonian fluids subjected to high advection flows. The cavity layout consists of a bi-unity cavity with a moving lid with impermeability and no-slip boundary condition at its walls, except for the lid at the superior edge which moves with a steady horizontal prescribed velocity – see, for instance Neofytou (2005), for problem statement details. A uniform finite element partitioning of the computational domain in 130x130 quadrilateral bilinear elements (Q1/Q1) has been considered in the usual way (Ciarlet, 1978), generating 17,161 nodal points. Three distinct Reynolds number were investigated, namely Re=1 (for simulating creeping flows), 400 and 1000.

Figures 1 and 2 show the comparison between GLS approximations obtained in this work and those by Ghia et al. (1982) including inertia effects — namely Re=400 e Re=1000 — as well as with results from Jurjevic (1999) for creeping flow (Re=1). The horizontal   at x=0.5 (Fig. 1) and vertical at y=0.5 (Fig. 2) profiles have shown good agreement with both articles for the investigated Reynolds numbers.

B. Carreau flow through an axisymmetric abrupt contraction

This section presents results obtained by using the GLS approximation defined in Eq. (14)-(15) for a Carreau fluid (Eq. (9)), through an axisymmetric 4:1 sudden contraction. The geometry was built in by joining a cylinder of radius R1 and length L1 with another one of smaller radius R2 and length L2. The geometric configuration is mathematically described by the aspect ratio βR1/R2. All simulated cases have employed β=4, as shown in Fig. 3(a).


Figure 1. Lid-driven cavity flow — horizontal velocity profiles at x=0.5: (a) Re=1; (b) Re=400; (c) Re=1000.

The imposed boundary conditions were: no-slip and impermeability at walls, symmetry velocity condition at the centerline; at the inlet a flat velocity profile was imposed while at the outlet free traction condition was adopted. Besides, the pressure is fixed at the outlet in one point: pref=0.

After a mesh refining process ensuring mesh independence, the computational domain has been partitioned in the usual way (Ciarlet, 1978) into 7154 bilinear elements (Q1/Q1)- as depicted in Fig. 3(b).

Reynolds number for this flow was defined considering the smaller radius cylinder R2, its average velocity V2 and the zero-shear-rate viscosity η0:

(21)

Introducing the Carreau number as

(22)

with the parameter λ defined as previously, and setting the infinite-shear-rate viscosity η=0, the Carreau viscosity equation (Eq. (9)) may be redefined as

(23)

According to Eq. (23), as illustrated in Fig. 4, with the shear stress given by , it may be seen that a Newtonian behavior was recovered for Cu=0. Besides when Cu tended to zero the flow curves were shifted to the left.


Figure 2. Lid-driven cavity flow — vertical velocity profiles at y=0.5: (a) Re=1; (b) Re=400; (c) Re=1000.



Figure 3. Axisymmetric abrupt contraction flow: (a) problem statement; (b) detail of refined mesh at the contraction region.


Figure 4: Flow curves for n=0.5 and distinct Cu values: (a) ; (b)

Results have been obtained considering Reynolds number values varying from Re=2 to Re=100 the former simulating negligible inertia flows — and the shear-thinning has been investigated by considering distinct values of Carreau number — ranging from Cu=0 to Cu=100 — and by varying the power-law exponent from n=0.2 to n=1.0.

In order to analyze the influence of Carreau number at the flow dynamics, at Figs. 5 to 7 inertial effects have been neglected (Re=2) and the power-law coefficient was fixed as n=0.2. Figure 5 presents Carreau viscosity contours for Cu from 0 to 100. As expected from Eq. (23), the viscosity function for Cu=0 (Newtonian case) is a constant surface (Fig. 5(a)). Regarding Fig. 5(b), it may be observed that the pseudoplastic effects are restricted to regions of higher shear rates — chiefly near the wall of the smaller tube — for small Carreau numbers (Cu=10). As Cu increases the shear-thinning effects are also present at regions of lower shear-rates of the larger tube, as depicted in Figure 5(c)-(d), for Cu=50 and Cu=100, respectively.

Figure 6 presents pressure contours for Re=2 and n=0.2. From all figures it may be verified that an increase at Carreau number causes a monotonic decrease at pressure drop along the flow. This decrease is particularly accentuated for Cu between 0 (Fig. 6(a) for Cu=0) and 10 (Fig. 6(b) for Cu=10). Besides, for the highest values of Cu (Fig. 6(c) for Cu=50 and (Fig. 6(d) for Cu=100), the pressure drop at the smaller tube is attenuated by shear-thinning effects. These comments could be anticipated by observing the progressive decrease of fluid viscosity shown in Fig. 5. The increase of Carreau number diminishes the fluid's resistance to flow and, consequently, the pressure drop.

Figure 7 presents a detail of the flow streamlines at the sudden contraction region, once again varying Cu from 0 to 100 and considering Re=2 and n=0.2. For Cu=0 — Newtonian case — a well-defined vortex is present at the contraction corner. This vortex tends to be collapsed as the shear-thinning effect becomes stronger, for Cu=50 and Cu=100 (Figs. 7(c)-7(d)). As the shear thinning effect increases with Carreau number growth, the decreasing viscosity forces the flow at the larger tube to be locally advective dominated, giving rise to the vortex collapse.





Figure 5: Viscosity contours for n=2 and Re=2: (a) Cu=0; (b) Cu=10; (c) Cu=50 (d) Cu=100.





Figure 6. Pressure contours for n=2 and Re=2: (a) Cu=0, (b) Cu=10, (c) Cu=50 and (d) Cu=100.





Figure 7: A detail of flow streamlines at the sudden contraction region for n=2 and Re=2: (a) Cu=0, (b) Cu=10, (c) Cu=50 and (d) Cu=100.

Figures 8 and 9 have also been obtained by neglecting inertial effects (Re=2). Figure 8 shows the influence of the power-law coefficient (ranging from n=0.2 to n=1.0) on the axial velocity profile at the contraction plane considering two distinct values of Carreau number — namely Cu=10 in Fig. 8(a) and Cu=100 in Fig. 8(b).

Comparing both Figs. 8(a) and 8(b), it may be noticed an increase of pseudoplasticity due to Carreau number augmentation, in which the power law index decrease causes the shear-thinning intensification. Figure 8 shows the evolution from the parabolic Newtonian profile (n=1) to an almost flat profile (n=0.2) — in both figures. This almost flat profile with very thin boundary layers near the walls — presenting severe velocity gradients, is characteristic of very strong shear-thinning, obtained by combining n=0.2 and Cu=100. The deviation from the classical Newtonian pattern (n=1) at the contraction plane results from the progressive decrease of viscosity related to the decrease of the power-law coefficient, being intensified by Carreau number augmentation. The concavity at the boundary layer edge for the strongest shear-thinning case (n=0.2 and Cu=100 — in Fig. 8(b)) may be eliminated by employing shock-capture strategies (Galeão and Carmo, 1988).

Further, Fig. 8(b) presents a comparison between results obtained in the present work for n=0.4 and those by Kim et al. (1983). These authors have employed a classical Galerkin finite element approximation with a very coarse mesh, consisting of 116 Q2/Q1 elements. As observed from this figure both results show an excellent agreement.

Figure 9 shows the influence of the power-law coefficient (ranging from n=0.2 to n=1.0) on the axial velocity profiles along symmetry axis, for Re=2 and Cu=100. At the larger tube (a low shear rate region) the velocity profiles are almost independent from the power-law coefficient; otherwise, at the smaller one, the velocity profile suffers a strong influence of the power-law coefficient. As may be noticed, at the latter region the velocity values decrease as n decreases. A comparison between the present work results for n=0.4 and those by Kim et al. (1983) is also presented in Fig. 9. Both results show a very good agreement, but in the vicinity of the contraction plane.



Figure 8. Axial velocity profiles at contraction plane for Re=2 and distinct values of n: (a) Cu=10 and (b) Cu=100.


Figure 9. Axial velocity profiles along symmetry axis for Re=2, distinct power-law coefficients and Cu=100.

Figures 10 and 11 investigate the influence of the inertia effects on the flow dynamics. Figure 10 presents radial velocity at the contraction plane for Reynolds numbers varying from 2 to 100. A fully-developed parabolic profile at the contraction plane is obtained for Re=2. The Reynolds number growth gives rise to quasi-uniform velocity profiles at the symmetry line subject to severe boundary layers. The most advective profile (Re=100) presents an overshoot at boundary layer edge, which may be eliminated by employing a shock-capture strategy (Galeão and Carmo, 1988).


Figure 10. Radial velocity profiles for disinct Reynolds numbers and Cu=0 at contraction plane.


Figure 11. Axial velocity profiles along symmetry axis for disinct Reynolds numbers, considering Cu=100 and n=0.2.

Figure 11 illustrates the axial velocity profiles along the symmetry axis, for the most shear-thinning case (Cu=100 and n=0.2) with Reynolds number varying from 2 to 55. This figure reveals an almost negligeable influence of Reynolds number on the velocity profiles. In the larger tube, even for the least Reynolds value (Re=2), the flow is advective-dominated due to the high value of the shear-thinning coefficient (n=0.2) — as indicated by the overlaping of all curves. Yet in the smaller tube one may notice a differentiation of the curves with the increase of Reynolds number. As Reynolds grows one may observe only a small augmentation of the centerline velocites due to the very low value of the shear-thinning coefficient of the fluid (n=0.2).

V. FINAL REMARKS

A Galerkin least-squares finite element formulation has been employed to approximate an isochoric Carreau flow through an axisymmetric sudden contraction, with the results being successfully compared to those by Kim et al. (1983). It was verified that the increase of shear-thinning effect causes a strong influence on the flow dynamics. The velocity profiles at the contraction plane become high advective-dominated with the increase of the shear-thinning coefficient and the increase of the Carreau number, with the presence of a strong boundary layer near the pipe wall. As shear-thinning increases, three distinct behaviors were verified in this article: the velocity profiles at the contraction plane become flatter; the vortex size at the contraction corner decreases until its collapse and the pressure drop along the flow suffers a strong reduction.

ACKNOWLEDGMENTS
The authors M. L. Martins Costa and S. Frey gratefully acknowledge the financial support provided by the Brazilian agency CNPq. The author F.S. Silveira acknowledges the scholarship provided by CNPq.

REFERENCES
1. Bird, R.B., R.C. Armstrong and O. Hassager, Dynamics of polymeric liquids, John Wiley & Sons (1987).         [ Links ]
2. Ciarlet, P.G., The finite element method for elliptic problems, North-Holland, Amsterdam (1978).         [ Links ]
3. Dahlquist, G. and A. Bjorck, Numerical Methods, Prentice-Hall, Englewood Cliffs, NJ (1969).         [ Links ]
4. Fang, P., R.M. Manglik and M.A. Jog, "Characteristics of laminar viscous shear-thinning fluid flows in eccentric annular channels", J. Non-Newt. Fluid Mech., 84, 1-17 (1999).         [ Links ]
5. Franca, L.P., and S. Frey, "Stabilized finite element methods: II. The incompressible Navier-Stokes equations", Comput. Methods Appl. Mech. Engrg., 99, 209-233 (1992).         [ Links ]
6. Galeão, A.C., and E.G.D. do Carmo, "A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems", Comput. Method Appl. Mech. Engrg., 68, 83-95 (1988).         [ Links ]
7. Ghia, U., K.N. Ghia and C.T. Shin, "Hi-Re solution for incompressible flow using the Navier-Stokes equationsand the multigrid method", J. Comput Physics., 48, 387-411 (1982).         [ Links ]
8. Jurjevic, R., "Modelling of two-dimensional laminar flow using finite element method", Int. Journal for Numerical ethods in Fluids, 31, 601-626 (1999).         [ Links ]
9. Kim, M.E., R.A. Brown and R.C. Armstrong, "The roles of inertia and shear-thinning in flow of an inelastic liquid through an axisymmetric sudden contraction", J. Non-Newt. Fluid Mech., 13, 341-363 (1983).         [ Links ]
10. Manica, R., and A.L. De Bortoli, "Simulation of sudden expnsion flows for power-law fluids", J. Non-Newt. Fluid Mech., 121, 35-40 (2004).         [ Links ]
11. Neofytou, P., "A 3rd order upwind finite volume method for generalized Newtonian fluid flows", J. Non-Newt. Fluid Mech., 36, 664-680 (2005).         [ Links ]
12. Reis Junior, L.A., and M.F. Naccache, "Analysis of non-Newtonian flows trough contractions and expansions", XVII Congresso Brasileiro de Engenharia Mecanica, Sao Paulo (2003)         [ Links ]
13. Roberts, G.P., H.A. Barnes and P. Carew, "Modelling the flow behaviour of very shear- thinning liquids", Chem. Engn. Sciences, 56, 5617-5623 (2001).         [ Links ]
14. Zinani, F. and S. Frey, "Galerkin Least-Squares Finite Element Approximations for Isochoric Flows of Viscoplastic Liquids", Journal of Fluids Engn.— Transactions of the ASME, 128, 853-863 (2006).
        [ Links ]

Received: November 23, 2007.
Accepted: Jannary 31, 2008.
Recommended by Subject Editor: Eduardo Dvorkin.