The yield zone concept and its application on a 4:1 abrupt contraction for an apparent-yield-stress fluid

G.G. Ramos^†, E.J. Soares^† and R.L. Thompson^‡

^† Mechanical Engng. Department, UFES, Vitória, ES 29075-910, Brazil. ramos@ct.ufes.br, edson@ct.ufes.br
^‡ Mechanical Engng. Department, UFF, Niterói, RJ 24210-240, Brazil. rthompson@mec.uff.br

Abstract— The yield surface, a concept widely used in the literature that was born together with the Bingham model is discussed from the perspective of other models that are more related to apparent-yield-stress materials. As a consequence, it is necessary to define a yield zone, an intermediate transition region in the flow domain where plasticity manifests itself. A Galerkin Finite Element Method is used to investigate the performance of a SMD viscoplastic material through a 4:1 abrupt contraction. The influence of dimensionless parameters, like the Jump number and an equivalent to the Bingham number, on the size of the yield zone and on the pressure loss of the contraction are investigated.

Keywords— Viscoplastic Material; Abrupt Contraction Flow; Finite Element Method; Yield Surface; Yield Zone.

I INTRODUCTION

A. Motivation
The importance of non-Newtonian materials is becoming recognized in many important industries such as: the petroleum industry, the food industry, the cosmetic and pharmaceutical industry, etc. One first reasons is the increasing necessity of determining more accurately certain global quantities that are related to the efficiency of processes that have a traditional Newtonian calculation counterpart. It was not uncommon, in the past, to use Newtonian results as an approximation for non-Newtonian flows. A second important reason is that non-Newtonian materials are much more complex, and can exhibit features that are not present in Newtonian fluids world such like pseudoplasticity, viscoplasticity, elasticity, thixotropy. These features can be used to understand and produce new synthetized fluids and to conceive new methodologies in order to optimize standard processes (Sifontes et al., 2010). In the recovery of oil and drilling mud removal, among other procedures, non-Newtonian fluids are widely used (Morales et al., 2011). An approach to compute friction loss coefficients for some non-Newtonian fluids is given by Kfuri et al. (2011). The present work deals with one important class of non-Newtonian materials, namely the apparent-yield-stress fluid. These fluids are microstructure and can change significantly their viscosity from a higher level one, when the material is fully structured, to a lower one, when the material is almost unstructured. The concept of yield zone as an apparent-yield-stress fluid counterpart of the yield surface associated to yield-stress materials is presented and analyzed for the particular case of the flow in a 4:1 abrupt contraction.

B. Yield-stress materials
There is a large variety of materials that have a viscoplastic-like behavior (Pereira et al., 2010). Conceptually, a viscoplastic material is a material that possesses a yield-stress, τ₀, a stress limit below which the material does not flow. The first model that could predict this kind of behavior was proposed by Bingham (1922) as

(1)

where τ is the norm of the extra-stress tensor and is the norm of the deformation rate. The quantity η is the viscosity function associated to the Generalized Newtonian Liquid equation, . The relation presented in Eq. (1) cannot be written as a function for the stress. When , there are an infinity of values, in the interval [0,τ₀] that τ can assume and the viscosity is infinite. In order to represent the Bingham model numerically, it is common to use a bi-viscosity function using a very high viscosity value for low values of the shear rate.

C. Apparent-yield-stress fluids
There is a controversy in the literature, with respect to the yield stress. This controversy was triggered by Barnes and Walters (1985) where the existence of this entity was questioned. By conducting a series of experiments benefiting from the improvement of the accuracy of a new generation of rheometers, they claimed that some materials, considered as yield-stress materials, do flow below its apparent yield stress. In this sense, this critical stress, would mark a transition between two states of the structure level of the material. Below the (apparent) yield stress, the material is fully structured with a very high viscosity. Above this limit, the material starts a breakdown process that ends with a very low viscosity. An important review on the subject was made by Barnes (1999). A recent discussion on the existence or non-existence of the yield-stress is found in de Souza Mendes and Thompson (2013). As suggested by Hartnett and Hu (1989), Astarita (1990) the relevant issue is related to how we model the material for engineering purposes. Therefore, the yield stress can be considered "an engineering reality." Independently of the existence or not of yield-stress materials, one cannot ignore the existence of apparent-yield-stress fluids. In this case, there are not many models that can handle this kind of behavior, i.e. a very high level of viscosity below a critical stress (or shear rate) and a very low viscosity above this critical point.

