Pore-scale modeling of viscoelastic flow in porous media using a Bautista–Manero fluid

This article examines the extensional flow and viscosity and the converging–diverging geometry as the basis of the peculiar viscoelastic behavior in porous media. The modified Bautista–Manero model, which successfully describes elasticity, thixotropic time dependency and shear-thinning, was used for modeling the flow of viscoelastic materials which also show thixotropic attributes. An algorithm, originally proposed by Philippe Tardy, that employs this model to simulate steady-state time-dependent flow was implemented in a non-Newtonian flow simulation code using pore-scale modeling. The simulation results using two topologically-complex networks confirmed the importance of the extensional flow and converging–diverging geometry on the behavior of non-Newtonian fluids in porous media. The analysis also identified a number of correct trends (qualitative and quantitative) and revealed the effect of various fluid and flow parameters on the flow process. The impact of some numerical parameters was also assessed and verified.


Introduction
The study of the flow of non-Newtonian fluids in porous media is of immense importance and serves a wide variety of practical applications in processes such as enhanced oil recovery from underground reservoirs, filtration of polymer solutions, and soil remediation through the removal of liquid pollutants. Viscoelasticity and thixotropy are two of the main features of non-Newtonian behavior. They are usually associated with the polymeric substances that are widely used in petroleum industry, chemical engineering systems, and many other scientific and industrial applications. Despite the fact that a massive amount of literature on these subjects do exist, there are few attempts to model these phenomena in connection with the flow through porous media. This is partly due to the mathematical complexity of these fluid models. Further difficulties are usually encountered in modeling the flow through porous media especially morphologically-complex ones.
In the recent years, a number of researchers have used porescale network modeling to describe the flow of complex fluids in porous media. In this context, Lopez et al. (2003), Lopez and Blunt (2004) investigated single-and two-phase flow of shear-thinning fluids in porous media using Carreau model in conjunction with network modeling. Sochi and Blunt (2008) followed a similar approach in their investigation of single-phase flow of Ellis and Herschel-Bulkley fluids through porous media. Thompson (2004, 2005) used a network model extracted from a computergenerated random sphere packing to investigate various aspects of non-Newtonian flow in packed beds. In this article we adopt pore-scale network modeling using two random networks to simulate the flow of Bautista-Manero fluids in porous media. Bautista-Manero is a reasonably-sophisticated model that is capable of describing various aspects of viscoelasticity and thixotropy among other non-Newtonian attributes.

Non-Newtonian fluids
Non-Newtonian fluids are commonly divided into three broad groups: time-independent, viscoelastic and time-dependent (Skelland, 1967;Chhabra and Richardson, 1999). The first group is characterized by the fact that the strain rate at a given point is solely dependent upon the instantaneous stress at that point. The second group includes the fluids that show partial elastic recovery upon the removal of a deforming stress. The third group consists of those fluids for which the strain rate is a function of both the magnitude and the duration of stress and possibly of the time lapse between consecutive applications of stress. A large number of models have been proposed in the literature to model all types of non-Newtonian fluids under various flow conditions. Most these models are basically empirical in nature and arise from curve-fitting exercises (Barnes et al., 1993).

