Flows of Bingham Materials Through Ideal Porous Media: an Experimental and Theoretical Study

Flows of Bingham Materials Through Ideal Porous Media: an Experimental and Theoretical Study

P. R. S. Mendes

M. F. Naccache

C. V. M. Braga

A. O. Nieckele

F. S. Ribeiro

Department of Mechanical Engineering

Pontifícia Universidade Católica do Rio de Janeiro

22453-900 Rio de Janeiro, RJ. Brazil

The flow of Bingham liquids through porous media has been studied. Experiments have been performed to determine the flow rate / pressure drop relationship for the flow of a grease of Binghamian rheological behavior through an array of rods of circular cross section. The yield stress and plastic viscosity of the grease have been determined with the aid of a controlled stress rotational rheometer. To investigate a wider range of the flow parameters, the mass and momentum conservation equations have been solved numerically, in conjunction with the generalized Newtonian constitutive law and the bi-viscosity model. The finite volume method has been employed to obtain the numerical solution. These numerical results also yielded a flow rate / pressure drop relationship, which is in very good agreement with the experimental results. A capillaric theory has been developed to determine an analytical relationship between the flow rate and pressure drop for flows of Bingham liquids through porous media. It is shown that the predictions of this theory are in good agreement with the experimental and numerical results.

Keywords: Porous media, Bingham material, viscoplastic

Introduction

When a petroleum well is drilled, the pressure in the porous rock near the well decreases, altering the thermodynamic state of the petroleum. Light hydrocarbons volatilize, which causes the temperature and solubility of the heavy hydrocarbons to decrease. In the case of waxy crudes, paraffins tend to precipitate off the oil.

The presence of sub-micron precipitated paraffin crystals cause the oil to become a colloidal dispersion, whose viscosity is high, and whose mechanical behavior is often non-Newtonian. Typically the oil acquires a non-zero yield stress and a shear thinning viscosity. Elasticity is also possible, although less common.

Flows of non-Newtonian fluids through porous media do not obey Darcy's law. When the fluid is elastic, the converging-diverging passages found in porous media, which impose extensional kinematics, cause normal stresses to reach high values, becoming the main source of flow resistance (Durst et al, 1987; Souza Mendes and Naccache, 1994a).

Purely viscous non-Newtonian fluids are well represented by the generalized Newtonian constitutive relation. If the fluid has a yield stress, often the Bingham plastic viscosity function is representative of rheological data. For oils with this type of rheological behavior, Wu et al. (1990) proposed a seepage velocity/pressure gradient relation where a "yield pressure gradient" G appears, which is a function of both the fluid and the porous medium:

where u is the seepage velocity, m_p the plastic viscosity, K_N the Newtonian permeability and p the pressure. This equation has been proven useful in reproducing experimental data, via least-squares fits to the data for the determination of G, for each combination of porous medium and fluid. However, it will be shown that Wu's equation has a qualitative behavior that is not in agreement with the one observed experimentally for low flow rates.

Nomenclature

b = one of the cross section dimensions (Fig.3), m

D = rod diameter (Fig. 1), m

D* = dimensionless rod diameter

D = rate-of-deformation tensor, s^-1

G, G' = yield pressure gradients, Pa/m

H = one of the cross section dimensions (Fig. 3), m

L = row spacing (Fig. 1), m

L* = dimensionless row spacing

L_t = tube length, m

K_N = Newtonian permeability, m²

p = pressure, Pa

Q = flow rate, m³/s

R = tube radius, m

Re = Reynolds number, dimensionless

S = column spacing, m

u = average seepage velocity, m/s

u = local seepage velocity vector, m/s

v = velocity vector field, m/s

= average velocity, m/s

W = spacing between pressure taps (Fig. 3), m

x = cartesian coordinate in the main flow direction (Fig.1), m

y = cartesian coordinate normal to the main flow direction (Fig.1), m

Y = yield number, dimensionless

Greek Symbols

= deformation rate, s^-1

= small deformation rate