A curious issue related to viscoplastic models is that when simulating these fluids, it is often convenient to employ a regularization parameter, m, since the discontinuity on the derivative of the viscosity function can be a source of numerical divergence. In fact, there is a finite zero-shear-rate viscosity associated to each value of regularization parameter chosen. A widely used regularization model when trying to capture the behavior of a Bingham material is the one proposed by Papanastasiou (1987) given by

(2)

transforming the stress-strain relation into a function. As m→∞ Papanastasiou's model approaches Bingham's one. In the manner it was conceived, i.e. to reproduce Bingham-like behavior, m has to have a high value. The literature recommends a dimensionless form m*=1000 (You et al., 2008). However, from the apparent-yield-stress fluids perspective, the zero-shear-rate viscosity, , of the Papanastasiou model is given by and so, its viscosity function can be rewritten as

(3)

giving to m a physical interpretation through its relation with η₀, the high viscosity level at low shear rates. Although there is a quantity in the model represented by the symbol τ₀, it is worth noticing that, in fact, there is not a true yield stress, since the material flows with a level of stress below this value. This quantity is, the so-called apparent yield stress. Therefore, this model is more aligned to the second way of thinking. Papanastasiuos model was used by several authors as: Dimakopoulos and Tsamopoulos (2003), Sousa et al. (2007), You et al. (2008), Thompson et al. (2010), Freitas et al. (2011), among others.

In order to construct a regularized version of the Hershell-Bulckley model, the corresponding adjustment, commonly used in the literature (Mitsoulis et al., 1993) is given by

(4)

However, there is a great difference between the adjusted Hershell-Bulckley model and the adjusted Bingham model (Papanastasiou), in the → 0 limit, when the fluid presents shear-thinning behavior (n<1). Although the inclusion of the regularization parameter in the case of the Hershell-Bulckley model do transform the stress-strain relation into a function, it does not present a bounded value for η₀.

D. The SMD equation
The previous result motivate de Souza Mendes and Dutra (2004) to propose another equation for the stress-strain relation, in consonance to the idea of the existence of a finite high viscosity value for low shear-rates, as

(5)

where a rheological regularization parameter m=η₀/τ₀ is used with the property that the identity is obtained. In this case, m is a characteristic time of the material with a very clear physical interpretation. The reciprocal of m, ₀, is, therefore, the value of the shear rate below which we have a Newtonian plateau.

For a better understanding of the main features of the SMD equation, we analyze in detail its main aspects of the model. For a deeper comprehension, the reader is referred to the works de Souza Mendes and Dutra (2004). The SMD equation can be written as a stress function of the shear-rate,

(6)

As shown in Fig. 1, the stress function can be roughly divided into three parts. The first one, for low values of the shear-rate, where stresses are below the "yield-stress," is a Newtonian region with a viscosity plateau of value h₀. The intermediate part is defined by a region where the stress has achieved a critical value, t₀, and remains close to this value for a certain range of shear-rate. In the third and last part the material behaves as a power-law fluid of parameters K, n.

Figure 1. Typical flow curve of an apparent-yield-stress fluid.

Another characteristic shear rate of the model is ₁, the value of the shear rate that marks the transition from the end of the microstructure breakdown to the beginning of the power-law behavior. It is defined as

(6)