Viscoelastic fluids
Viscoelastic substances exhibit a dual nature of behavior by showing signs of both viscous fluids and elastic solids. Polymeric fluids often show strong viscoelastic effects, which can include shear-thinning, extension thickening, viscoelastic normal stresses, and time-dependent rheology. The equations describing the flow of viscoelastic fluids consist of the basic laws of continuum mechanics and the rheological equation of state describing a particular fluid and relates the viscoelastic stress to the deformation history. Many differential and integral viscoelastic constitutive models have been proposed in the literature to describe the the observed viscoelastic phenomena. What is common to all these is the presence of at least one characteristic time parameter to account for the fluid memory (Larson, 1988(Larson, , 1999Keunings, 2004;Owens and Phillips, 2002).
The behavior of viscoelastic fluids is drastically different from that of Newtonian and inelastic non-Newtonian fluids. This includes the presence of normal stresses in shear flows, sensitivity to deformation type, and memory effects such as stress relaxation and time-dependent viscosity. These features underlie the observed peculiar viscoelastic phenomena such as rod-climbing (Weissenberg effect), die swell and open-channel siphon (Boger, 1987;Larson, 1988).
The behavior of viscoelastic fluids at any time is dependent on their recent deformation history, that is they possess a fading memory of their past. Indeed a material that has no memory cannot be elastic, since it has no way of remembering its original shape. Many materials are viscoelastic but at different time scales that may not be reached. Therefore the concept of a natural time of a material is important in characterizing the material as viscous or elastic. The ratio between the material time scale and the time scale of the flow is indicated by a non-dimensional number: the Deborah or the Weissenberg number (Barnes et al., 1993;Boger, 1987;Bird et al., 1987;Larson, 1988).
A common feature of viscoelastic fluids is stress relaxation after a sudden shearing displacement where stress overshoots to a maximum then starts decreasing exponentially and eventually settles to a steady-state value. This phenomenon also takes place on cessation of steady shear flow where stress decays over a finite measurable length of time. A defining characteristic associated with stress relaxation is the relaxation time which may be defined as the time required for the shear stress in a simple shear flow to return to zero under constant strain condition (Bird et al., 1987;Deiber, 1978;Larson, 1988).

Important aspects for viscoelastic flow in porous media
In the last few decades, a general consensus has emerged that in the flow of viscoelastic fluids through porous media elastic effects should arise, though their precise nature is unknown or controversial. Strong experimental evidence indicates that the flow of viscoelastic fluids through packed beds can exhibit rapid increases in the pressure drop, or an increase in the apparent viscosity, above that expected for a comparable purely viscous fluid. This increase has been attributed to the extensional nature of the flow field in the pores caused by the successive expansions and contractions that a fluid element experiences as it traverses the pore space. Even though the flow field at pore level is not an ideal extensional flow due to the presence of shear and rotation, the increase in flow resistance is referred to as an extension thickening effect (Phan-Thien and Khan, 1987;Plog, 2002;Deiber andSchowalter, 1981, Pilitsis andBeris, 1989).
There are two major interrelated aspects that have strong impact on the flow through porous media. These are extensional flow and converging-diverging geometry.

Extensional flow
One complexity in the flow in general and through porous media in particular usually arises from the coexistence of shear and extensional components; sometimes with the added complication of inertia. Pure shear or elongational flow is very much the matrix transpose x l network lower boundary in the non-Newtonian code x u network upper boundary in the non-Newtonian code exception in practical situations, especially in the flow through porous media. By far the most common situation is for mixed flow to occur where deformation rates have components parallel and perpendicular to the principal flow direction. In such flows, the elongational components may be associated with the converging-diverging flow paths (Barnes et al., 1993;Sorbie, 1991). A general consensus has emerged recently that the flow through packed beds has a substantial extensional component and typical polymer solutions exhibit strain hardening in extension, which is mainly responsible for the reported dramatic increases in pressure drop. Thus in principle the shear viscosity alone is inadequate to explain the observed excessive pressure gradients. It is therefore essential to know the relative importance of elastic and viscous effects or equivalently the relationship between normal and shear stresses for different shear rates (Carreau et al., 1997;Chhabra et al., 2001;Deiber, 1978).
Elongational flow is fundamentally different from shear, the material property characterizing the flow is not the viscosity, but the elongational viscosity. The behavior of the extensional viscosity function is very often qualitatively different from that of the shear viscosity function. For example, highly elastic polymer solutions that possess a viscosity which decreases monotonically in shear often exhibit an extensional viscosity that increases dramatically with strain rate. Thus, while the shear viscosity is shear-thinning the extensional viscosity is extension thickening (Larson, 1999;Barnes et al., 1993;Wapperom, 1996).

Converging-diverging geometry
An important aspect that characterizes the flow in porous media and makes it distinct from bulk is the presence of convergingdiverging flow paths. This geometric factor significantly affects the flow and accentuate elastic responses. Several arguments have been presented in the literature to explain the effect of convergingdiverging geometry on the flow behavior. Despite this diversity, there is a general consensus that in porous media the converging-diverging nature of the flow paths brings out both the extensional and shear properties of the fluid. The principal mode of deformation to which a material element is subjected as the flow converges into a constriction involves both a shearing of the material element and a stretching or elongation in the direction of flow, while in the diverging portion the flow involves both shearing and compression. The actual channel geometry determines the ratio of shearing to extensional contributions. In many realistic situations involving viscoelastic flows the extensional contribution is the most important of the two modes. As porous media flow involves elongational flow components, the coil-stretch phenomenon can also take place. Consequently, a suitable model pore geometry is one having converging and diverging sections which can reproduce the elongational nature of the flow, a feature that is not exhibited by straight capillary tubes (Marshall and Metzner, 1967;Deiber, 1978;Denys, 2003;Durst et al., 1987;Pilitsis and Beris, 1989).
For long time, the straight capillary tube has been the conventional model for porous media and packed beds. Despite the general success of this model with the Newtonian and inelastic non-Newtonian flow, its failure with elastic flow was remarkable. To redress the flaws of this model, the undulating tube and converging-diverging channel were proposed in order to include the elastic character of the flow. Various corrugated tube models with different simple geometries have been used as a first step to model the effect of converging-diverging geometry on the flow of viscoelastic fluids in porous media (e.g. Marshall and Metzner, 1967;Deiber and Schowalter, 1981;Souvaliotis and Beris, 1992). Those geometries include conically shaped sections, sinusoidal corrugation and abrupt expansions and contractions. Similarly, a bundle of converging-diverging tubes forms a better model for a porous medium in viscoelastic flow than the normally used bundle of straight capillary tubes, as the presence of diameter variations makes it possible to account for elongational contributions.
Many investigators have attempted to capture the role of the successive converging-diverging character of packed bed flow by numerically solving the flow equations in conduits of periodically varying cross-sections. Different opinions on the success of these models can be found in the literature. Examples include (Denys, 2003;Deiber, 1978;Talwar and Khomami, 1992;Chhabra et al., 2001;Pilitsis and Beris, 1989;Podolsak et al., 1997). With regards to modeling viscoelastic flow in regular or random networks of converging-diverging capillaries, very little work has been done.

Modeling the flow in porous media
The basic equations describing the flow of fluids consist of the basic laws of continuum mechanics which are the conservation principles of mass, energy and linear and angular momentum. These governing equations indicate how the mass, energy and momentum of the fluid change with position and time. The basic equations have to be supplemented by a suitable rheological equation of state, or constitutive equation describing a particular fluid, which is a differential or integral mathematical relationship that relates the extra stress tensor to the rate-of-strain tensor in general flow condition and closes the set of governing equations. One then solves the constitutive model together with the conservation laws using a suitable method to predict velocity and stress fields of the flows (Bird et al., 1987;Hulsen, 1990Hulsen, , 1986Carreau et al., 1997;Keunings, 2003Keunings, , 2004.
The mathematical description of the flow in porous media is extremely complex task and involves many approximations. So far, no analytical fluid mechanics solution to the flow through porous media has been found. Furthermore, such a solution is apparently out of reach for the foreseeable future. Therefore, to investigate the flow through porous media other methodologies have been developed, the main ones are the macroscopic continuum approach and the pore-scale numerical approach.
The widely used continuum models represent a simplified macroscopic approach in which the porous medium is treated as a continuum. All the complexities and fine details of the microscopic pore structure are absorbed into bulk terms like permeability that reflect average properties of the medium. Semi-empirical equations such as Darcy's law, Blake-Kozeny-Carman or Ergun equation fall into this category. Several continuum models are based in their derivation on the capillary bundle concept. The advantage of the continuum method is that it is simple and easy to implement with no computational cost. The disadvantage is that it does not account for the detailed physics at the pore level. One consequence of this is that in most cases it can only deal with steady-state situations with no time-dependent transient effects.
In the numerical approach, a detailed description of the porous medium at pore-scale level is adopted and the relevant physics of flow at this level is applied. To find the solution, numerical methods, such as finite volume and finite difference, usually in conjunction with computational implementation are used. The advantage of the numerical method is that it is the most direct approach to describe the physical situation and the closest to full analytical solution. It is also capable, in principle at least, to deal with timedependent transient situations. The disadvantage is that it requires a detailed pore space description. Moreover, it is usually very complex and hard to implement and has a huge computational cost with serious convergence difficulties. Due to these complexities, the flow processes and porous media that are currently within the reach of numerical investigation are the most simple ones.
Pore-scale network modeling is a relatively novel method developed to deal with flow through porous media. It can be seen as a compromise between these two extreme approaches as it partly accounts for the physics and void space description at the pore level with reasonable and generally affordable computational cost. Network modeling can be used to describe a wide range of properties from capillary pressure characteristics to interfacial area and mass transfer coefficients. The void space is described as a network of pores connected by throats. The pores and throats are assigned some idealized geometry, and rules which determine the transport properties in these elements are incorporated in the network to compute effective transport properties on a mesoscopic scale. The appropriate pore-scale physics combined with a geologically representative description of the pore space gives models that can successfully predict average behavior (Blunt et al., 2002;Blunt, 2001).
In our investigation to the flow of Bautista-Manero fluids in porous media we use network modeling. Our model uses threedimensional networks built from a topologically-equivalent three-dimensional voxel image of the pore space with the pore sizes, shapes and connectivity reflecting the real medium. Pores and throats are modeled as having triangular, square or circular cross-section by assigning a shape factor which is the ratio of the area to the perimeter squared and obtained from the pore space image. Most of the network elements are not circular. To account for the non-circularity when calculating the volumetric flow rate analytically or numerically for a cylindrical capillary, an equivalent radius R eq is defined: where the geometric conductance, G, is obtained empirically from numerical simulation. Two networks obtained from Statoil and representing two different porous media have been used: a sand pack and a Berea sandstone. These networks are constructed by Øren and coworkers (Øren et al., 1997;Øren and Bakke, 2003) from voxel images generated by simulating the geological processes by which the porous medium was formed. The physical and statistical properties of the networks are given in Sochi and Blunt (2008). Assuming a laminar, isothermal and incompressible flow at low Reynolds number, the only equations that need to be considered are the constitutive equation for the particular fluid and the conservation of volume as an expression for the conservation of mass. Because initially the pressure drop in each network element is not known, an iterative method is used. This starts by assigning an effective viscosity to each network element. The effective viscosity is defined as that viscosity which makes Poiseulle's equation fit any set of laminar flow conditions for time-independent fluids (Skelland, 1967). By invoking the conservation of volume for incompressible fluid, the pressure field across the entire network is solved using a numerical solver (Ruge and Stüben, 1987). Knowing the pressure drops, the effective viscosity of each element is updated using the expression for the flow rate with a pseudo-Poiseuille definition. The pressure field is then recomputed using the updated viscosities and the iteration continues until convergence is achieved when a specified error tolerance in total flow rate between two consecutive iteration cycles is reached. Finally, the total volumetric flow rate and the apparent viscosity in porous media, defined as the viscosity calculated from the Darcy's law, are obtained.
With regards to modeling the flow in porous media of complex fluids which have time dependency due to thixotropic or elastic nature, there are three major difficulties: The difficulty of tracking the fluid elements in the pores and throats and identifying their deformation history, as the path followed by these elements is random and can have several unpredictable outcomes.
The mixing of fluid elements with various deformation history in the individual pores and throats. As a result, the viscosity is not a well-defined property of the fluid in the pores and throats. The change of viscosity along the streamline since the deformation history is constantly changing over the path of each fluid element.
In the current work, we deal only with one case of steady-state viscoelastic flow, and hence we did not consider these complications in depth. Consequently, the tracking of fluid elements or flow history in the network and other dynamic aspects are not implemented in the non-Newtonian code. However, time-dependent effects in steady-state conditions are accounted for in the Tardy algorithm which is implemented in the code and will be presented in detail in Section 3.1.
In our modeling approach, to solve the pressure field across a network of n nodes we write n equations in n unknowns which are the pressure values at the nodes. The essence of these equations is the continuity of flow of incompressible fluid at each node in the absence of source and sink. We solve this set of equations subject to the boundary conditions which are the pressures at the inlet and outlet. This unique solution is 'consistent' and 'stable' as it is the only mathematically acceptable solution to the problem, and, assuming the modeling process and the mathematical technicalities are correct, should mimic the unique physical reality of the pressure field in the porous medium.

Bautista-Manero model
This is a relatively simple model that combines the Oldroyd-B constitutive equation for viscoelasticity and the Fredrickson's kinetic equation for flow-induced structural changes usually associated with thixotropy. The model requires six parameters that have physical significance and can be estimated from rheological measurements. These parameters are the low and high shear rate viscosities, the elastic modulus, the relaxation time, and two other constants describing the build up and break down of viscosity.
The Oldroyd-B model is a simplification of the more elaborate and rarely used Oldroyd 8-constant model which also contains the upper convected, the lower convected, and the corotational Maxwell equations as special cases. Oldroyd-B is the second simplest nonlinear viscoelastic model and is apparently the most popular in viscoelastic flow modeling and simulation. It is the nonlinear equivalent of the linear Jeffreys model, and hence it takes account of frame invariance in the nonlinear regime. Consequently, in the linear viscoelastic regime the Oldroyd-B model reduces to the linear Jeffreys model. The Oldroyd-B model can be obtained by replacing the partial time derivatives in the differential form of the Jeffreys model with the upper convected time derivatives (Bird et al., 1987) s þ k 1 s where s is the extra stress tensor, k 1 is the relaxation time, k 2 is the retardation time, l 0 is the low-shear viscosity, _ c is the rate-of-strain tensor, and s 5 is the upper convected time derivative of the stress tensor: where t is time, v ¼ ðv x ; v y ; v z Þ is the fluid velocity vector, ðÁÞ T is the transpose of the tensor and rv is the fluid velocity gradient tensor defined by rv ¼ @vx @x @vx @y @vx @z @vy @x @vy @y @vy @z @vz @x @vz @y @vz @z Similarly, _ c 5 is the upper convected time derivative of the rateof-strain tensor given by: The kinetic equation of Fredrickson that accounts for the destruction and construction of structure is given by where l is the non-Newtonian viscosity, t is the time of deformation, k is the relaxation time upon the cessation of steady flow, l 0 and l 1 are the viscosities at zero and infinite shear rates, respectively, k is a parameter that is related to a critical stress value below which the material exhibits primary creep, s is the stress tensor and _ c is the rate-of-strain tensor. In this model, k is a structural relaxation time, whereas k is a kinetic constant for structure break down.
The elastic modulus G o is related to these parameters by G o ¼ l=k (Bautista et al., 1999(Bautista et al., , 2000Manero et al., 2002;Tardy and Anderson, 2005 of viscoelastic systems that also exhibit thixotropy and rheopexy under shear flow. The model predicts creep behavior, stress relaxation and the presence of thixotropic loops when the sample is subjected to transient stress cycles. The Bautista-Manero model has also been found to fit steady shear, oscillatory and transient measurements of viscoelastic solutions (Bautista et al., 2000(Bautista et al., , 1999Tardy and Anderson, 2005;Manero et al., 2002).

Tardy algorithm
This algorithm is proposed by Philippe Tardy to compute the pressure drop-flow rate relationship for the steady-state flow of a Bautista-Manero fluid in simple capillary network models. The bulk rheology of the fluid and the dimensions of the capillaries making up the network are used as inputs to the models (Tardy and Anderson, 2005). In the following paragraphs we outline the basic components of this algorithm and the logic behind it. This Table 1 Some values of the worm-like micellar system studied by Anderson et al. (2006) which is a solution of surfactant concentrate (a mixture of EHAC and 2-propanol) in an aqueous solution of potassium chloride.

Parameter Value
Go will be followed by some mathematical and technical details related to the implementation of this algorithm in our non-Newtonian code.
The flow in a single capillary can be described by the following general relation where Q is the volumetric flow rate, G 0 is the flow conductance and DP is the pressure drop. For a particular capillary and a specific fluid, G 0 is given by For a network of capillaries, a set of equations representing the capillaries and satisfying mass conservation should be solved simultaneously to produce a consistent pressure field, as presented in Section 2. For Newtonian fluid, a single iteration is needed to solve the pressure field since the conductance is known in advance as the viscosity is constant. For purely viscous non-Newtonian fluid, we start with an initial guess for the viscosity, as it is unknown and pressure-dependent, and solve the pressure field iteratively updating the viscosity after each iteration cycle until convergence is reached. For memory fluids, the dependence on time must be taken into account when solving the pressure field iteratively. Apparently, there is no general strategy to deal with such situation. However, for the steady-state flow of memory fluids a sensible approach is to start with an initial guess for the flow rate and iterate, considering the effect of the local pressure and viscosity variation due to converging-diverging geometry, until convergence is achieved. This approach is adopted by Tardy to find the flow of a Bautista-Manero fluid in a simple capillary network model. In general terms, the Tardy algorithm can be summarized as follows (Tardy and Anderson, 2005) For fluids without memory the capillary is unambiguously defined by its radius and length. For fluids with memory, where going from one section to another with different radius is important, the capillary should be modeled with contraction to account for the effect of converging-diverging geometry on the flow. The reason is that the effects of fluid memory take place on going through a radius change, as this change induces a change in shear rate and generation of extensional flow fields with a consequent viscoelastic and thixotropic effects. Examples of the converging-diverging geometries are given in Fig. 18. Each capillary is discretized in the flow direction and a discretized form of the flow equations is used assuming a prior knowledge of the stress and viscosity at the inlet of the network. Starting with an initial guess for the flow rate and using iterative technique, the pressure drop as a function of the flow rate is found for each capillary. Finally, the pressure field for the whole network is found iteratively until convergence is achieved. Once the pressure field is found the flow rate through each capillary in the network can be computed and the total flow rate through the network can be determined by summing and averaging the flow through the inlet and outlet capillaries.
A modified version of the Tardy algorithm was implemented in our non-Newtonian code. In the following, we outline the main steps of this algorithm and its implementation.
From a one-dimensional steady-state Fredrickson equation in which the partial time derivative is written in the form @ @t ¼ V @ @x , the following equation can be obtained  Fig. 9. The Tardy algorithm Berea results for Go ¼ 1:0 Pa; l 1 ¼ l 0 ¼ 0:1 Pa s; k ¼ 1:0 s; k ¼ 10 À5 Pa À1 ; f e ¼ 1:0; m ¼ 10 slices, with varying fm (0.6, 0.8, 1.0, 1.2 and 1.4). In this formulation, the fluid speed V is given by Q =pr 2 and the average shear rate in the tube with radius r is given by Q =pr 3 .
By similar argument, another simplified equation can be obtained from the Oldroyd-B model The converging-diverging feature of the capillary is implemented in the form of a parabolic profile as outlined in Appendix A. Each capillary in the network is discretized into m slices, each with width dx ¼ L=m where L is the capillary length.
From Eq. (10), applying the simplified assumption ds dx ¼ ðs 2 Às 1 Þ dx where the subscripts 1 and 2 stand for the inlet and outlet of the slice, respectively, we obtain an expression for s 2 in terms of s 1 and l 2 From Eqs. (9) and (11) The algorithm starts by assuming a Newtonian flow in a network of straight capillaries. Accordingly, the volumetric flow rate Q, and consequently V and _ c, for each capillary are obtained.
Starting from the inlet of the network where the viscosity and the stress are assumed to have known values of l 0 and RDP=2L for each capillary, respectively, the computing of the non-Newtonian flow in a converging-diverging geometry takes place in each capillary independently by calculating l 2 and s 2 slice by slice, where the values from the previous slice are used for l 1 and s 1 of the current slice.
For the capillaries which are not at the inlet of the network, the initial values of the viscosity and stress at the inlet of the capillary are found by computing the Q-weighted average from all the capillaries that feed into the pore which is at the inlet of the corresponding capillary.
To find l 2 of a slice, a bisection numerical method is used. To eliminate possible non-physical roots, the interval for the accepted root is set between zero and 3l 0 with l 0 used in the case of failure. These conditions are logical as long as the slice is reasonably thin and the flow and fluid are physically viable.
In the case of a convergence failure, error messages are issued to inform the user. No failure has been detected during the many runs of this algorithm. Moreover, extensive sample inspection of the l 2 values has been carried out and proved to be sensible and realistic.
The value found for the l 2 is used in conjunction with Eq. (11) to find s 2 which is needed as an input to the next slice. Averaging the value of l 1 and l 2 , the viscosity for the slice is found and used with Poiseuille law to find the pressure drop across the slice.
The total pressure drop across the whole capillary is computed by summing up the pressure drops across its individual slices. This total is used with Poiseuille law to find the effective viscosity for the capillary as a whole.
Knowing the effective viscosities for all capillaries of the network, the pressure field is solved iteratively using an algebraic multi-grid solver, and hence the total volumetric flow rate from the network and the apparent viscosity are found.
A flow chart of the modified Tardy algorithm is presented in Fig. 1.

Initial results of the modified tardy algorithm
The modified Tardy algorithm was tested and assessed. Various qualitative aspects were verified and proved to be correct. Some general results and conclusions are outlined below with sample graphs and data for the sand pack and Berea networks using a calculation box with x l ¼ 0:5 and x u ¼ 0:95. We would like to remark that despite our effort to use typical values for the parameters and variables, in some cases we were forced, for demonstration purposes, to use eccentric values to accentuate the features of interest. In Table 1 we present some values for the Bautista-Manero model parameters as obtained from the worm-like micellar system studied by Anderson et al. (2006) to give an idea of the parameter ranges in real fluids. The system is a solution of a surfactant concentrate (a mixture of the cationic surfactant erucyl bis(hydroxyethyl)methylammonium chloride (EHAC) and 2-propanol in a 3:2 weight ratio) in an aqueous solution of potassium chloride.

Convergence-divergence
A dilatant effect due to the converging-diverging feature has been detected relative to a network of straight capillaries. The viscosity increase is a natural viscoelastic response to the radius tightening at the middle. As the corrugation feature of the tubes is exacerbated by narrowing the radius at the middle, the viscoelastic dilatant behavior is intensified. The natural explanation is that the increase in apparent viscosity should be proportionate to the magnitude of tightening. This feature is presented in Figs and 3 for the sand pack and Berea networks, respectively, on a log-log scale.

Diverging-converging
The investigation of the effect of diverging-converging geometry by expanding the radius of the capillaries at the middle revealed a thinning effect relative to the straight and convergingdiverging geometries, as seen in Figs. 4 and 5 for the sand pack and Berea networks, respectively, on a log-log scale. This is due to viscoelastic response in the opposite sense to that of converging-diverging by enlarging the radius at the middle.

Number of slices
As the number of slices of the capillaries increases, the algorithm converges to a stable and constant solution within acceptable numerical errors. This indicates that the numerical aspects of the algorithm are functioning correctly because the effect of discretization errors is expected to diminish by increasing the number of slices and hence decreasing their width. A sample graph of apparent viscosity versus number of slices for a typical data point is presented in Fig. 6 for the sand pack and in Fig. 7 for Berea.

Boger fluid
A Boger fluid behavior was observed when setting l 0 ¼ l 1 .
However, the apparent viscosity increased as the convergingdiverging feature is intensified. This feature is demonstrated in Fig. 8 for the sand pack and in Fig. 9 for Berea on a log-log scale. Since Boger fluid is a limiting and obvious case, this behavior indicates that the model, as implemented, is well-behaved. The viscosity increase is a natural viscoelastic response to the radius tightening at the middle, as discussed earlier.

Elastic modulus
The effect of the elastic modulus G o was investigated for shearthinning fluids, i.e. l 0 > l 1 , by varying this parameter over several orders of magnitude while holding the others constant. It was observed that by increasing G o , the high-shear apparent viscosities were increased while the low-shear viscosities remained constant and have not been affected. However, the increase at high-shear rates has almost reached a saturation point where beyond some limit the apparent viscosities converged to certain values despite a large increase in G o . A sample of the results for this investigation is presented in Fig. 10 for the sand pack and in Fig. 11 for Berea on a log-log scale. The stability at low-shear rates is to be expected because at low flow rate regimes near the lower Newtonian plateau the non-Newtonian effects due to elastic modulus are negligible. As the elastic modulus increases, an increase in apparent viscosity due to elastic effects contributed by elastic modulus occurs. Eventually, a saturation will be reached when the contribution of this factor is controlled by other dominant factors and mechanisms from the porous medium and flow regime.
3.2.6. Shear-thickening A slight shear-thickening effect was observed when setting l 0 < l 1 while holding the other parameters constant. The effect of increasing G o on the apparent viscosities at high-shear rates was similar to the effect observed for the shear-thinning fluids though it was at a smaller scale. However, the low-shear viscosities 1.0E-01 1.0E+00

1.0E+01
Darcy Velocity / (m/s) Apparent Viscosity / (Pa.s) Decreasing k Fig. 16. The Tardy algorithm sand pack results for Go ¼ 1:0 Pa; l 1 ¼ 0:001 Pa s; l 0 ¼ 1:0 Pa s; k ¼ 10 s; fe ¼ 1:0; f m ¼ 0:5; m ¼ 10 slices, with varying k (10 À3 , 10 À4 , 10 À5 , 10 À6 and 10 À7 Pa À1 ). are not affected, as in the case of shear-thinning fluids. A convergence for the apparent viscosities at high-shear rates for large values of G o was also observed as for shear-thinning. A sample of the results obtained in this investigation is presented in Figs. 12 and 13 for the sand pack and Berea networks, respectively, on a log-log scale. It should be remarked that as the Bautista-Manero is originally a shear-thinning model, it should not be expected to fullypredict shear-thickening phenomenon. However, the observed behavior is a good sign for our model because as we push the algorithm beyond its limit, the results are still qualitatively reasonable.

Relaxation time
The effect of the structural relaxation time k was investigated for shear-thinning fluids by varying this parameter over several orders of magnitude while holding the others constant. It was observed that by increasing the structural relaxation time, the apparent viscosities were steadily decreased. However, the decrease at high-shear rates has almost reached a saturation point where beyond some limit the apparent viscosities converged to certain values despite a large increase in k. A sample of the results is given in Fig. 14 for the sand pack and Fig. 15 for Berea on a loglog scale. The effects of relaxation time are expected to ease as the relaxation time increases beyond a limit such that the effect of interaction with capillary constriction is negligible. Despite the fact that this feature requires extensive investigation for quantitative confirmation, the observed behavior seems qualitatively reasonable considering the flow regimes and the size of relaxation times of the sample results.

Kinetic parameter
The effect of the kinetic parameter for structure break down k was investigated for shear-thinning fluids by varying this parameter over several orders of magnitude while holding the others constant. It was observed that by increasing the kinetic parameter, the apparent viscosities generally decreased. However, in some cases the low-shear viscosities has not been substantially affected. A sample of the results is given in Figs. 16 and 17 for the sand pack and Berea networks, respectively, on a log-log scale. The decrease in apparent viscosity on increasing the kinetic parameter is a natural response as the parameter quantifies structure break down and hence reflects thinning mechanisms. It is a general trend that the low flow rate regimes near the lower Newtonian plateau are usually less influenced by non-Newtonian effects. It will therefore be natural that the low-shear viscosities in some cases experienced minor changes.

Conclusions
A non-Newtonian flow simulation code based on pore-scale network modeling has been extended to account for viscoelastic behavior using a Bautista-Manero fluid model. The basis of the implementation is a numerical algorithm originally suggested by Philippe Tardy. A modified version of this algorithm was used to investigate the effect of converging-diverging geometry and extensional flow on the steady-state flow. The implementation of this model has been examined using a sand pack and a Berea networks. Various physical and computational aspects that have qualitative and quantitative implications were investigated and assessed. These include converging-diverging geometry, numerical convergence, limiting cases (e.g. Boger fluid and shear-thickening), elastic modulus, relaxation time and kinetic parameter. The results confirmed the importance of the extensional flow and convergingdiverging geometry on the behavior of non-Newtonian fluids in porous media, and identified a number of correct trends. A conclusion was reached that the current implementation is a proper basis for investigating the steady-state viscoelastic effects and some thixotropic effects due to converging-diverging geometry in porous media. This implementation can be used as a suitable foundation for further development. The future work can include elaborating the modified Tardy algorithm and thoroughly investigating its viscoelastic and thixotropic predictions in more rigorous quantitative terms. The code (called Non-Newtonian Code 2) with complete documentation and the networks can be obtained from this URL: http://www3.imperial.ac.uk/earthscienceandengineering/research/perm/porescalemodelling/software/non-newtonian% 20code.

Appendix A. Converging-diverging geometry and tube discretization
There are many simplified converging-diverging geometries for corrugated capillaries which can be used to model the flow of viscoelastic fluids in porous media; some suggestions are presented in Fig. 18. In this study we adopted the paraboloid and hence developed simple formulae to find the radius as a function of the distance along the tube axis, assuming cylindrical coordinate system, as shown in Fig. 19.
For the paraboloid depicted in Fig. 19, we have Since rðÀL=2Þ ¼ rðL=2Þ ¼ R max and rð0Þ ¼ R min ð14Þ the paraboloid is uniquely identified by these three points. On substituting and solving these equations simultaneously, we obtain that is In the modified Tardy algorithm for steady-state viscoelastic flow as implemented in our non-Newtonian code, when a capillary is discretized into m slices the radius rðzÞ is sampled at m z-points given by z ¼ À L 2 þ k L m ðk ¼ 1; . . . ; mÞ ð 17Þ More complex polynomial and sinusoidal converging-diverging profiles can also be modeled using this approach.