of bi-viscosity model, s^-1

D p = pressure drop along tube, Pa

h() = viscosity function, Pa.s

f = porosity, dimensionless

m_p = plastic viscosity, Pa.s

m_¥ = large viscosity of bi- viscosity model, Pa.s

t = extra-stress tensor, Pa

r = mass density, kg/m³

t = von Mises extra-stress, Pa

t ₀ = yield stress, Pa

t _R = wall shear stress, Pa

The present work reports a performance evaluation of an equation firstly proposed by Souza Mendes and Naccache (1994b) for flows of Bingham liquids through porous media. The goal of this evaluation is to assess its capability of describing the influence of rheological parameters.

To this end, the flow through an ideal porous medium is studied numerically and experimentally. The geometry of this ideal porous medium differs significantly from the one assumed in the theory that led to the proposed equation. The flow rate/pressure drop relationship obtained from these studies is then compared with the predictions of the proposed equation.

Theory

The simple theory presented in this section is aimed at obtaining a relationship between the flow rate and the pressure drop for the flow of a Bingham liquid through porous media (Souza Mendes and Naccache, 1994b).

The fluid considered here has a Binghamian behavior, which is accounted for through the use of the Generalized Newtonian Liquid (GNL) model:

where t is the extra-stress, D º [grad v + (grad v) ^T] / 2 the rate of deformation, and is a scalar measure of its intensity.

For the Bingham liquid, the viscosity function is

where t₀ is the fluid yield stress.

Starting with a differential momentum balance for the fully developed flow through a tube, it is possible to obtain an expression for the flow rate as a function of the pressure drop for this non-Newtonian fluid:

where t _R=R D p/(2L_t ) is the shear stress at the tube wall, R the tube radius, and Dp the pressure drop along the tube length L_t. Equation (4) is the so called Buckingham's equation, which applies only when t_R >t ₀.

For high flow rates, t_R is much larger than t₀, and the term (1/3) (t₀/t_R)⁴ is much smaller than the other terms in Eq. (4), and often can be neglected. The error in omitting this term is less than 1.8% when t₀/t_R < 0.4. However, for low flow rates, the fourth power term is of the same order of magnitude of the others, and cannot be neglected.

Classic capillaric theories assume that the porous medium is composed of several parallel capillary tubes (Scheidegger, 1974). Then Eq. (4) is applied to each of the tubes, and the average velocity through the medium is taken as the average velocity through each tube times the porosity. Equation (4) can be written in the following form after neglecting the above mentioned term and applying the ideas of the capillaric theories:

Equation (5) can be applied locally if it is re-written in vector form as

It is worth noting that this reduces to Eq. (1) if

and

Therefore, because the fourth power term is neglected in its derivation, Eq. (1) is expected to perform well for large flow rates only. This can be a major drawback. For example, in flows of petroleum through reservoirs toward the well, typical flow velocities are of the order of a few centimeters per day.

If the fourth power term in Eq. (4) is preserved, similar arguments and some more involved algebra lead to the following expression:

where G' is given by

In Eqs. (8), (7) and (10) the tube radius R should be regarded as a characteristic length of the pore.

Unfortunately, pore geometries are very complex and non-uniform. Then, the assumption that just one parameter (like the tube radius R) represents the morphology of typical porous media is seldom enough. Another assumption about the pore channel geometry, which introduces more geometrical parameters, is also considered in Souza Mendes and Naccache (1994b), where a convergent-divergent channel is considered and compared with the straight tube geometry. It is shown that a geometric parameter related to the ratio of throat to cavity diameters may allow a better representation of the effects of pore morphology on the permeability.

Numerical Study

Previous Work

The geometry of the ideal porous medium employed in the numerical and experimental studies is described with the aid of Fig. (1). This geometry is common in heat exchanger, nuclear reactor, and other engineering applications.