E. The Yield Zone concept
The yield surface is a surface in space that limits yielded from unyielded regions. This concept is widely used in the literature and is common, when dealing with viscoplastic materials to show as a post-processing result, the regions where the material yields, bounded by this yield surface (Mitsoulis and Zisis, 2001; Freitas et al., 2011). The accurate position of the yield surface when a regularization-parameter model is used is a matter of deep discussion in the literature. Liu et al. (2002) and Frigaard and Nouar (2005) compared the regularization-parameter model proposed by Bercovier and Englement (1980) and Papanastasiou (1987). It is difficult, if not impossible to determine the position of the interface independently of the value of the regularization parameter. Putz et al. (2009) investigated the numerical implementation of the yield surface on a lubricating problem and found that there is an intrinsic paradox on determining this surface.

The difficulties on determining the accurate position of the yield surface motivate the introduction of a yield zone. The yield zone is a region in space that confines the major structure changes of an apparent-yield-stress fluid or a yield-stress material modeled with a regularization parameter. It begins where the fluid viscosity departs from its zero-shear-rate value and ends when the material is fully unstructured, for high values of shear rate and stress, when the viscosity is less influenced by the (apparent or not) yield stress. The yield zone gives a map of the regions where the breaking down process is taking place and separates yielded and unyielded regions in a smooth fashion, closely related to apparent-yield-stress fluids and computationally amenable.

When η₀ is very large, the value of the characteristic shear rate leads to t≈2τ₀. Hence, we used this result to define the boundaries of the yield zone.

One of the few analysis of the SMD fluid in a complex flow was conducted by Naccache and Barbosa (2007) using the commercial software FLUENT for the case of a planar expansion followed by a contraction. The results are given mainly as a function of the so-called Jump number. This parameter is a dimensionless number that measures the relative range the deformation rate assume when the stress intensity is close to t₀. It is defined as

(3)

II. PHYSICAL FORMULATION
A. Governing equations and boundary conditions
The problem chosen to evaluate the performance of the SMD equation and to illustrate the yield zone concept is the 4:1 abrupt contraction. The scheme of this problem is depicted in Fig. 2.

Figure 2. The axisymmetric 4:1 abrupt contraction considered.

The velocity and pressure fields are defined by the governing equations that impose conservation of mass, given by Eq. (4) below and momentum for a non-inertial incompressible fluid, Eqs. (5) and (6), together with the appropriate boundary conditions.

	(4)
	(5)
	(6)

where u and v are, respectively the axial and radial components of the velocity field and the quantities T_xx, T_rx, T_xr, T_rr, and T_θθ are components of the stress tensor.

The boundary conditions are: at the inlet and outlet, the flow is considered fully-developed and the pressure is imposed; along the symmetry axis both the shear stress and the radial velocity vanish; the no-slip and impermeability conditions are imposed along the walls.

III. NUMERICAL FORMULATION
The use of Finite Element methods is an important tool used to solve complex physical problems. In the present work the Finite Element method is used with a Galerkin formulation to solve the differential equations that govern the problem. Biquadratic basis functions Φ_j are used to represent the velocity and nodal coordinates, while linear discontinuous functions χ_j are employed to expand the pressure field. The velocity and pressure are represented in terms of appropriate basis functions

(7)

The coefficients of the expansions are the unknown of the problem

(8)

The corresponding weighted residuals of the Galerkin method related to conservation of mass and momentum in x and r directions are:

	(9)
	(10)
	(11)

where T is the stress tensor, n is the normal vector, ||J|| is the Jacobian of the transformation of the mesh, Ω is the domain of the numerical problem and Γ is the boundary of this domain.

As indicated above, the system of partial differential equations, and boundary conditions is reduced to a set of simultaneous algebraic equations for the coefficients of the basis functions of all the fields. This set is non-linear and sparse. It is solved by Newton's method. The linear system of equations, corresponding to each Newton iteration, was solved using a frontal solver.

The mesh used to solve the present problem is depicted in Fig. 3. It has 3801 nodes with 25404 degrees of freedom. A mesh refinement analysis was conducted and a less than 1% of difference was obtained.

Figure 3. Mesh used to simulate the problem considered.

IV. RESULTS
The fields of stress intensity for a fixed Jump number, J=2500, for different values of power-law index are shown in Fig. 4.

Figure 4. Stress intensity τ for fixed jump number J=2500 and yield stress τ₀=0.5; and power-law index varying from n=0.7 to n=1.0.

