J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
A general transformation procedure for differential viscoelastic models G. Mompean a,∗ , R.L. Thompson b , P.R. Souza Mendes b a
b
Laboratoire de Mécanique de Lille, Polytech’Lille, Université des Sciences et Technologies de Lille, URM-CNRS 8107, Cité Scientifique, 59655 Villeneuve d’Ascq Cedex, France Department of Mechanical Engineering, Pontifícia Universidade Católica, Rio de Janeiro RJ22453-900, Brazil Received 22 April 2002; received in revised form 6 February 2003
Abstract A general transformation procedure (GTP) for viscoelastic differential constitutive equations is proposed. This procedure is devised to transform a differential model (e.g. Oldroyd-B, White–Metzner (WM) and Phan-Thien– Tanner (PTT)) into another model which has just one differential equation to be solved. The resulting model is suitable for three-dimensional and time-dependent flows. A new generalized objective time derivative is employed in the procedure, lending important features to the transformed model. For example, for the case of the Oldroyd-B model, the transformed model is capable of predicting non-zero second normal stress differences. Calculations of the stress field have been performed with an Oldroyd-B based GTP model (the GTP[OldB] model) for a 4:1 contraction plane flow. The results obtained have been analyzed with the aid of a flow-type classifier, and were found to be in good agreement with the ones obtained with the original model. © 2003 Elsevier Science B.V. All rights reserved. Keywords: Transformation; Differential; Viscoelastic
1. Introduction Numerical prediction of viscoelastic fluids flowing through three-dimensional complex geometries is an area of great interest, encompassing industrial and scientific applications. Viscoelastic liquids are commonly found in the industries of plastics, food, paints, petroleum, to name a few. In these industrial applications, viscoelastic liquids often flow through complex three-dimensional geometric configurations, and their viscoelastic nature typically requires multi-mode constitutive models to truthfully represent their mechanical behavior. Examples where time-dependent effects are present also abound in industrial processes. ∗
Corresponding author. E-mail addresses:
[email protected] (G. Mompean),
[email protected] (R.L. Thompson),
[email protected] (P.R. Souza Mendes). 0377-0257/03/$ – see front matter © 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0377-0257(03)00042-9
152
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
A great challenge while performing numerical simulations of these flows is to solve the mass and momentum conservation equations in conjunction with a constitutive equation for stress. The resulting system of equations is a mathematical description of a physical problem which involves complex phenomena and effects, such as elasticity, inertia and diffusion. For example, in the case of a three-dimensional flow of a viscoelastic fluid described by any differential single-mode constitutive model, the task is to solve 10 coupled partial differential equations (namely, six for the non-Newtonian extra-stress tensor, three for momentum conservation, and one for mass conservation). This number can increase rapidly for multi-mode fluid models (e.g. [3]). The computing time and the storage memory are important parameters that can reach huge values for three-dimensional numerical simulation of multi-mode viscoelastic fluids. Even with powerful modern computers and the most efficient algorithms, performing such simulations still constitutes a formidable task. Due to this reason, the quest for simpler viscoelastic models for numerical simulations has become an important area of research in the last years. The use of the general transformation procedure (GTP) presented in this work goes in this direction. The aim here is to propose a method for obtaining extra-stress models starting from differential constitutive equations, such as, for example, the Upper-Convected Maxwell (UCM), Oldroyd-B, White–Metzner (WM) and Phan-Thien–Tanner (PTT) models. The procedure is based on the transformation of the constitutive equation, keeping its elasticity prediction capability. The transformation reduces the set of six differential equations found in these models into one only. This class of models has the obvious advantage, for engineering applications, of requiring much less computational time and storage memory, which is especially important for three-dimensional simulations. This kind of approach has been extensively employed by the turbulence modeling community to simplify the differential Reynolds stress equation [26], while keeping the main features of the physics of turbulence. Inspired in the analogy between viscoelastic fluids and turbulence modeling mentioned by Rivlin [28], this approach was recently applied to flows of viscoelastic liquids by Mompean et al. [21] and Mompean [19], to obtain constitutive models via polynomial expansions [26]. The approach is interesting because of its capability to predict viscoelastic effects. However, the models proposed in these recent articles are not frame indifferent, since the vorticity tensor W has been used as a basis tensor in the polynomial expansion. Objectivity is attained in the present work by the modification of a previously used hypothesis for the advection of non-Newtonian extra-stress. In this connection, two new kinematic tensors (see [1]) are needed, namely (i) Ω, the rate-of-rotation of the principal directions of S, the rate-of-deformation tensor, ¯ = W −Ω. The latter is the relative-rate-of-rotation tensor, which is objective and measures the and (ii) W rate-of-rotation of a material particle as seen by an observer fixed to the principal axes of S. Inspired in the ¯ early work by Schunk and Scriven [34], Souza Mendes et al. [35] and Thompson et al. [37] employed W to propose algebraic constitutive models with enhanced capability for predicting rheological steady-state material functions, but not capable of predicting transient elasticity effects. Flows of viscoelastic liquids inside or around complex geometries having singularities are rather difficult to simulate numerically. For this reason, the flow through abrupt planar or axisymmetric contractions is often employed as a test flow for numerical simulations in the field of non-Newtonian fluid mechanics, and will also be used here to evaluate the present formulation. Even for the Newtonian case, the flow originated by this geometry is complex as the fluid is submitted to a mix of shear, extension and solid body motion. The singularity present in the sharp corner of this geometry is responsible for many difficulties that arise in numerical simulations of viscoelastic flows. The strong gradients of pressure and extra-stress components near the singularity can introduce numerical errors.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
153
Several authors have studied numerical simulations of viscoelastic flows through contractions. Keunings [16], using a finite element method, made simulations in an axisymmetric 4:1 contraction for Newtonian, Maxwell and Leonov-like fluids and showed the existence of a critical value for the Deborah number. The critical value reached depends on the grid employed. In Purnode and Crochet [27], the viscoelastic flow of a polyacrylamide solution was simulated using a single-mode FENE-P model, which predicts shear thinning and a finite extensional viscosity. The results were compared with the experimental data of Evans and Walters [7,8]. These comparisons showed that numerical results showing streamlines for different flow conditions are in qualitative agreement with experimental results. The influence of several parameters (ratio of contraction, sharpness of the corner, polymer concentration and inertia), was analyzed for the mechanism of vortex enhancement and the occurrence of lip vortex in planar contractions. Because only qualitative agreement was achieved between numerical simulations and experimental results, one of the conclusions of these studies was the need of more appropriate constitutive equations. Sato and Richardson [33] used a mixed numerical approach combining the finite-element and finitevolume methods in order to save computing time and storage memory. This method, which consist of solving the momentum and mass conservation by the finite-volume method and the extra-stress tensor by a finite-element method, was applied to a 4:1 axisymmetric contraction. Sasmal [32] used a finite-volume method with a pseudo-transient elastic stress term. In his calculations, a creeping flow is assumed, neglecting three terms: (i) the non-linear term, (ii) the time-dependent term in the momentum equation, and (iii) the time-dependent term in the constitutive equation for the UCM. Matallah et al. [18] showed stress-splitting schemes for viscoelastic flows through planar contractions with rounded and sharp re-entrance corners. They presented a large review of numerical solutions and vortex behavior applied to planar contractions. Recently, Phillips and Williams [25] presented a semi-Lagrangian finite-volume method applied to a planar contraction using the Oldroyd-B fluid. The results concerning vortex enhancement with increasing Deborah number are presented and analyzed. In the papers mentioned above, the authors employed differential constitutive equations in order to save computational effort, whenever an axis or a plane of symmetry could be identified in a given geometrical configuration, flow symmetry boundary conditions were applied to solve the flow equations. While this type of boundary conditions is physically correct for laminar flows of Newtonian fluids, it is not always true for flows of viscoelastic fluids, as shown by Boger [4]. For example, in the case of planar contractions, numerical simulations are often done by assuming one symmetry plane for two-dimensional flow, and two symmetry planes for the three-dimensional situation [20,29–31,38]. These assumptions can conceal physical solutions that would be found if considering the entire flow domain, which is often prohibitive due to the usual complexity of the constitutive equations. In this study, we present solutions for the flow through a planar contraction, using a new model to be described shortly. The paper is organized as follows. In Section 2, the conservation equations are presented. The general transformation procedure is introduced in Section 3 and applied to the Oldroyd-B and PTT models. The finite-volume numerical method used to discretize the equations is discussed in Section 4. Numerical results for the flow through the 4:1 contraction are presented in Section 5, together with a discussion of some important features of the proposed model. Finally, conclusions are given in Section 6.
154
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
Fig. 1. Geometrical parameters of the 4:1 contraction.
2. Governing equations The mass and momentum conservation equations coupled with the constitutive relation for the nonNewtonian extra-stress components of an Oldroyd-B fluid are considered in this section. In order to test the performance of a GTP derived model, we chose to obtain numerical solutions for flow through an abrupt contraction, as illustrated in Fig. 1. The velocity field as well as the independent variables are made non-dimensional using a characteristic velocity scale U and length scale L, taken respectively as the average velocity in the downstream half channel and the width of the downstream half channel (H ). The pressure, and non-Newtonian extra-stress are scaled with η∗ U/L. The viscosity η∗ is equal to the sum of the Newtonian viscosity η0∗ and the polymeric viscosity η1∗ . Therefore, the Reynolds number is defined as Re = ρUL/η∗ , where ρ is the density. The viscosities are made dimensionless as ηi = ηi∗ /η∗ , i = 0, 1. The Deborah number is given by De = λγ˙ , where λ is the relaxation time of the viscoelastic fluid, and γ˙ = U/L is the characteristic shear rate in the downstream channel for the contraction flow. The partial differential equations to be solved in dimensionless form are (i) mass conservation ∇ · v = 0,
(1)
(ii) momentum conservation Re
dv = ∇ ·T, dt
(iii) Oldroyd-B constitutive equation dτ T τ + De − ∇v τ − τ ∇v = 2η1 S, dt
(2)
(3)
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
155
where τ = T + pI − 2η0 S.
(4)
The velocity vector is v = vj ej , where ei (i = 1, 2, 3) are orthonormal cartesian basis vectors. The quantity T is the total (Cauchy) stress tensor, p is the pressure, and τ is the polymeric part of the extra-stress tensor. The operator d/dt is the material time derivative, I the identity tensor, and S the (symmetric) rate-of-deformation tensor defined as S = 21 (∇v + ∇v T ).
(5)
In Eqs. (3) and (5), the velocity gradient was defined as: (∇v)ij = ∂vj /∂xi , and ∇v T is its transpose. For future reference, it is convenient to define the vorticity tensor, W : W = 21 (∇v − ∇v T ).
(6)
3. The GTP formulation This section is divided into five subsections. In the first one, the general equations and the hypotheses to obtain the extra-stress model are presented and applied to the Oldroyd-B model, yielding the GTP[OldB] model. In the second subsection, a second law analysis is presented illustrating the capability of the formulation to predict anisotropic stress relaxation upon cessation of flow. In the third one, the GTP[OldB] model is applied to simple shear flow, illustrating its ability to predict non-zero second normal stress differences N2 . In the fourth subsection, the GTP formulation is extended to the PTT differential model, yielding the GTP[PTT] model. In the last subsection, a flow-type classifying parameter for plane flows is introduced. This parameter will be useful in the discussion of the results. 3.1. General equations The constitutive equation, Eq. (3), is rewritten in terms of S and W : dτ 1 2η1 =− τ+ S + (Sτ + τ S) − (W τ − τ W ). dt De De We now introduce the extra-stress traceless tensor Γ , defined as Γ = τ − 13 Iτ I , where Iτ = {τ } is the first invariant (the trace) of τ . Using this definition, Eq. (7) becomes 1 1 dIτ 1 1 η1 dΓ = − Γ + SΓ + Γ S − + Iτ I − (W Γ − Γ W ) + 2 + Iτ S. dt De 3 dt De De 3
(7)
(8)
(9)
An evolution equation for Iτ is obtained by taking the trace of Eq. (7): d{τ } 1 2η1 = − {τ } + {S} + {(Sτ + τ S)} − {(W τ − τ W )}. (10) dt De De This equation is significantly simplified as follows. Firstly we notice that {S} = 0 for isochoric motions. Then we multiply Eq. (8) by S and take the trace to obtain that {τ S} = {Γ S} + (Iτ /3){S} = {Γ S}.
156
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
Therefore, {(Sτ + τ S)} = 2{Γ S}. Moreover, the trace of W τ (and of τ W ) is zero, and hence we get the following equation for Iτ : dIτ 1 = − Iτ + 2{Γ S}. dt De Combining Eqs. (9) and (11), we obtain 1 2 η1 1 dΓ = − Γ + Γ + Γ S − {Γ S}I − (W Γ − Γ W ) + 2 + Iτ S. dt De 3 De 3
(11)
(12)
To obtain an explicit relation for Γ , we firstly need to obtain an algebraic form for dΓ /dt. In the turbulence model proposed by Gatski and Speziale [15], it was assumed that the extra-stress anisotropy tensor, Γ /Iτ , reaches an equilibrium state and therefore its time derivative is zero. This hypothesis was firstly used for viscoelastic fluids by Mompean et al. [21] and Mompean [19]. Following these ideas but being consistent with objectivity, we assume that the objective time derivative of Γ /Iτ is zero, or, equivalently, DΓ Γ dIτ , = Iτ dt Dt
(13)
where the operator D/Dt is any objective time derivative. This simplifying hypothesis preserves memory prediction capability, and is exactly verified for simple shear and elongational flows. The following general form for DΓ /Dt encompasses several objective time derivatives available in the literature [36]: DΓ dΓ = − c1 (SΓ + Γ S) − c2 (Γ W − W Γ ) − (1 − c2 ) (Γ Ω − ΩΓ ), Dt dt
(14)
where c1 and c2 are constants, and the tensor Ω is the rate-of-rotation of the eigenvectors of S(˙s i = Ω · s i , i = 1, 2 or 3, where s i is the ith eigenvector of S). This expression preserves symmetry for all values of c1 and c2 , and particular choices of these coefficients correspond to different objective time derivatives: the contravariant convected time derivative when c1 = c2 = 1; the covariant convected time derivative when c1 = −c2 = −1; the corotational or Jaumann derivative when c1 = 0 and c2 = 1; and the Harnoy derivative [11] when c1 = c2 = 0. This last one is based on Ω (DΓ /Dt = dΓ /dt − (Γ Ω − ΩΓ )). Eqs. (13) and (14) are now combined with Eq. (12) to yield Γ 2 η1 1 ¯ ¯ 0 = − 2{Γ S} + α1 (SΓ + Γ S) − {Γ S}I − α2 (W Γ − Γ W ) + 2 (15) + Iτ S, Iτ 3 De 3 ¯ = W − Ω is the relative-rate-of-rotation tensor, the parameters α1 = 1 − c1 and α2 = 1 − c2 where W have been introduced for compactness. It is clear that different values of α1 and α2 correspond to different choices of objective time derivative: α1 = α2 = 0 ⇒ contravariant convected time derivative; α1 = 2 and α2 = 0 ⇒ covariant convected time derivative; α1 = 1 and α2 = 0 ⇒ Jaumann derivative; and α1 = α2 = 1 ⇒ Harnoy derivative. ¯ , Iτ ). In Our task now is to use Eq. (15) to obtain an explicit expression for Γ of the form Γ (S, W the frame of Newtonian turbulent flows, Pope [26] has used a tensor polynomial expansion to obtain an explicit algebraic representation for the Reynolds stress tensor equation. With an analogous procedure, Gatski and Speziale [9] have extended this analysis to three-dimensional flows, obtaining an algebraic
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
157
stress model from a differential second-moment closure for turbulent flows. Since Γ is symmetric and traceless, the same ideas can be applied here, i.e. we can expand Γ as Γ =
N
βn T (n) ,
(16)
n=1
¯ , and the scalar where the basis tensors T (n) are known symmetric traceless tensor functions of S and W ¯ coefficients βn are functions of invariants involving Γ , S and W . In particular, choosing the three-term ¯ −W ¯ S and T (3) = S 2 − (1/3){S 2 }I , we get [14]: basis {T (1) , T (2) , T (3) }, where T (1) = S, T (2) = S W ¯ S} 6{Γ S 2 } {Γ W 1 2 {Γ S} 2 ¯ ¯ S+ S − {S }I . (17) SW − W S + Γ = 3 ¯ 2} {S 2 } {S 2 }2 {S 2 }{W For two-dimensional flows, this three-term basis gives an exact representation of Γ . For three-dimensional flows, N ≥ 5 provides an exact representation (see [14]). As a final step towards an explicit expression ¯ S} and {Γ S 2 } as functions of the invariants of for Γ , we need to write the scalar invariants {Γ S}, {Γ W ¯ . These are obtained from manipulations of Eq. (15). To obtain {Γ S}, we the kinematic tensors S and W multiply Eq. (15) by S and take the trace, to get 1 η1 1 2 2 ¯ 0 = − {Γ S} + α1 {Γ S } + α2 {Γ W S} + (18) + Iτ {S 2 }. Iτ De 3 The appropriate choice of the root of this equation gives (see [21]): (α1 )2 α1 1 η1 (α2 )2 2 ¯ 2 + − + Iτ Iτ {S 2 } + I {W }. {Γ S} = + De 2 3 3 2 τ
(19)
¯ S} and {Γ S 2 } are obtained similarly by multiplying Eq. (15) by S 2 The other two scalar invariants {Γ W ¯ and W S, respectively, and taking the trace. The result is ¯ S} = 1 α2 Iτ {W ¯ 2 }, {Γ W 2 and {Γ S 2 } = Iτ
1
α 2 1
−
1 3
{S 2 },
(20)
(21)
where the Cayley–Hamilton theorem has been applied to simplify some tensor products (see [21], Appendix A). Replacing the above expressions for the scalar invariants in Eqs. (11) and (17)), we get the two following expressions, which represent the GTP[OldB] model: S η1 (α1 )2 α1 1 (α2 )2 2 ¯ 2 Γ= 2 + − + Iτ Iτ {S 2 } + I {W } De 2 3 3 2 τ {S } α2 1 2 Iτ 2 ¯ ¯ , (22) (S W − W S) + (3α1 − 2) S − {S }I + 2 3 {S } 2 where Iτ is obtained from the solution of 1 dIτ η1 (α2 )2 2 ¯ 2 (α1 )2 α1 1 = − Iτ + 2 + − + Iτ Iτ {S 2 } + I {W }. dt De De 2 3 3 2 τ
(23)
158
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
It is shown shortly that, with appropriate choices of α1 and α2 , the GTP[OldB] model can predict non-zero second normal stress differences. 3.2. Second law analysis For isothermal processes of non-Newtonian fluids, the second law of thermodynamics reduces to (e.g. [2]): dA T · ∇v ≥ , (24) dt where A is the dimensionless Helmholtz free energy. The quantity on the left-hand side of the above expression is the so-called stress power. For incompressible inelastic materials, such as the Generalized Newtonian Liquid (GNL), whose thermodynamic state is defined by two intensive thermodynamic properties only, dA/dt = 0, since there is neither temperature nor density gradients. This means that these materials cannot store free energy and hence the stress power is always non-negative. For the Oldroyd-B fluid, inequality (24) can be re-written as dA 2η0 {S 2 } + {Γ S} ≥ . (25) dt To evaluate {Γ S}, we use Eq. (11), and the above inequality becomes 1 dA 1 dIτ 2 + Iτ > , (26) 2η0 {S } + 2 dt λ dt where the strict inequality has been adopted due to the existing irreversibilities (dissipation and relaxation). This inequality holds for both the Oldroyd-B fluid and the GTP[OldB] model, since Eq. (11) holds for both. For viscoelastic fluids, an important feature that the second law should carry is the capability of predicting relaxation upon cessation of flow. In this case, we have 1 dIτ 1 + Iτ = 0, (27) 2 dt λ and inequality (26) becomes dA . 0> dt The solution for Iτ in Eq. (27) is Iτ (t) = Iτ0 e−t/λ ,
(28)
(29)
where Iτ0 is the trace of τ at the time of cessation of the flow. By replacing Eq. (29) into Eq. (22) it can be seen that an anisotropic relaxation process of finite duration can be predicted by the GTP[OldB] model. More specifically, an inspection of Eq. (22) shows that the terms related to the basis tensors ¯ −W ¯ S) and (S 2 − (1/3){S 2 }I ), do not vanish, in general, upon a cessation of flow. From the (S W second law point of view, inequality (28) shows that during this process the stored Helmholtz free energy decreases (or, equivalently, entropy increases). As relaxation is an irreversible process, this behavior is physically expected.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
159
3.3. Predictions for simple shear flow For simple shear flow, the components of Γ can in general be written in terms of the normal stress differences. These are defined as N1 = TXX − TYY and N2 = TYY − TZZ , where X is the direction of the flow, Y is the direction of the velocity gradient and Z in the neutral direction. Accordingly, they are written in the present work as N1 = T11 − T33 and N2 = T33 − T22 . Therefore, for simple shear flow the components of Γ are 1 (2N1 + N2 ) 0 2η1 S13 3 [Γ ] = (30) 0 − 13 (N1 + 2N2 ) 0 , 2η1 S13
− 13 (N1 − N2 )
0
where S13 is the constant shear rate. The Oldroyd-B constitutive equation predicts null second normal stress difference in shear flows, i.e. N2 = 0. Therefore, for the Oldroyd-B model the components of Γ are: 2 N 0 2η1 S13 3 1 [Γ ] = 0 (31) − 13 N1 0 . 2η1 S13
− 13 N1
0
For the same flow, Eq. (22) gives for the GTP[OldB] model:
[Γ ] = η1 Iτ 2De
1 2 Iτ α1 + α2 − 2 3
0 +
α12 − α22 1 − α1 + 4 6
1/2
Iτ2
0
Iτ −α1 +
2 3
0
η1 Iτ + 2De
α12 − α22 1 − α1 + 4 6 0
1 2 Iτ α1 − α2 − 2 3
Iτ2
1/2 . (32)
Therefore, the prediction of the GTP[OldB] model for N1 is N1 = Γ11 − Γ33 = α2 Iτ .
(33)
It is interesting to note that, when α1 = α2 = 1, then the Oldroyd-B and GTP[OldB] models yield the same prediction for simple shear flow. A comparison between Eqs. (30) and (32) shows that, to keep the predictions for the shear stress unchanged, the following condition must be satisfied: 1 (α12 4
− α22 ) + 16 (1 − α1 ) = 0.
It can be shown that this condition restricts the range for the coefficient α2 to √ |α2 | ≥ 31 5. Another important result is that, if α1 and α2 are chosen to satisfy the following relation, N2 α2 − 2 = 0, 3α1 − 1 + 2 N1
(34)
(35)
(36)
160
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
then the GTP[OldB] model is able to predict non-vanishing second normal stress difference. For null second normal stress difference, α1 = α2 = 1 satisfy Eq. (36). The discussion of this section shows that the GTP can transform a differential model into another whose prediction capabilities may be better in some aspects than the ones of the original model. Moreover, different choices of α1 and α2 lead to different predictions of rheological functions. Hence, if the flow of interest is predominantly extensional, for example, then other guidelines (good prediction of extensional rheological properties) should be followed for a proper choice of α1 and α2 . 3.4. Transformation of other differential models: obtaining the GTP[PTT] model The same steps discussed in the previous subsections to apply the GTP on the Oldroyd-B constitutive equation can be applied to other constitutive differential models as well. In this section, we show how to apply the transformation to the PTT constitutive equation, thus obtaining the GTP[PTT] model. The PTT equation (see [24,23]) can be written in a general form as dτ T f ({τ })τ + De (37) − (∇v) τ − τ ∇v = 2η1 S. dt The function f depends on the elongational behavior, on the viscosity and on the relaxation time of the viscoelastic fluid. Using the previous definitions, Eq. (37) can be rewritten as dτ f 2η1 ¯ τ − τW ¯ ). =− τ+ S + (Sτ + τ S) − (W dt De De Taking the trace of Eq. (38), we can write: dIτ η1 (α2 )2 2 ¯ 2 f (α1 )2 α1 1 = − Iτ + 2 + − + Iτ Iτ {S 2 } + I {W }. dt De De 2 3 3 2 τ
(38)
(39)
The same procedure used to obtain Eq. (23), as shown in the last section, has been used here to obtain Eq. (39). It is important to note that, when applying the assumption made in Eq. (13), using now Eq. (39) for the PTT fluid, and rewriting Eq. (37) in terms of the symmetric traceless tensor Γ , the obtained equation will be exactly the same as Eq. (15). Then the final expression is identical to the one obtained for the GTP[OldB] model, namely, Eq. (22). The contribution of the function f of the PTT model appears only in the first right-hand side term of the Eq. (39). In other words, the GTP[PTT] model is thus given by Eqs. (22) and (39), which are to be solved simultaneously to give the stress field as a function of the flow kinematics. 3.5. Flow-type sensitive classifier As viscoelastic fluids exhibit different rheological behaviors when submitted to viscometric or extensional flows, it is important to identify regions where the fluid is being sheared, elongated or is moving as a rigid body. Generally, in a complex flow, all these kinds of motion are present and there are regions where the material is being exposed to a combination of the kinematics of these limiting cases. In order to aid our later discussion about the performance of the GTP in the contraction flow, it is interesting to have a means to locally classify the flow. To this end, the criterion proposed by Astarita [1]
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
161
based on the relative-rate-of-rotation tensor was adopted. This criterion is applicable to plane flows only [5,6,12,29–31,35–37], as is the case of the plane contraction. Astarita’s criterion is based on the parameter R, defined as ¯ } {W 2
R=−
{S } 2
=−
W¯ mn W¯ nm , Smn Snm
(40)
which will take values between 0 and +∞. It can be shown [1] that for extensional flow, R = 0 and as the flow approaches rigid body motion, R → ∞. Other types of plane flow lie in between, including simple shear flow, for which R = 1. This parameter is a measure of how much the material avoids stretching through rotation relative to the principal directions of S. In addition to classifying flows, this parameter has also been used in the construction of constitutive models [34,35,37]. In this work, this parameter was normalized to avoid numerical problems when R → ∞: D=
1−R . 1+R
(41)
The scalar D take values between −1 and +1 only. For simple shear flow, D = 0; for extensional flow, D = +1 and as the flow approaches a rigid body motion, D → −1. The parameter D is used here as an important tool to analyze which types of flow the GTP formulation gives the same results as the corresponding constitutive equations in a differential form. 4. Numerical method 4.1. Solution of governing equations Viscoelastic numerical simulations were carried out using two different sets of equations: (i) the Oldroyd-B fluid and (ii) the GTP[OldB] extra-stress model. To model the Oldroyd-B fluid flow the following system has to be solved: (a) mass conservation, Eq. (1), (b) momentum conservation, Eq. (2) and (c) the constitutive equation for the non-Newtonian extra-stress, Eq. (3). When using the GTP[OldB] extra-stress model the equations to be solved are: (a) mass conservation, Eq. (1), (b) momentum conservation, Eq. (2) and (c) the traceless extra-stress expression, Eq. (22) coupled with the equation for the trace of the non-Newtonian extra-stresses, Eq. (23). As it will be discussed in the section of results, in this work, we chose to compare the predictions for the stress field of different models for a given kinematic field. To this end, the Newtonian solution of a flow inside the 4:1 contraction at Reynolds number 0.01 has been chosen. In this manner, the conservation and the Newtonian constitutive equations were first solved, and then the constitutive models (Oldroyd-B and GTP[OldB]) were used to evaluate the stress field with the velocity field. This was one in order to analyze the effect of the proposed procedure on the stress field only, with no indirects effects of changes in the velocity field. To solve these differential equations, a finite-volume numerical method has been successfully employed. The spatial discretization is performed in a staggered grid [22]. The non-linear terms (convective flux) of the momentum equations, the advection terms of the constitutive equations for the Oldroyd-B differential model, and the advection terms for the trace of the extra-stress tensor (when using the GTP[OldB] formulation), are obtained by the second-order accurate scheme, quadratic upstream interpolation scheme
162
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
for convective kinematics (QUICK) proposed by Leonard [17]. The diffusion terms of the momentum equations are calculated with the second-order accurate centered difference scheme. The pressure and the normal stress components of the viscoelastic tensor are evaluated at the center of the control volumes. The velocities are staggered and evaluated at the center of the faces; the off-diagonal components of the viscoelastic tensor are attached to nodes at the mid-edges. In order to evaluate the shear extra-stress component (off-diagonal term), for example (τ13 (i, k)), its equation is solved at the center of the cell with the normal components, and then the values are extra-polated linearly at the corner of the cell. The decoupling procedure employed for the pressure is derived from the work of Harlow and Welch [10]. Details of this algorithm used to solve viscoelastic flows are well documented in Mompean and Deville [20], and are briefly mentioned next. From the momentum conservation equation, a discrete Poisson equation is obtained for the pressure by enforcing mass conservation implicitly. In the momentum equation, the diffusion terms due to the Newtonian stresses (solvent) and due to the non-Newtonian extra-stress components (viscoelastic fluid) are treated explicitly. The non-linear terms are also evaluated explicitly. The differential constitutive equation for the Oldroyd-B fluid is solved using an Euler explicit scheme. With such information, the right-hand side (B) of the linear system for the pressure is evaluated at time n, and the new pressure values are calculated at time n + 1 [A]pn+1 = B n .
(42)
The matrix [A] of the system, used to obtain the pressure (p) at n + 1, is symmetric and positive definite. It can be solved by a direct Cholesky factorization or by a preconditioned conjugate gradient method. This scheme represents a real evaluation in time, no pseudo-transient algorithm being used here. The steady-state solution is computed until the transient of the time-dependent process has vanished. 4.2. Mesh A uniformly spaced mesh with 3.000 nodes (75 nodes in the streamwise direction (x or x1 ) and 40 nodes in the normal direction (z or x3 )) has been used to solve the flow in the contraction. This mesh has 60 nodes upstream of the entrance (Le = 16H ) and 15 nodes downstream (L0 = 4H ), as shown in Fig. 1. In the region after the abrupt contraction, 10 nodes are used to discretize the outflow channel in the normal direction. 4.3. Boundary conditions For the momentum equation, the same boundary conditions are used whether using the Oldroyd-B differential constitutive equations or the GTP[OldB] model. These conditions are: (1) at the inlet a parabolic profile is imposed for the v1 velocity component; (2) at the walls, the no-slip condition is employed (v1 = 0; v3 = 0); (3) at the outlet of the domain a Neumann condition is used for all variables, except for the pressure, which is set to zero as a reference value. When solving the Oldroyd-B constitutive equations, the boundary conditions for the non-Newtonian extra-stress components are deduced from a plane Poiseuille flow, and hence at the inlet of the domain, the following conditions are imposed:
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
163
• τ11 = 2η1 De(∂v1 /∂x3 )2 , • τ33 = 0, • τ13 = η1 (∂v1 /∂x3 ). For the GTP[OldB] model, only one boundary condition has to be given (for the trace of the extra-stress tensor) at the inlet. For a fully developed Poiseuille flow, it was adopted, for the trace of τ , the same value that was used in the Oldroyd-B, and therefore, it takes the form: • Iτ = (1/α2 ) 2η1 De(∂v1 /∂x3 )2 . Accordingly, the above expression is employed as boundary condition at the inlet for the Iτ equation (Eq. (23)). 5. Results In order to better assess the differences between the predictions of the original and transformed constitutive equations, it is interesting to compare the corresponding stress fields for a fixed flow field. The velocity field chosen here was the Newtonian solution of the flow inside the 4:1 contraction. 5.1. The newtonian flow field For the numerical simulation of the flow of a Newtonian fluid through the 4:1 contraction, the final grid was selected after different grids were tested in order to ensure grid-independent results. The influence of the upstream channel length (Le ) on the results was verified; to allow the development of the flow in the entry region, a ratio of 4 between the upstream channel length (Le ) and the downstream channel length (L0 ) was used (see Fig. 1), which corresponds to a corner step position of x = 16. The Reynolds number based on the average exit velocity and on the downstream height was 0.01. The Newtonian velocity vector field with the streamlines is presented in Fig. 2. The dimensionless size of the recirculation, Xv = Lv /2Hu , is equal to 0.18, which is in good agreement with previously published results (e.g. [25]). The field of the second eigenvector of the rate-of-deformation tensor for the flow of a Newtonian fluid through the contraction is shown in Fig. 3. We note in this figure that the principal directions of the rate-of-deformation tensor are essentially constant in the fully-developed region of the upstream and downstream channels, as expected for a simple shear flow. The regions where significant changes in the principal directions of S occur are the entry region of the downstream channel (15 < x < 17; 0 < z < 1), and the vortex zone. These changes in principal directions are further illustrated with the aid of Fig. 4, where contours of the dimensionless angular velocity of the eigenvectors, scaled by (2H /U ), are shown. Only the region close to the contraction plane (14 ≤ z ≤ 18) is shown, because elsewhere the angular velocity of the eigenvectors of S is null. There are two neighboring regions where large angular velocities occur. The upstream region has positive angular velocities, while the downstream one has negative angular velocities. In the close vicinity of the contraction plane (x = 15.8 and z = 0.9), a positive value of order one is observed, while near the symmetry plane (x = 16.4 and z = 0.2), a negative value of −2.0 is found. 5.1.1. The flow-type classifier D Isobands of the flow-type classifier D are shown in Fig. 5 for the flow discussed above. The three limiting types of flow (extensional, shear and rigid body) are clearly identified, as discussed next. In the
164
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
Fig. 2. Velocity field and streamlines of the flow of a Newtonian fluid. Re = 0.01.
regions away from the contraction plane, the flow is fully developed and hence it is simple shear flow, which corresponds to D = 0. Just upstream of the contraction plane, a region of elongational flow is found, which is represented in Fig. 5 by the clear contours, where D 1. In this region, the streamlines converge as the fluid flows into the downstream channel, and the material particles are stretched.
Fig. 3. Field of the second eigendirection of the rate-of-deformation tensor, for the flow of a Newtonian fluid. Re = 0.01.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
165
Fig. 4. Angular velocity of the eigenvectors of the rate-of-deformation tensor near the contraction plane for the flow of a Newtonian fluid. Re = 0.01.
Concerning elliptical flows (flows “between” shear and solid body motion), there are three regions where the values of D become negative: (a) just after the sharp corner; (b) near the corner vortex, by the contraction plane; and (c) just before the extensional region at the entrance of the downstream channel.
Fig. 5. Isobands of the D flow-type classifier.
166
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
The occurrence of the first two regions are easy to understand. At the corner there are vortices which are expected to rotate at rates larger than the local deformation rate; near the lip, the fluid makes a sharp turn to enter the downstream channel, and the rate-of-rotation is again larger than the local deformation rate. The third region cannot be seen as a rotation motion but as a motion close to rigid body motion (plug flow). At first glance, it may seem strange that such region of rigid body motion (D → −1) appears between a region of (expected) shear (D = 0) and of (expected) extension (D = 1). However, the existence of such region becomes clear when we examine separately the behaviors of the numerator and denominator of Eq. (40) (for details, see [36]). The deformation rate in the shear region vanishes at the centerline and has its maximum at the wall. Therefore, a particle that travels near the centerline comes from a low deformation rate region, but, as it approaches the contraction plane, the eigenvectors of S have to change their direction, from ± π/4 with respect to the x-axis (shear flow), to zero radians (extension). This rotation of the eigenvectors in this region can also be observed in Fig. 3. Thus, the existence of rigid body motion in this region indicates that, in that region, the rates of rotation of the eigenvector directions are larger than the rates of deformation. It can be seen that the opposite is true for regions far from the centerline, where transition from shear to extension occurs smoothly as far as D is concerned. 5.2. Range of Deborah number for calculations with the Oldroyd-B model In order define the range of admissible values for the Deborah number when using the Oldroyd-B constitutive equation, we performed a simplified analysis which gives approximate expressions for the non-Newtonian extra-stresses. In this analysis, we assumed that the flow is a motion with constant stretch history (MWCSH, e.g. [13]) and hence we can neglect the advection terms from the constitutive Eq. (3). Therefore, the components of the non-Newtonian extra-stress become τ11 =
2η1 S11 + 2Dev1,3 τ13 , 1 − 2DeS11
(43)
τ13 =
2η1 [S13 + DeS11 (v3,1 − v1,3 )] , 2 1 − 4De2 (S11 + v3,1 v1,3 )
(44)
τ33 =
−2η1 S11 + 2Dev3,1 τ13 , 1 + 2DeS11
(45)
where v3,1 = ∂v3 /∂x1 and v1,3 = ∂v1 /∂x3 . Because the denominators of Eqs. (43) and (44) can go to zero depending on the value of De, τ11 and τ13 have a singular behavior. This imposes a limit that can be expressed in terms of the Deborah number, for a given velocity field. The limit imposed by the non-Newtonian extra-stress τ11 is Demax < 1/2S11 . The limit related to τ13 is not as simple to calculate, but to avoid the singularities we moni2 tored the limits imposed by Eqs. (43) and (44) by plotting contours of 1 − 4De2 (S11 + v3,1 v1,3 ) and 1 − 2DeS11 to observe possible changes of sign. For the flow studied here (Re = 0.01), we calculated the velocity gradients and then the limiting Deborah numbers, which turned out to be De = 1.13 and De = 0.38 corresponding to Eqs. (43) and (44), respectively. Accordingly, to avoid these singularity problems of the Oldroyd-B constitutive equation, De ≤ 0.38 for all the cases studied in the research reported here.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
167
5.3. Performance of the GTP[OldB] model in the contraction flow 5.3.1. Prediction of non-zero second normal stress differences in shear flow As mentioned before, the GTP[OldB] model is capable of predicting non-vanishing second normal stress differences in simple shear flow (N2 = 0), depending upon the choices of the parameters α1 and α2 . In fact, we can choose a ratio of N2 to N1 that is representative to a given material, and then obtain the values of α1 and α2 with the aid of Eq. (36). For illustration purposes, we have chosen a value for N2 /N1 and obtained the corresponding values of α1 and α2 . Our choice of N2 /N1 was based on the result of an analysis of the rod-climbing phenomenon for a second order fluid found in Joseph [15] which yields −0.25N1 ≤ N2 ≤ 0 for vanishingly small shear rates. This range is in agreement with data for many polymeric liquids, and, accordingly, we have assumed here N2 = −0.1N1 . With this assumption, we get α1 = 0.92 and α2 = 0.95. Fig. 6 shows the distribution of N2 along the normal axis for position P3 (where the flow is viscometric). Because Ψ2 = N2 /γ˙ 2 is assumed to be constant, this curve is a parabola, as expected. 5.3.2. Prediction of the stress field and comparisons with other models Using the flow field obtained for a Newtonian fluid, we have evaluated stress fields as predicted by the following constitutive models: • the Oldroyd-B model, with η0 = 1/9 and η1 = 8/9, • the non-objective AESM-W model,
Fig. 6. Second normal stress difference for simple shear flow (N2 ) as predicted by the GTP[OldB] model (α1 = 0.92, α2 = 0.95).
168
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
• the GTP[OldB] model with α1 = α2 = 1 (with these values of α1 and α2 , the GTP[OldB] model reproduces the steady material function predictions of the Oldroyd-B model for simple shear flow), • the GTP[OldB] model with α1 = 0.92 and α2 = 0.95 (with these values of α1 and α2 , the GTP[OldB] model predicts, in shear flow, second normal stress differences such that N2 = −0.1N1 ). Three representative axial positions, namely, P1, P2 and P3, are marked with dashed lines in Fig. 5. These correspond respectively to the dimensionless axial coordinate values of x = 12, x = 15.5 and x = 19.6, and were chosen for analyzing the components of Γ because different types of kinematics occur along each one. The first position, P1, is a position along all kinds of motion exist, and the changes of flow-type occur smoothly across the channel. Position P2 lies where very strong gradients of D occur. It encompasses the vortex region, the region near the lip and the extensional region close to the centerline. The third one, P3, is a position where the flow is simple shear all along. Position P3 was also used to validate the results obtained by the numerical code, because the solution for viscometric flow can also be obtained analytically with the aid of Eqs. (31) and (32). Figs. 7 and 8 show the predictions given by the different models for the Γ11 and Γ13 tensor components, respectively. For both components, it is excellent the agreement of the analytically obtained Oldroyd-B predictions with the numerically obtained Oldroyd-B predictions and also with the predictions of the GTP[OldB] model with α1 = α2 = 1, as expected. Moreover, the AESM-W model also agrees well with the theoretical model, as also expected, because at this position the vorticity and relative-rate-of-rotation tensors are equal due to the absence of rotation of the eigenvectors of S.
Fig. 7. Numerical results for Γ11 as a function of the transversal coordinate z at axial position P3 (x = 19.6) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). The continuous line pertains to the predictions of the Oldroyd-B model obtained analytically for simple shear flow (th-Old-B). Re = 0.01 and De = 0.36.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
169
Fig. 8. Numerical results for Γ13 as a function of the transversal coordinate z at axial position P3 (x = 19.6) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). The continuous line pertains to the predictions of the Oldroyd-B model obtained analytically for simple shear flow (th-Old-B). Re = 0.01 and De = 0.36.
Fig. 9. Numerical results for Γ11 as a function of the transversal coordinate z at axial position P1 (x = 12) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). Re = 0.01 and De = 0.36.
170
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
Fig. 10. Numerical results for Γ13 as a function of the transversal coordinate z at axial position P1 (x = 12) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). Re = 0.01 and De = 0.36.
Fig. 11. Numerical results for Γ11 as a function of the transversal coordinate z at axial position P2 (x = 15.5) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). Re = 0.01 and De = 0.36.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
171
Fig. 12. Numerical results for Γ13 as a function of the transversal coordinate z at axial position P2 (x = 15.5) as predicted by the Oldroyd-B, AESM-W and two GTP[OldB] models (α1 = 1, α2 = 1 and α1 = 0.92, α2 = 0.95). Re = 0.01 and De = 0.36.
In Fig. 7, it is also seen that the curve for Γ11 given by the GTP[OldB] model with α1 = 0.92 and α2 = 0.95 is slightly displaced from the curves predicted by the other models. This is also expected, because the GTP[OldB] model predicts Γ11 = (N1 /2α2 )(α1 + α2 − 2/3) (see Eqs. (32) and (33)) for simple shear flow. Then, Γ11 = (2/3)N1 for α1 = α2 = 1, whereas for α1 = 0.92 and α2 = 0.95, Γ11 = (1.9/3)N1 . Figs. 9 and 10 show Γ11 and Γ13 , respectively, for the four models at position P1. The predictions of the four models are in good agreement for both components, but especially for Γ11 . For Γ13 (Fig. 10), there is little difference among the predictions of the two GTP[OldB] models and the AESM-W model, a particularly good agreement being observed close to the wall and near the centerline. The results predicted by the Oldroyd-B model for Γ13 differ a little from the others in the region of extensional flow. At position P2 (Figs. 11 and 12), the predictions for Γ11 and Γ13 of the four models are again in good agreement, although not as good as in positions P3 and P1. As also observed in position P1, there is not much difference among the predictions of the two GTP[OldB] models and the AESM-W model. The results predicted by the Oldroyd-B model again depart from the others in the region of extensional flow.
6. Concluding remarks This work presents a general procedure to obtain a model for the extra-stress tensor, starting from a differential constitutive model. It reduces the set of differential equations for the components of the non-Newtonian extra-stress to one differential equation for its trace only, plus an explicit equation for the traceless stress tensor. This procedure is inspired in the literature of turbulent flows.
172
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
The procedure was applied to the Oldroyd-B constitutive model with a three-term basis of tensors that gives an exact representation for the traceless tensor in two-dimensional flows. The resulting constitutive model, namely, the GTP[OldB] model, was employed to the steady, two-dimensional planar flow through a 4:1 abrupt contraction. The Newtonian velocity field was used to compute the stress tensor for the AESM-W model and for two other GTP[OldB] models produced by particular choices of the general objective time derivative. The first is a kind of corotational derivative based on the rotation rate of the eigenvectors of the rate-of-deformation tensor, and leads to the same stress tensor of the original model for shear flow. The other was defined by keeping the non-zero rheological functions of the Oldroyd-B models in shear flow, but including second normal stress difference prediction capability. A flow-type indicator, D, was used as a tool in the analysis of the differences and similarities of the models. The predictions of the GTP[OldB] model are in qualitative agreement with the ones of the Oldroyd-B model, and, even quantitatively, the agreement is quite reasonable. The possibility of prediction of non-zero N2 is an interesting result. Although the procedure is essentially a method for simplifying a differential model (in the sense described in Section 3), it leads to an enhancement in the prediction capability, because another rheological property is captured by the transformed model. This result suggests that the GTP model based on the White–Metzner model (the GTP[WM] model) would probably be able to predict non-zero and shear rate dependent second normal stress coefficients, in addition to the prediction capabilities of the original model. A thermodynamic analysis has shown that the GTP[OldB] model is able to predict Helmholtz free energy storage. This was exemplified with the aid of the sudden cessation of flow situation. Other constitutive models that perform better in extensional flows can be transformed to obtain other GTP models, as it was exemplified for the PTT model. By allowing the parameters α1 and α2 to be ¯ , realistic predictions can be obtained not only for one of the functions of the second invariants of S and W limiting cases (shear flow or extensional flow) but for both at the same time, and better quality predictions are also expected for complex flows. However, it is important to emphasize that other differential models may lead to difficulties in the application of the GTP, especially in the process of obtaining explicit expressions for the coefficients βn that appear in Eq. (16). Acknowledgements The authors wish to acknowledge the financial support from Coordenação de Aperfeiçoamento de Pessoal de N´ıvel Superior (CAPES, Brazil) and Comité Français d’Evaluation de la Coopération Universitaire avec le Brésil (COFECUB, France) under contract 354/01. The second author thanks the Fundação de Amparo Pesquisa do Estado do Rio de Janeiro (FAPERJ) for its support. The authors would also like to thank Dr. L. Thais (EUDIL) and Dr. A. Rivens (Lille III) for their constructive remarks during the preparation of the manuscript. References [1] G. Astarita, Objective and generally applicable criteria for flow classification, J. Non-Newtonian Fluid Mech. 6 (1979) 69–76. [2] G. Astarita, M.E. Mackay, The generalised engineering Bernoulli equation (gebe), and the first and second laws of thermodynamics for viscoelastic fluids, J. Rheol. 40 (1986) 335–346.
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
173
[3] F.T.P. Baaijens, Numerical analysis of start-up planar and axisymmetric contraction flows using multi-mode differential constitutive models, J. Non-Newtonian Fluid Mech. 48 (1993) 147–180. [4] D.V. Boger, Viscoelastic flows through contractions, Annu. Rev. Fluid Mech. 19 (1987) 157–182. [5] P.O. Brunn, E. Ryssel, The w-D fluid: general theory with special emphasis on stationary two-dimensional flows, Continuum Mech. Thermodyn. 9 (1997) 73–82. [6] P.O. Brunn, E. Ryssel, The Giesekus fluid in w-D form for steady two-dimensional flows. Part I. Theory, Rheol. Acta 38 (1999) 415–422. [7] R.E. Evans, K. Walters, Flow characteristics associated with abrupt changes in geometry in the case of highly elastic liquids, J. Non-Newtonian Fluid Mech. 20 (1986) 11–29. [8] R.E. Evans, K. Walters, Further remarks on the lip-vortex mechanism of vortex enhancement in planar-contraction flows, J. Non-Newtonian Fluid Mech. 32 (1989) 95–105. [9] T.B. Gatski, C.G. Speziale, On explicit algebraic stress models for complex turbulent flows, J. Fluid Mech. 254 (1993) 59–78. [10] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182–2189. [11] A. Harnoy, Stress relaxation effect in elastico-viscous lubricants in gears and rollers, J. Fluid Mech. 76 (3) (1976) 501–517. [12] R.R. Huilgol, Comments on objective and generally applicable criteria for flow classification, by G. Astarita, J. Non-Newtonian Fluid Mech. 7 (1980) 91–95. [13] R.R. Huilgol, N. Phan-Thien, Fluid Mechanics of Viscoelasticity, No. 6 in Rheology Series, Elsevier, Amsterdam, 1997. [14] T. Jongen, T. Gatski, General explicit algebraic stress relations and best approximation for three-dimensional flows, Int. J. Eng. Sci. 36 (1998) 739–763. [15] D.D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Applied Mathematical Sciences, Springer-Verlag, New York, 1990. [16] R. Keunings, On the high Weissenberg number problem, J. Non-Newtonian Fluid Mech. 20 (1986) 209–226. [17] B.P. Leonard, A stable accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Eng. 19 (1979) 59–98. [18] H. Matallah, P. Townsend, M.F. Webster, Recovery and stress-splitting schemes for viscoelastic flows, J. Non-Newtonian Fluid Mech. 75 (1998) 139–166. [19] G. Mompean, On predicting abrupt contraction flows with differential and algebraic viscoelastic models, Comput. Fluids 31 (2002) 935–956. [20] G. Mompean, M. Deville, Unsteady finite-volume simulation of Oldroyd-B fluid through a three-dimensional planar contraction, J. Non-Newtonian Fluid Mech. 72 (1997) 253–279. [21] G. Mompean, T. Jongen, T. Gatski, M. Deville, On algebraic extra-stress models for the simulation of viscoelastic flows, J. Non-Newtonian Fluid Mech. 79 (1998) 261–281. [22] V.S. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, New York, 1980. [23] N. Phan-Thien, A nonlinear network viscoelastic model, J. Rheol. 22 (1978) 259–283. [24] N. Phan-Thien, R.I. Tanner, A new constitutive equation derived from network theory, J. Non-Newtonian Fluid Mech. 2 (1977) 353–365. [25] T.N. Phillips, A.J. Williams, Viscoelastic flow through a planar contraction using a semi-Lagrangian finite-volume method, J. Non-Newtonian Fluid Mech. 87 (1999) 215–246. [26] S.B. Pope, A more general effective-viscosity hypothesis, J. Fluid Mech. 72 (1975) 331–340. [27] B. Purnode, M.J. Crochet, Flows of polymer solutions through contractions. Part 1. Flows of polyacrylamide solutions through planar contractions, J. Non-Newtonian Fluid Mech. 65 (1996) 269–289. [28] R.S. Rivlin, The relation between the flow of non-Newtonian fluids and turbulent Newtonian fluids, Q. Appl. Math. 15 (1957) 212–215. [29] E. Ryssel, P.O. Brunn, Comparison of a quasi-Newtonian fluid with a viscoelastic fluid in planar contraction flow, J. Non-Newtonian Fluid Mech. 86 (1999a) 309–335. [30] E. Ryssel, P.O. Brunn, Flow of a quasi-Newtonian fluid through a planar contraction, J. Non-Newtonian Fluid Mech. 85 (1999) 11–27. [31] E. Ryssel, P.O. Brunn, The Giesekus fluid in w-D form for steady two-dimensional flows. Part II. Numerical simulation, Rheol. Acta 38 (1999) 423–436. [32] G.P. Sasmal, A finite-volume approach for calculation of viscoelastic flow through an abrupt axisymmetric contraction, J. Non-Newtonian Fluid Mech. 56 (1995) 15–47.
174
G. Mompean et al. / J. Non-Newtonian Fluid Mech. 111 (2003) 151–174
[33] T. Sato, S. Richardson, Explicit numerical simulation of time-dependent viscoelastic flow problems by finite-element/finite-volume method, J. Non-Newtonian Fluid Mech. 51 (1994) 249–275. [34] P.R. Schunk, L.E. Scriven, Constitutive equation for modeling mixed extension and shear in polymer solution processing, J. Rheol. 34 (7) (1990) 1085–1117. [35] P.R. Souza Mendes, M. Padmanabhan, L.E. Scriven, C.W. Macosko, Inelastic constitutive equations for complex flows, Rheol. Acta 34 (1995) 209–214. [36] R.L. Thompson, Performance of a new constitutive equation for non-Newtonian liquids (in Portuguese), Ph.D. thesis, Pontif´ıcia Universidade Católica, PUC-Rio, Brazil, 2001. [37] R.L. Thompson, R.P. Souza Mendes, M.F. Naccache, A new constitutive equation and its performance in contraction flows, J. Non-Newtonian Fluid Mech. 86 (1999) 375–388. [38] S.C. Xue, N. Phan-Thien, R.I. Tanner, Three-dimensional numerical simulations of viscoelastic flows through planar contractions, J. Non-Newtonian Fluid Mech. 74 (1998) 195–245.