If the fluid flowing externally to the rods is Newtonian, pressure drop and heat transfer information abounds in the literature, and is reviewed by Zukauskas (1972). For non-Newtonian liquids, some work is found for flow of pseudoplastic power-law liquids around a cylinder (e.g., Shah et al., 1962, Mizushina et al.,} 1978, and Mizushina & Usui, 1978). For flow past rod bundles of viscoelastic polymeric solutions, pressure drop data are reported by Chmielewski and Jayaraman (1992). However, no information is available in the open literature for Bingham liquids, as far as the authors know.

Analysis

In the present analysis, the liquid is assumed to be incompressible, so that the equation of mass conservation reduces to

where v is the velocity vector field.

The momentum equation is, for steady flow and negligible external forces,

where t the extra-stress and r the mass density.

The extra-stress t is a function of the flow kinematics. For the Generalized Newtonian Liquid model (GNL), it is given by Eq. (2). For a Bingham liquid, the viscosity function, h(), is given by Eq. (3).

The viscosity function as given by Eq. (3) is not simple to model numerically, and the usual approach is to replace it by another viscosity function, the so called bi-viscosity model (Beverly and Tanner, 1992):

where

Equations (12), (2) and (13), plus the continuity equation (Eq. (11)) complete the problem formulation, together with appropriate periodically-developed (cyclic) boundary conditions, i.e., same inlet and outlet velocity profiles. Assuming two-dimensional flow, these equations are written in Cartesian coordinates and then integrated numerically in the domain shown in Fig. 1.

The physical dimensionless parameters that govern the problem are the Reynolds number, Re, and the yield number, Y, defined respectively as

In addition, there are two geometrical parameters, the ratios

where D is the rod diameter, L the row spacing, and S the column spacing (Fig. 1).

Because typical flows of petroleum in reservoirs are rather slow, in the present paper attention is focused toward situations of negligible inertial forces (Re® 0), and therefore the effect of the Reynolds number is not examined.

Numerical Method

The conservation equations were discretized by the finite volume described by Patankar (1980), by using an non-orthogonal curvilinear system of coordinates which adapts to the boundaries of the domain. Staggered velocity components were employed to avoid unrealistic pressure field and the contra-variant velocity component was selected as the dependant variable in the momentum conservation equations (Pires and Nieckele, 1994a, 1994b).

The contra-variant velocity component conservation (U_x, V_h ) equations are obtained from the cartesian velocity component (u, v) equations by

The resulting discretized equation for U_x , at the east face, as shown in Fig. 2, can be written as

where the neighboring coefficients a_nb are determined based on the Power-law approximation. A_e and A_e' are geometric parameters. The latter is zero for orthogonal mesh. b_{U x ,e} is due to the curvature of the coordinate system, ant it is given by

For example, is the velocity component at the face w but parallel to the contravariant component at face e

The cartesian velocity components needed for the above equation are determined from

where Ja is the Jacobian.

The source term is obtained from the source term of the momentum equations for the cartesian velocity component b_u and b_v, as

The pressure-velocity coupling was solved by an algorithm based on the SIMPLEC (Van Doormaal and Raithby, 1984). The resulting algebraic system was solved by the TDMA line-by-line algorithm (Patankar, 1980) with the block correction algorithm (Settari and Aziz, 1973) to increase convergence rate.

A transfinite interpolation scheme was employed to generate a mesh of 40 x 40 control volumes in the computational domain.

Further details of this study, and results for a wide range of the governing parameters, are given in Nieckele et al. (1995).

Experimental Study

Before Eq. (9) can be employed in engineering design, some experimental evidence of its reliability is desirable. In order to assess its description capability regarding the fluid rheology, perhaps the first step is to check its performance against experimental data pertaining to an ideal porous medium. In this manner, influences other than the ones related to rheological parameters are kept under full control.

Due to the previously mentioned scarcity in the literature of experimental data related to flows of Bingham liquids through porous media, experiments have been performed to generate some data to be used in the desired assessment. These experiments are now described.

Test Fluid

A commercial grease has been employed in the experiments. It has been rheologically characterized in the facilities of the Laboratory for Characterization of Fluids (LCF) at PUC-Rio, and the following was observed:

1. Its density is r =884,kg/m³.

2. Its yield stress is t₀ = 450 Pa. The measuring method of this quantity was a stress sweep, with the aid of a rotational rheometer.

3. Its viscosity is rather dependent on the shearing time (thixotropy). However, a representative value of its plastic viscosity for the experiments was m_p=110 Pa.s. The rotational rheometer was also employed in these viscosity measurements.

Test Section

The test section is described with the aid of Fig. 3.

Inside a duct of square cross section (100 x 100 mm) and total length C = 700 mm, twenty seven cylindrical rods (D = 9.3 mm) are assembled in the staggered arrangement shown in Fig. 3. This rod assembly is located at half-length of the square duct, to eliminate entrance and exit effects on the measurements. To reduce the influence of the side walls, three semi-cylindrical rods are attached to each of two of the side walls, following the same staggered pattern. The longitudinal pitch is S=20 mm, while the transversal one is L=10 mm (Fig. 3). There are two pressure taps, one located just upstream of the rod assembly and other just downstream.

The Experimental Apparatus

Upstream and downstream of the square duct which contains the test section, there are two identical cylinder/piston assemblies (Fig. 4). Prior to the runs, the grease is carefully loaded into the system, in such a way as to eliminate the presence of air bubbles.

During a run, the upstream piston is pushed with the aid of compressed air. The air pressure is kept at a selected level with the aid of an appropriate constant pressure valve. Because the upstream and downstream cylinder/piston assemblies are identical, the total resistance to motion is invariant during a given run, which implies a constant flow rate through the test section.

The grease flow rate is easily and accurately determined by measuring the upstream piston displacement rate. This is done with the aid of stopwatch and reference notches, which were made on the surface of the piston rod.

The pressure drop across the test section is determined with a differential manometer, which full scale is 2 kgf/cm².

This assembly yielded good repeatability and low scattering of the data, as it will be seen later.

Results and Discussion

It is important to emphasize that the theory used to obtain Eq. (9) employs simplifying assumptions that might imply limitations in its range of applicability. For instance, the geometry of a porous rock is by far more complex than the assumed parallel straight tubes. Also, molecular flow effects, which might be present in flows through porous rocks, are also neglected in the theory when the no-slip condition is assumed at solid-liquid interfaces. It should be emphasized also that both the experiments and the theory do not consider effects inherent to micro-size pores.

Because of these possible limitations of the theory, it is interesting to start the validation process by comparing its predictions with results that are free from influences of physical mechanisms other than the ones taken into account by the theory. Results obtained in the present research for the porous medium of idealized geometry fall into this category. Therefore, if good agreement is observed, then it can be concluded that the theory is capable of well representing the main effect of interest here, viz., the dependence of the u x ê grad p ê relationship on the fluid rheological parameters (t₀ and r_p).

In addition, the uncertainty of results of experiments with ideal porous media like the one performed in the present research is typically much lower than the uncertainty related to experiments using natural porous rocks. Therefore, comparisons with theoretical predictions tend to be more conclusive.

Furthermore, the geometrical simplicity of the ideal porous medium chosen (Fig. 1) allows numerical solutions of the conservation equations that govern the flow past the rod bundle. Because the phenomena taken into account in these governing equations are known, good agreement between numerical and experimental results indicates that the meaningful experiment parameters were kept under control in the course of each data run. Moreover, the availability of numerical results gives another tool for theory evaluation.

Comparisons: Qualitative Behavior

Figure 5 presents numerical and experimental results for the flow of the grease in the geometrical configuration shown in Fig. 3. The predictions of Eqs. (9) and (6) for this situation are also shown. In this figure, the open circles represent the experimental data obtained with the apparatus shown in Figs. 3 and 4, while the black squares pertain to the numerical results, obtained as explained before.