On the left, we can see an unyielded blue region and a yielded red region separated by a green region, the yield zone. Interestingly, the yield zone occupies a significant portion of the domain near the contraction. This fact justifies the use of the yield zone concept, since the attempt to define the yield surface would be subjected to a great uncertainty. This result can explain the problems and discussions in the literature of the subject. There is a small influence of the power-law index on the sizes of each region. However, comparing the result obtained with n=7 to the Newtonian result (n=1), we can see that the unyielded region (τ<τ₀) is smaller, while the red region (τ >2τ₀) is bigger. While the unyielded region is a part of the flow where shear is dominant, the red region just before the entrance of the small tube has a strong extensional component. On the right, we show the same result, but refining the stress level, so as to understand how the stresses evolve in the yield zone.

Figure 5 shows the field of deformation for different values of the Jump number. Here we can see that this dimensionless quantity has a strong influence on the field, as expected. The increase of J increases the portion of the field which is between ₀ and ₁. This change is accompanied by a decrease in the region of deformation rates below ₀ while the region above ₁becomes unaltered. This difference also shows the importance of the yield zone concept.

Figure 5. Deformation rate for different value of the Jump number, J=10, 100, 1000, and 2500.

The Couette correction, ΔP, a dimensionless number defined as the ratio of the pressure loss in the contraction and two times the stress at the wall of the small tube is given in Fig. 6, as a function of the critical shear stress. The critical shear stress is the stress evaluated at the critical shear rate, ₁, above which the material is mostly unstructured. The Couette correction gives a quantification of the loss of energy due to the abrupt contraction. For small values of the critical stress, there is significant variation of the Couette correction. However, after a certain value of critical stress, there is no change detectable. Another important point is that, differently from yield zone, this result is not sensitive to changes in the jump number.

Figure 6. Couette correction as a function of the critical shear stress for different values of the Jump number.

V. CONCLUSIONS
The present work analyzes the SMD equation, conceived to model viscoplastic behavior, in a 4:1 abrupt contraction from the numerical point of view using a Glerkin/Finite Element method. The objective of the investigation is to study the influence of the Jump number and the power-law exponent of the model on the size and shape of the yield zone. The yield zone, in contrast to the yield surface, is a region of the domain where stresses are close to the apparent yield stress. In this region, there can be a large range of deformation rate, depending on the value of the Jump number. Therefore, this region captures, in a more physical basis, the part of the domain where the major collapse of the structure of the material happens. This intermediate region bounds a region where the viscosity of the material is very high, corresponding to a high structural level of the material, and another region where the viscosity is very low, corresponding to a low structural level of the material.

In the range analyzed, namely 0.7≤n≤1, the power-law index had a timid influence on the general results. The variation of the Jump number has a strong influence on the size of the yield zone, but does not affect significantly the Couette-correction.

REFERENCES
1. Astarita, G., "The engineering reality of the yield stress," J. Rheol., 34, 275-277 (1990).         [ Links ]
2. Barnes, H., "The yield stress, a review," J. Non-Newtonian Fluid Mech., 81, 133-178 (1999).         [ Links ]
3. Barnes, H. and K. Walters, "The yield stress mith?," Rheol. Acta, 24, 323-326 (1985).         [ Links ]
4. Bercovier, M. and M. Engleman, "A finite-element method for incompressible non-Newtonian flows," J. Comp. Phys., 36, 313-326 (1980).         [ Links ]
5. Bingham, E., Fluidity and plasticity, McGraw-Hill, Boston (1922).         [ Links ]
6. de Souza Mendes, P.R. and E.S. Dutra, "Viscosity function for yield stress fluids," Appl. Rheol., 525, 1-26 (2004).         [ Links ]
7. de Souza Mendes, P.R. and R.L. Thompson, "A unified approach to model elasto-viscoplastic thixotropic yield stress materials and apparent-yield-stress fluids," Rheol. Acta, In Print (2013).         [ Links ]
8. Dimakopoulos, Y. and J. Tsamopoulos, "Transient displacement of a viscoplastic material by air in a straight and suddenly constricted tubes," J. Non-Newtonian Fluid Mech., 112, 43-75 (2003).         [ Links ]
9. Freitas, J.F., E.J. Soares and R.L. Thompson, "Immiscible Newtonian displacement by a viscoplastic material in a capillary plane channel," Rheol. Acta, 50, 403-422 (2011).         [ Links ]
10. Frigaard, I.A. and C. Nouar, "On the usage of viscosity regularization methods for viscoplastic fluid flow computation," J. Non-Newtonian Fluid Mech., 127, 1-26 (2005).         [ Links ]
11. Harnett, J.P. and R.Y.Z. Hu, "The yield stress - an engineering reality," J. Rheol., 34, 275-277 (1990).         [ Links ]
12. Kfuri, S.L.D., J.Q. Silva, E.J. Soares and R.L. Thompson, "Friction losses for power-law and voscoplastic materials in a entrance of a tube and an abrupt contraction," J. Petr. Sci. Eng., 76, 224-235 (2011).         [ Links ]
13. Liu, B., S. Muller and M.M. Denn, "Convergence of a regularization method for creeping flow of a Bingham material around a sphere," J. Non-Newtonian Fluid Mech., 102, 179-191 (2002).         [ Links ]
14. Mitsoulis, E., S.S. Abdali and N.C. Markatos, "Flow simulations of a Herschel-Bulkley fluids through extrusion dies," Can. J. Chem. Eng., 71, 147-160 (1993).         [ Links ]
15. Mitsoulis, E., and T. Zisis, "Flow of a Bingham plastic in a lid-driven cavity," J. Non-Newtonian Fluid Mech., 101, 173-180 (2001).         [ Links ]
16. Morales, C.J., U. Riebel, N.M. Guzmán and M. Guerra, "Formulation of water in paraffin emulsions," Latin Am. Appl. Res., 41, 105-112 (2011).         [ Links ]
17. Naccache, M.F. and R. Barbosa, "Creeping flow of viscoplastic materials through a planar expansion followed by a contraction," Mech. Research Comm., 34, 423-431 (2007).         [ Links ]
18. Papanastasiou, T., "Flows of materials with yield," J. Rheol., 81, 385-404 (1987).         [ Links ]
19. Pereira, F.A.R., C.H. Ataíde and M.A.S. Barrozo, "CFD approach using discrete phase model for annular flow analysis," Latin Am. Appl. Res., 40, 53-60 (2010).         [ Links ]
20. Putz, A., I.A. Frigaard and D.M. Martinez, " On the lubrication paradox and the use of regularization methods for lubrication flows," J. Non-Newtonian Fluid Mech., 163, 62-77 (2009).         [ Links ]
21. Sifontes, A.B., M. Urbina, F. Fajardo, L. Melo, L. García, M. Mediavilla, N. Carrión, J.L. Brito, P. Hernández, R. Solano, G. Mejias and A. Quintero, "Preparation of gamma-Alumina foams of high surface area employing the polyurethane sponge replica method," Latin Am. Appl. Res., 40, 185-191 (2010).         [ Links ]
22. Sousa, D.A., E.J. Soares, R. Queiroz, R.L. Thompson, "Numerical investigation on gas displacement of a shear-thinning liquid and a viscoplastic material in capillary tubes," J. Non-Newtonian Fluid Mech., 144, 149-159 (2007).         [ Links ]
23. Thompson, R.L., E.J. Soares and R.D.A. Bacchi, "Further remarks on numerical investigation on gas displacement of a shear-thinning liquid and a viscoplastic material in capillary tubes," J. Non-Newtonian Fluid Mech., 165, 448-452 (2010).         [ Links ]
24. You, Z., R.R. Huilgol and E. Mitsoulis, "Application of the Lambert W function to steady shearing flows of the Papanastasiou model," Int. J. Eng. Sci., 46, 799-880 (2008).         [ Links ]

Received: March 25, 2012
Accepted: March 2, 2013.
Recommended by Subject Editor: Adrian Lew.