The curve shown in Fig. 5 for the Eq. (9) was obtained by a procedure of least-squares fitting to the experimental data which assumes the functional relationship determined by the theory, viz., Eq. (9) (also shown in the insert of Fig. 5). In this fitting procedure, R and G' are the parameters to be determined. Fig. 5 also shows the curve for the Eq. (6), with the same values for the parameters R and G'.

All the other parameters appearing in Eq. (9), and the ones needed to obtain the numerical results, were taken to be equal to the ones that actually occurred in the experiments, namely: 0.83, m_P=110 Pa.s, t₀ = 450 Pa, r = 884 kg/m³.

It is observed that the qualitative behaviors of the numerical, experimental and theoretical predictions are the same:

1. At low-pressure gradients there is no flow.

2. At a certain critical value of the pressure gradient, the onset of motion is observed.

3. In the range of pressure gradients just beyond this critical value, the velocity/pressure gradient dependence is clearly non-linear. It is important to notice that Eq. (9) doesn't predict this behavior.

4. For larger pressure gradients there is an asymptotic trend toward a linear behavior.

The trends described above are interpreted as follows:

• For low pressure gradients, the stresses that occur are lower than the yield stress of the fluid everywhere in the flow field, and hence there is no motion.

• For pressure gradients above the critical value, the fluid stresses exceed the yield stress, and motion takes place.

• In the range of pressure gradients just beyond the observed threshold, the fluid viscosity decreases fast with the deformation rate (and hence with the pressure gradient), because in this range it is highly influenced by the yield stress, t₀ (see Eq. (3)). Thus, in this range the permeability ¾ which is the slope of the curve in Fig. 5 ¾ increases.

• In the range of higher pressure gradients, the deformation rates are also higher, and, consequently, the viscosity asymptotically approaches a constant value, namely, m_P. Thus, the permeability also approaches a constant value asymptotically, as seen in Fig. 5.

Comparisons: Quantitative Evaluation

In the insert of Fig. 5, R and G' values are given as determined by the least-squares procedure mentioned above, namely, R = 16.5 mm and G' = 3.76397 ' 10⁵ Pa.

It is interesting to compare the above R value with a characteristic length related to the flow passages found in the rod-bundle porous medium. A possible choice for this characteristic length is the rod spacing in the narrowest cross section. It can be observed in Fig. 3 that, if this characteristic length is chosen, for this case it would be equal to 10.7 mm. Therefore, the R value obtained via least squares is of the same order of magnitude as the characteristic length of the porous medium employed in the experiments differing by a factor of 1.5. This discrepancy is acceptable if it is recalled that the geometry assumed in the theory is completely different from the one related to the data.

Another interesting comparison consists in employing Eq. (10) for evaluating G', using the values for t₀ and R relative to the experiments:

If this value is compared with the one obtained via least squares, it can be seen that they differ by a factor of 4.5. Again, the discrepancy is large, but the order of magnitude provided by the theory (Eq. 10) is correct.

These discrepancies are surely due to the dramatic difference between the morphology assumed in the theory and the one employed in the experiments and computation. It seems reasonable to expect that the qualitative behavior observed and discussed above be independent of the morphology, and thus, it should also be observed in more complex morphologies like the ones found in petroleum reservoirs. Quantitatively, however, the results should be a function of the porous medium morphology.

Conclusions

The comparisons discussed above suggest that:

• The proposed theory describes correctly the effect of a Binghamian rheology in flows through porous media, in the same way as Darcy's law describes correctly the flow of Newtonian fluids through porous media. However, more data are needed to fully validate the theory.

• The qualitative behavior observed for Bingham liquids should not depend on the porous medium morphology, but, quantitatively, the morphological characteristics should matter.

• In the same way as Darcy's law has its limitations, since it groups two many important effects into one parameter only (the permeability), the same is true for Eq. (9).

• The theory should be useful in fittings to experimental data, as well as in oil reservoir simulations. It is worth noting that, for flows at low rates of Bingham liquids through porous media, Eq. (1) will not yield good fittings because it is a linear relation. Equation (9) has a correct qualitative behavior at all flow rates.

Acknowledgements

Financial support for the present research was provided by Petrobras S. A., CNPq and MCT.

Manuscript received: September, 2000. Technical Editor: Átila P. Silva Freire.

• Beverly, C. R., and Tanner, R. I., 1992, "Numerical analysis of three-dimensional Bingham plastic flow",' J. Non-Newtonian Fluid Mechanics, Vol. 42, pp. 85-115.
• Chmielewski, C., and Jayaraman, K., 1992, "The effect of polymer extensibility on crossflow of polymer solutions through cylinder arrays", J. Rheology, Vol. 36(6), pp. 1105-1126.
• Dullien, F. A. L., 1992, "Porous Media: Fluid Transport and Pore Structure'" (2nd ed.), Academic Press.
• Durst, F., Haas, R., and Interthal, W., 1987, "The nature of flows through porous media", J. Non-Newtonian Fluid Mechanics, Vol. 22, pp. 169--189.
• Mizushina, T., and Usui, H., 1978, "Approximate solution of the boundary layer equations for the flow of a non-Newtonian fluid around a circular cylinder", Heat Transfer-Jpn. Res., Vol. 7(2), pp. 83--92.
• Mizushina, T., Usui, H., Veno, K., and Kato, K., 1978, "Experiments of pseudoplastic fluid cross flow around a circular cylinder", Heat Transfer-Jpn. Res., Vol. 7(3), pp. 92-101.
• Nieckele, A. O., Naccache, M. F., Souza Mendes, P. R., and Pires, L. F. G., 1995, "Flow of Bingham liquids through arrays of tubes", XIII Congresso Brasileiro de Engenharia Mecânica, in CD-ROM.
• Patankar S. V., 1980, "Numerical Heat Transfer and Fluid Flow", Hemisphere Publishing Corporation.
• Pires, L. F. G., and Nieckele, A. O., 1994a, "Hydrodynamics in tubes with curvilinear obstructions", Revista Brasileira de Ciência Mecânicas, vol XVI, No.4,pp. 407-414 (in portuguese).
• Pires L. F. G., and Nieckele, A. O., 1994b, "A numerical scheme for flow solution, based on contra-variant velocity components in non orthogonal coordinate", Proc. V Encontro Nacional de Ciências Térmicas, São Paulo, São Paulo, December 7-9, pp. 343-346 (in portuguese).
• Scheidegger, A. E.,1974, "The Physics of Flow Through Porous Media" (3rd ed.), University of Toronto Press.
• Settari, A., and Aziz, K., 1973, "A Generalization of the additive correction methods for the iterative solution of matrix equations", SIAM J. Numerical Analyses, Vol. 10, pp. 506-521.
• Souza Mendes, P. R., and Naccache, M. F., 1994a, "A modified Darcy's Law which accounts for extension-thickening rheological fluid behavior", V Encontro Nacional de Ciências Térmicas, São Paulo, São Paulo, December 7-9.
• Souza Mendes, P. R., and Naccache, M. F., 1994b, "A modified Darcy's Law which accounts for Binghamian rheological fluid behavior", XXII Encontro sobre Escoamento em Meios Porosos, Florianópolis, Santa Catarina, October 19-21.
• Van Doormaal, J. P., and Raithby, G. D., 1984, "Enhancements of the SIMPLE Method for prediction incompressible fluid flows", Numerical Heat Transfer, Vol. 7, pp. 147-163.
• Wu Y-S, Pruess K., and Whitherspoon, P. A., 1990, Flow and displacement of Bingham non-Newtonian fluids in porous media, 60th California Regional Meeting, Society of Petroleum Engineers, paper no. SPE 20051.
• Zukauskas, A., 1972, "Heat transfer from tubes in crossflow", Adv. Heat Transfer , Vol. 8, pp. 93-160. Academic, New York.

Publication Dates

• Publication in this collection
  18 Sept 2002
• Date of issue
  Mar 2002

History

• Received
  Sept 2000