13
Computers and Chemical Engineering 28 (2004) 2363–2375 An adaptive remeshing strategy for shear-thinning fluid flow simulations A. Fortin a,, F. Bertrand b , M. Fortin a , É. Chamberland a , P. E. Boulanger-Nadeau a , A. El. Maliki a , N. Najeh a a GIREF, Départment de Mathématiques et de Statistique, Université Laval, Québec, Canada G1K 7P4 b Départment de Génie Chimique, École Polytechnique de Montréal, Montréal, Canada, H3C 3A7 Received 20 May 2004; accepted 2 June 2004 Abstract In the last few years, a team effort has allowed the development of a very general adaptive finite element solution procedure which can be applied to a large variety of problems. In this paper, this complete strategy is presented and applied to the solution of three-dimensional shear-thinning fluid flow problems. Carreau–Yasuda and power-law models are used to model shear-thinning viscosity effects. © 2004 Elsevier Ltd. All rights reserved. Keywords: Shear-thinning fluids; Iterative methods; Error estimation; Adaptive remeshing; Finite element method 1. Introduction Accurate computations of more and more complex flow problems are necessary in the design of almost all poly- mer processes. Two-dimensional computations are no longer sufficient for most industrial problems. Three-dimensional problems are, however, much more difficult and require im- portant computer time and huge memory if care is not taken. Direct solvers (Gaussian elimination) are not suitable and iterative solvers must be implemented. It is also important to assess the accuracy of the re- sults. In two-dimensional problems, it is still possible to refine the mesh (almost) everywhere in the domain to test mesh independence. This is generally not possible for three-dimensional computations. Adaptive remeshing strategies based on error estimation as described in Habashi and co-workers (Ait Ali Yahia et al., 2002; Dompierre, Vallet, Bourgault, Fortin, & Habashi, 2002; Habashi, Dompierre, Bourgault, Fortin, & Vallet, 1998; Habashi et al., 2000) and Hecht and Mohammadi (1997) provide a tool to refine the mesh only where necessary, allowing very accurate solutions at minimum cost. Large linear or non linear systems requires efficient iter- ative solvers. The literature on iterative solvers is abundant Corresponding author. Tel.: +1-418-656-3489; fax: +1-418-656-3404. E-mail address: [email protected] (A. Fortin). and the reader is referred to Bertrand and Tanguy (2002) and Guénette, Fortin, Marcotte, and Labbé (2004) for a com- plete review of iterative methods for the Stokes problem. In this last paper, a very efficient Stokes solver suitable for quadratic elements was developed allowing computations on very fine meshes. The same solver will be used in this work. An efficient strategy for the three-dimensional simulation of shear-thinning fluid flows such as those encountered in polymer processes comprises: a second-order mixed finite element; a preconditioned iterative method for solving linear and non linear systems; an appropriate adaptive remeshing strategy. The aim of this work is to present such a procedure. In order to illustrate the potential of the proposed strategy, four problems will be considered, each with its own particular- ities: the three-dimensional Poiseuille flow, a 18 to 1 con- traction flow, the flow of a sphere falling in a tube and the flow in a static mixer. 2. Mathematical formulation In many applications, polymer melts can be modelled without viscoelasticity effects. In these cases, the Cauchy stress tensor has the form: 0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2004.06.007

An adaptive remeshing strategy for shear-thinning fluid flow simulations

Embed Size (px)

Citation preview

Page 1: An adaptive remeshing strategy for shear-thinning fluid flow simulations

Computers and Chemical Engineering 28 (2004) 2363–2375

An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortina,∗, F. Bertrandb, M. Fortina, É. Chamberlanda, P. E. Boulanger-Nadeaua,A. El. Maliki a, N. Najeha

a GIREF, Départment de Mathématiques et de Statistique, Université Laval, Québec, Canada G1K 7P4b Départment de Génie Chimique, École Polytechnique de Montréal, Montréal, Canada, H3C 3A7

Received 20 May 2004; accepted 2 June 2004

Abstract

In the last few years, a team effort has allowed the development of a very general adaptive finite element solution procedure which canbe applied to a large variety of problems. In this paper, this complete strategy is presented and applied to the solution of three-dimensionalshear-thinning fluid flow problems. Carreau–Yasuda and power-law models are used to model shear-thinning viscosity effects.© 2004 Elsevier Ltd. All rights reserved.

Keywords:Shear-thinning fluids; Iterative methods; Error estimation; Adaptive remeshing; Finite element method

1. Introduction

Accurate computations of more and more complex flowproblems are necessary in the design of almost all poly-mer processes. Two-dimensional computations are no longersufficient for most industrial problems. Three-dimensionalproblems are, however, much more difficult and require im-portant computer time and huge memory if care is not taken.Direct solvers (Gaussian elimination) are not suitable anditerative solvers must be implemented.

It is also important to assess the accuracy of the re-sults. In two-dimensional problems, it is still possible torefine the mesh (almost) everywhere in the domain totest mesh independence. This is generally not possiblefor three-dimensional computations. Adaptive remeshingstrategies based on error estimation as described in Habashiand co-workers (Ait Ali Yahia et al., 2002; Dompierre,Vallet, Bourgault, Fortin, & Habashi, 2002; Habashi,Dompierre, Bourgault, Fortin, & Vallet, 1998; Habashiet al., 2000) andHecht and Mohammadi (1997)provide atool to refine the mesh only where necessary, allowing veryaccurate solutions at minimum cost.

Large linear or non linear systems requires efficient iter-ative solvers. The literature on iterative solvers is abundant

∗ Corresponding author. Tel.:+1-418-656-3489;fax: +1-418-656-3404.

E-mail address:[email protected] (A. Fortin).

and the reader is referred toBertrand and Tanguy (2002)andGuénette, Fortin, Marcotte, and Labbé (2004)for a com-plete review of iterative methods for the Stokes problem.In this last paper, a very efficient Stokes solver suitable forquadratic elements was developed allowing computationson very fine meshes. The same solver will be used in thiswork.

An efficient strategy for the three-dimensional simulationof shear-thinning fluid flows such as those encountered inpolymer processes comprises:

• a second-order mixed finite element;• a preconditioned iterative method for solving linear and

non linear systems;• an appropriate adaptive remeshing strategy.

The aim of this work is to present such a procedure. Inorder to illustrate the potential of the proposed strategy, fourproblems will be considered, each with its own particular-ities: the three-dimensional Poiseuille flow, a 18 to 1 con-traction flow, the flow of a sphere falling in a tube and theflow in a static mixer.

2. Mathematical formulation

In many applications, polymer melts can be modelledwithout viscoelasticity effects. In these cases, the Cauchystress tensor has the form:

0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved.doi:10.1016/j.compchemeng.2004.06.007

Page 2: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2364 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

σ = −pI + τ,

wherep is the pressure and the extra stress tensorτ can bewritten as:

τ = 2η(|γ(u)|)γ(u), (1)

expressing the viscosity dependence on the rate of straintensorγ(u) of the velocity fieldu = (u1, u2, u3) defined by:

(γ(u))ij = 1

2

(∂ui

∂xj+ ∂uj

∂xi

), (2)

or more precisely on its second invariant|γ(u)|:|γ(u)|2 = 2γ(u) : γ(u) = 2(γ(u))ij (γ(u))ij .

The viscosity can be modelled in different ways. Thepower-law model can be expressed as:

η = η0|γ(u)|n−1, (3)

wheren is the pseudoplasticity (or power-law) index andη0 is the consistency index. The power-law is often usedto model the viscosity behavior of molten metals. Its mainweakness is that the viscosity tends towards∞ for smallshear rates. The Carreau equation:

η = η0(1+ λ2|γ(u)|2)(n−1)/2, (4)

or its generalization, the Carreau–Yasuda model, given by:

η = η0(1+ λm|γ(u)|m)(n−1)/m, (5)

alleviate this problem and are well adapted to polymer melts.In the following, the development of the numerical method-ology will be done with the Carreau–Yasuda model, but thepower-law or Carreau models would follow the same lines.

Since polymers are incompressible and inertia is ne-glected, the non linear Stokes-like problem can be written:−∇ · (2η(|γ(u)|)γ(u)− pI) = 0,

∇ · u = 0.(6)

To simplify the notation, homogeneous Dirichlet bound-ary conditions will be assumed (u = (0, 0, 0) on∂Ω) but thegeneralization to non homogeneous Dirichlet or Neumannconditions is straightforward.

The associated variational formulation is obtained by mul-tiplying the above system by test functionsw ∈ V andq∈ QwhereV andQ are appropriate functional spaces. The nonlinear Stokes-like problem takes the classical form:

∫Ω(2η0(1+ λm|γ(u)|m)(n−1)/mγ(u) :

γ(w)− p∇ · w)dv = 0,∫Ωq∇ · u dv = 0.

(7)

In this work, the discretization of the above system wasperformed using the second-order (O(h2)) Taylor–Hood el-ement (seeBrezzi and Fortin, 1991) illustrated inFig. 1.

Fig. 1. Taylor–Hood element.

The velocity components are thus approximated by piece-wise quadratic polynomials while the pressure is piecewiselinear.

Newton’s method can be applied to solve the non linearsystem (7). Starting from a first solution approximation (u0,p0), a correction vector (u, p) is computed as the solutionof the linear system:

au0(δu,w)− ∫Ωδp∇ · w dv = −∫

ΩR1((u0, p0),w)dv,∫

Ωq∇ · δu dv = −∫

ΩR2((u0, p0), q)dv,

(8)

where the bilinear formau0(δu,w) is:

au0(δu,w)

=∫Ω

(2η0(1+ λm|γ(u0)|m)(n−1)/mγ(δu) : γ(w))dv

+∫Ω

(α(1+ λm|γ(u0)|m)(n−m−1)/m(γ(u0) : γ(δu))

× (γ(u0) : γ(w)))dv,

and whereα = 4η0(n − 1)λm|γ(u0)|m−2 and the residualcomponentsRi((u0, p0), w) are given by:

R1((u0, p0),w) = ∫Ω

2η0(1+ λn|γ(u0)|n)(m−1)/nγ(u0) :

γ(w)dv− ∫Ωp0∇ · w dv,

R2((u0, p0), q) =∫Ωq∇ · u0 dv.

Once a solution of this system has been obtained, theinitial estimate of the solution is updated:

u0← u0+ δu and p0← p0+ δp.

The procedure is repeated until the norm of the correction(u, p) is small enough.

3. Iterative solution of Stokes-like problems

The linear system (8) must be solved at each iterationof Newton’s method. In three-dimensional applications,this system can be huge, involving hundreds of thousandsequations. An efficient iterative solver is thus extremelyimportant. A three-dimensional non linear Stokes solverwas presented inGuénette et al. (2004). This solver wasdesigned for highly viscous non Newtonian fluid flow prob-lems and was used in this work. A brief description will

Page 3: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2365

now be given and the reader is referred toGuénette et al.(2004)for more details.

The Stokes-like problem for each Newton’s iteration canbe written in matrix form as:[A(uk) Bt

B 0

] [δu

δp

]= −

[R1(uk, pk)

R2(uk, pk)

]. (9)

The solver is based on the Arrow-Hurwicz scheme (seeFortin & Glowinski, 1983), which is obtained by artificiallyinserting time derivatives of bothu and p in the Stokesproblem. This gives:

∂u

∂t= f − Au− Btp,

∂p

∂t= Bu,

which in turn, using an explicit–implicit scheme, can bewritten in matrix form:

Mu(uk+1− uk) = f − Auk − Btpk (explicit),

Mp(pk+1− pk) = Buk+1 (implicit),

whereMu andMp are symmetric positive definite mass ma-trices. Usinguk+1 = uk + u and pk+1 = pk + p, theabove system can finally be written as:[Mu 0

B −Mp

] [δu

δp

]=

[f

0

]−

[A Bt

B 0

] [uk

pk

]

= −[R1(uk, pk)

R2(uk, pk)

].

In practice, this scheme does not converge very quickly.One could use the left hand side matrix as a triangular pre-conditioner for the system (9). Unfortunately, this wouldlead to a very expensive scheme due to the size of matrixMu. Consequently, a better choice would be to consider thefollowing family of triangular preconditioners:[Pu 0

B −Pp

] [A(uk) Bt

B 0

] [δu

δp

]

= −[Pu 0

B −Pp

] [R1(uk, pk)

R2(uk, pk)

]. (10)

Several choices are possible for the preconditioning ma-trices Pu and Pp ranging from the expensive choicesPu

= A(u) andPp = Mp, to the cheapestPu = diag(A) andPp

= diag(Mp) where diag(E) stands for the diagonal part ofthe matrixE. In this work, the preconditionerPu consistsin a few (less than 10) SSOR iterations with the matrixAwhile Pp = Mp since the size of the pressure–mass matrixis comparatively small.

This choice of preconditioner yields a symmetric but in-definite matrix so that conjugate gradient or conjugate resid-ual methods cannot be used (seeSaad, 1996for details

on Krylov subspace methods). Therefore, the BiCGSTABmethod proposed byvan der Vorst (1992)was chosen. Thelost of symmetry is largely compensated by the nice conver-gence properties of this solver and by the fact that it extendsto the case of the non symmetric Navier–Stokes equations.

4. Brief description of the adaptive method

This method has been abundantly described in the litera-ture and the reader is referred toHabashi et al. (1998, 2000)andHecht and Mohammadi (1997)for a complete presen-tation. Only a very brief description of the method will begiven.

4.1. Error estimator based on a metric

Starting from a linear approximationuh of a solutionu,the error can be written as:

e = u− uh,

whereh is the element length. In the one-dimensional case,it is well-known from elementary numerical analysis thatthe maximum error on an element verifies:

emax= h2

8

d2u

dx2(ξ),

for someξ in the element. Following the terminology ofHabashi et al. (2000), the equidistribution of the error isachieved on a mesh if:

h2

8

∣∣∣∣d2u

dx2(xi)

∣∣∣∣ = C, (11)

at every nodexi of the domain for some prescribed toleranceC. Such a mesh is then called optimized. To achieve this goal,second-order derivatives are approximated at every nodexiof the domain andEq. (11)determines an element lengthmap. A new mesh satisfying as much as possible the elementlength map is then produced and a new solution is computed.

The same analysis is possible for three-dimensional prob-lems when working on element edges. An optimized meshmust satisfy:

l2

8

∣∣∣∣d2u

dv2(xi)

∣∣∣∣ = C,

wherel is the length of an edge starting at nodexi andv isa unit vector tangent to this edge. The tangential derivativeis defined as:

d2u

dv2(xi) = vtH(xi)v.

where H is the Hessian matrix of functionu containingsecond-order derivatives and is obviously unknown. It can,however, be reconstructed using the method described in

Page 4: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2366 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

Section 4.2. The optimized mesh must then satisfy:

l2

8vtH(xi)v = C. (12)

If nodesxi andxj define an edge of lengthl, then the unittangential vector is:

v = xi − xj

l,

and from (12), the optimized mesh must satisfy:

(xj − xi)tH(xi)(xj − xi) = C.

If H is positive definite, this relation defines a metric:

||xi − xj||H =√(xj − xi)tH(xi)(xj − xi) = L, (13)

where L is a prescribed target value for the edge lengththroughout the computational domain. The length of theedges is now related to the error. Consequently, when a valueof L is fixed, the error level is also fixed. If the Hessian isnot positive definite, one uses its absolute value.

The new mesh is then obtained from the previous one bytrying to equidistribute the error on the mesh or more pre-cisely on the edges of the mesh. This can be achieved byproducing a mesh with edges satisfying relation (13). Thealgorithm proceeds by trying to equilibrate the error by edgeswapping and node displacement, and control its level byedge refinement or node suppression. This must of coursebe done in the respect of the geometrical boundaries of thecomputational domain, which is perhaps the main difficultyin three-dimensional problems. The method is restricted tolinear approximations since, otherwise, the error is no longerrelated to the Hessian. When using a quadratic approxima-tion, as in this work, one can neglect the values of the so-lution on element edges. This definitely leads to a loss ofinformation in the adaptive process.

Eq. (13)requires approximations of second-order deriva-tives of the (unknown) solutionu. In the case of a linear ap-proximation, it is not possible to replaceu by the numericalsolution uh whose second-order derivatives vanish. In thenext section, an accurate approximation of these derivativesat mesh nodes is constructed.

4.2. Approximation of first- and second-order derivatives

The method will be explained for a two-dimensional meshbut its extension to the three-dimensional case is straight-forward. Starting from a standardC0 finite element solu-tion uh of any degree in elementK, an approximation of the

first- and second-order derivatives at the element vertices isneeded. Sinceuh is not differentiable at the element vertices,one possibility is to fit a local quadratic interpolation on ev-ery vertexP∗ = (x∗1, x

∗2) of the mesh. To do so, one has to

define a patch of elements adjacent to the vertex (seeFig. 2).The patch usually consists of adjacent elementsK but it isalso possible (and in some cases necessary) to include ele-ments not directly connected to the vertex (see second casein Fig. 2). Over this patch, the solution is approximated bya quadratic (truncated) Taylor series:

uapp(x1, x2)

= uh(x∗1, x∗2)+ u

appx1 (x∗1, x

∗2)(x1− x∗1)

+ uappx2 (x∗1, x

∗2)(x2− x∗2)+ u

appx1x1(x

∗1, x∗2)

12(x1− x∗1)

2

+ uappx1x2(x

∗1, x∗2)(x1− x∗1)(x2− x∗2)

+ uappx2x2(x

∗1, x∗2)

12(x2− x∗2)

2,

where onlyuh(x∗1, x∗2) is known. If the patch of elements

adjacent to vertex(x∗1, x∗2) containsN nodesP i = (xi1, x

i2),

then the idea is to determine the best values (in a least squaresense) of first and second-order derivatives, that is valuesminimizing:

uapp(xi1, xi2)− uh(x

i1, x

i2), i = 1, 2, 3, . . . , N,

which leads to the following minimization problem:

miny||By − b||2, (14)

where

y =

uappx1 (x∗1, x

∗2)

uappx2 (x∗1, x

∗2)

uappx1x1(x

∗1, x∗2)

uappx1x2(x

∗1, x∗2)

uappx2x2(x

∗1, x∗2)

, b =

uh(x11, x

12)− uh(x

∗1, x∗2)

uh(x21, x

22)− uh(x

∗1, x∗2)

uh(x31, x

32)− uh(x

∗1, x∗2)

...

uh(xN1 , xN2 )− uh(x

∗1, x∗2)

,

andB is theN × 5 matrix (N × 9 in three dimensions):

B =

(x11− x∗1)

(x21 − x∗1)

(x31− x∗1)

...

(xN1 − x∗1)

(x12− x∗2)

(x22 − x∗2)

(x32− x∗2)

...

(xN2 − x∗2)

12(x

11− x∗1)

2

12(x

21 − x∗1)

2

12(x

31− x∗1)

2

...

12(x

N1 − x∗1)

2

(x11− x∗1)(x

12− x∗2)

(x21 − x∗1)(x

22 − x∗2)

(x31− x∗1)(x

32− x∗2)

...

(xN1 − x∗1)(xN2 − x∗2)

12(x

12− x∗2)

2

12(x

22 − x∗2)

2

12(x

32− x∗2)

2

...

12(x

N2 − x∗2)

2

.

Once this procedure has been done for each node of thetriangulation, it is quite easy to interpolate the first andsecond-order derivatives in an element. For example, whenmoving a node, derivatives must be evaluated at its newlocation. First, one has to determine the element contain-ing this new position. In this element, values of the first-and second-order derivatives on the vertices are now known

Page 5: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2367

Fig. 2. Nodes adjacent toP∗.

from the solution of the minimization problem (14). Con-sequently, a linear interpolation can be used to approximatethe values of second-order derivatives at the new node po-sition. Similarly, a linear interpolation can also be used todetermine first-order derivatives. For better results, however,and since the values of second-order derivatives at the ele-ment nodes are already known, a third-order Hermite inter-polation can be used to interpolate first-order derivatives.

5. Numerical simulations

Four problems were considered to assess the efficiencyand the accuracy of the proposed strategy. One of the diffi-culties is the appropriate choice of the value ofL which isstrongly linked to the total number of elements. Too smalla value ofL would quickly result in an exceedingly largenumber of elements. Consequently, the value ofL in Eq. (13)was chosen in order to obtain accurate solutions but also tocontrol the total number of elements. Moreover, in all simu-lations, the variable of adaptation (denotedu in Section 4.1)will be the velocity norm:

||u||2 =√u2x + u2

y + u2z .

This means that the adaptive method will try to detect thevariations in the velocity norm and adapt the mesh accord-ingly. The Hessian matrix of the velocity norm will thendefine the new metric (13) on which the whole adaptationprocess is based.

5.1. Poiseuille flow in a cylinder

In this first example, the objective is to show how meshadaptation yields meshes that depend upon the rheology ofthe polymer. The Poiseuille flow in a cylinder provides avery simple illustration of this behavior. A pressure dropis imposed at both ends of a circular cylinder and no slipboundary conditions are imposed along its walls.

In this example, a power-law fluid is considered and theadaptive strategy is applied to two values of the consistencyindex: n = 1 (Newtonian) andn = 0.4. It is well knownthat whenn = 1, the velocity profile is parabolic but whenn decreases, the velocity profile flattens out in the middle ofthe cylinder resulting in steep gradients close to the cylinderwall. This is clearly illustrated inFig. 3 where iso-values

Fig. 3. Axial velocity ux in a cross-section of the cylinder:n = 1 (top)and n = 0.4 (bottom).

of the axial velocity on a plane perpendicular to the flowdirection are shown. The influence of the power-law indexn is easily seen and the adapted meshes should reflect thedifferences in the solutions.

5.1.1. Newtonian caseLet us first consider the Newtonian case (n = 1). The

value ofL which controls the edge lengths and thus the er-ror level was set arbitrarily (by trial and error) so that thenumber of elements remains smaller than 10 000. The pro-cedure was started with a uniform mesh ( 6000 elements)given inFig. 4. The adaptive strategy was then applied twice(using the same value ofL) resulting in two other mesheswith 8221 and 9813 elements. On the last mesh, the averageedge length was close enough to the target valueL.

The flow being uniaxial, the mesh adaptation techniquetends to produce elements with very large edge lengths alongthe flow direction as can be seen inFig. 5. It is also usefulto look at the exit section of the cylinder to appreciate thedifferences between the initial and adapted meshes.Fig. 6

Page 6: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2368 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

X

YX

YVU

Z

Z

Fig. 4. Initial mesh.

shows the initial and adapted meshes at the exit section of thecylinder where the velocity profile is completely developed.The adapted meshes provide uniformly distributed elementsin the cross-section in accordance with the parabolic velocityprofile (seeFig. 3) as expected.

5.1.2. Non Newtonian case: n= 0.4It is expected in this case that the mesh should be more

refined near the cylinder boundary while preserving the el-ement elongation in the axial direction. Here again, twoadapted meshes were computed starting from the same initialmesh and value ofL as in the Newtonian case. The first ma-jor difference is the number of elements (17 752 and 17 313)which is significantly higher due to the presence of a steepvariation of the axial velocity along the walls (see the sec-ond picture inFig. 3). Fig. 7 illustrates the adapted mesheswhich are refined along the walls and coarser in the middleof the cylinder where the solution is nearly constant. This isa nice example that shows how the mesh, the accuracy of thesolution and the fluid rheology are intricately interrelated.

5.2. 18 to 1 contraction

In this problem, the flow of a Carreau–Yasuda fluid ina die with a 18 to 1 contraction is considered. The fluid

Y

Z

VU

XYZ12

3X

Fig. 5. Element elongation in the flow direction for the 9813 elementmesh.

parameters are:η0 = 39 275 Pa s, n = 0.0797,

λ = 0.4 m, m = 0.539,

and the dimensions of the geometry are the following:Length of the reservoir= 28 mm,Length of the die= 25 mm (−28≤ x ≤ 25),Height of the reservoir= 18.0 mm (−9 ≤ y ≤ 9),Height of the die= 1.0 mm,Width of reservoir and die exit= 10.0 mm (−5 ≤ z≤ 5).

At the reservoir entrance, a velocity profile with flow rate510 mm3/s is imposed:

u = (518 (1.0− (1

9y)2)(1.0− (1

5z)2),0,0),

while no slip is imposed on the reservoir and die walls. Theinitial (almost uniform) mesh is illustrated inFig. 8 andcontains 105 120 elements. This a typical situation wherea grid was obtained using a mesh generator with limitedadaptive functionalities. This mesh is probably appropriateto give satisfactory results but at a large computational price.Too many elements are present in the reservoir while in thedie, more elements should be added in the neighborhoodof the die walls. The adapted mesh ( 33 800 elements) hasexactly the right characteristics, as also illustrated inFig. 8.The adaptation variable was once again the velocity normand only one adapted mesh was built.

The diminution of the number of elements should notresult in a loss of accuracy.Fig. 9 presents a comparisonof the axial velocity profiles on the planey = 0 at the exitsectionx= 25 (the horizontal line passing right in the middleof the exit section). The steep variation of the velocity closeto the walls is again perfectly reproduced.

The other important variable in polymer processing isthe pressure which is illustrated inFig. 10for both mesheson the liney = z = 0 passing through the middle of thedomain in the flow direction. Though the pressure is notexplicitly taken into account in the adaptive process, it isnevertheless accurately predicted. As can be seen in the lasttwo figures, no accuracy was lost when decreasing from105 120 to 33 800 elements, the main difference being thecomputational time.

5.3. Drag on a falling sphere in a tube

In a recent paper,Missirlis, Assimacopoulos, Mitsoulis,and Chhabra (2001)presented an analysis of wall effects

Page 7: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2369

Fig. 6. Initial and adapted meshes at cylinder exit:n = 1.

on the drag on a sphere falling in a circular cylinder.The axisymmetry of the geometry (seeFig. 11) allowedtwo-dimensional computations on very fine meshes. It wasclearly shown that the shear thinning index (in a power-law

Fig. 7. Initial and adapted meshes at cylinder exit:n = 0.4.

viscosity model) and the ratioRc/Rs of cylinder to sphereradii are critical to the computed value of the drag. Walleffects become more important as the ratio of cylinder tosphere radiiRc/Rs becomes smaller.

The same problem is now considered but in a fullthree-dimensional simulation, which represents a muchgreater challenge. The objective is to show that evenwhen starting in the worst possible conditions (a verycoarse mesh and a three-dimensional approximation of atwo-dimensional problem) mesh adaptation provides veryaccurate results similar to those obtained on very finetwo-dimensional meshes.

A uniform velocity field is imposed on the top surfaceand on the cylinder walls. No slip was imposed on the sur-face of the sphere andux = uy = (σ·n)z = 0 at the bottomof the cylinder. A power-law fluid (η0 = 1) was consideredand drag values were computed for different values of thepower-law index and for different cylinder to sphere ratios.The computed values will now be compared to those ob-tained inMissirlis et al. (2001).

Page 8: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2370 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

VU

X

Y

Z

X

Y

Z

VU

X

Y

Z

X

Y

Z

Fig. 8. Meshes for the 18 to 1 contraction.

The drag force on the sphere is defined as:

D=∫Γs

· n ds =∫Γs

(−pI + τ) · n ds

=∫Γs

(−pn+ τ · n)ds,

Fig. 9. Axial velocity profile on the liney = 0 and x = 25 at the exitsection.

Fig. 10. Pressure on the liney = z= 0 for uniform and adapted meshes.

whereΓ s is the surface of the sphere,n its unit normal vectorandτ is given by relations (1–3). Due to the symmetry ofthe problem, the first two components of the drag vectorshould vanish so that:

D = (0,0,Dz).

In the following, results will be presented for power-lawindices between 1.0 and 0.4. The strategy is always the same:adapt the mesh using the velocity norm until the computeddrag is satisfactory in the sense that the first two componentsare close to 0 (<10−2). The drag is then computed for thenext value ofn starting from the last computed mesh. Thisstrategy allows us to benefit from all the previous adaptationwork.

5.3.1. Test case #1: Rc/Rs = 2In the first case, the ratioRc/Rs is set to 2 and a Newto-

nian fluid (power-law indexn= 1) is used. InMissirlis et al.(2001), a valueDz = 56.05 was found in a two-dimensional

Fig. 11. Drag on a sphere: geometry.

Page 9: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2371

Fig. 12. Rc/Rs = 2, n = 1: adapted meshes.

axisymmetric computation with very fine meshes. Inthis three-dimensional computation, three meshes illus-trated in Fig. 12 were necessary to stabilize the valueof the drag. The results are summarized in the followingtable.

Test case #1:Rc/Rs= 2,n= 1 (target value fromMissirliset al., 2001): −56.05)

Mesh Nb ofelements

Nb ofdof

Drag values

Mesh 1 1826 9194 (−14.62,−9.273,−48.86)Mesh 2 9067 41781 (−0.344, 1.3296,−57.40)Mesh 3 52730 235222 (0.0099, 0.0022,−55.55)

The initial mesh is very coarse and clearly inappropriatesince the sphere itself is poorly discretized. The numericalsolution is not symmetrical and the first two components ofthe drag vector, which should be 0, are respectively,−14.62and−9.27. The computed drag is quite far from the targetvalue 56.05.

As the mesh is refined and adapted to the computed so-lution, values of the first two components of the drag goto 0 as expected. The drag coefficientDz reaches−55.55which is within 0.9% of the value obtained inMissirliset al. (2001). The meshes together with iso-values of thevelocity norm and isobars are illustrated inFigs. 12–14on the planex = 0 passing through the center of thesphere.

The following tables provide the drag coefficients forn= 0.8, 0.6 and 0.4. As can be easily seen, the adaptation

Fig. 13. Rc/Rs = 2, n = 1: iso-values of the velocity norm.

work done forn = 1 proved beneficial, the results beingacceptable with the first mesh. When adaptation is per-formed at n = 0.8, the mesh is only slightly rearrangedand the resulting mesh proves sufficient for all other val-ues ofn. All computed drag values are within 1% of thetarget values.Fig. 15 shows the computed solution forn = 0.4.

Fig. 14. Rc/Rs = 2, n = 1: isobars.

Page 10: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2372 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

Fig. 15. Rc/Rs = 2, n = 0.4: mesh, iso-values of velocity norm andpressure.

Test case #1:Rc/Rs = 2, n = 0.8 (target value fromMissirlis et al., 2001): −41.17)

Mesh Nb ofelements

Nb ofdof

Drag values

Mesh 1 52730 235222 (0.0084,+0.0010,−40.82)Mesh 2 52833 235531 (0.0078,+0.0024,−40.82)

Test case #1:Rc/Rs = 2, n = 0.6 (target value fromMissirlis et al., 2001: −30.18)

Mesh Nb ofelements

Nb ofdof

Drag values

Mesh 1 52833 235531 (0.0063, 0.0011,−29.97)

Test case #1:Rc/Rs = 2, n = 0.4 (target value fromMissirlis et al., 2001: −22.07)

Mesh Nb ofelements

Nb ofdof

Drag values

Mesh 1 52833 235531 (0.0045, 0.0002,−21.89)

It may be surprising that a mesh adapted for the solution atn= 0.8 provides acceptable drag values for smaller values ofthe power-law index. This behavior was not expected but itcan be explained in the following way. The adaptation workperformed forn = 1 reaches two objectives: the respect ofthe geometry of the sphere and the accurate computation ofthe drag. Once this is done, very little work is needed forsmaller values ofn.

Fig. 16.Rc/Rs= 4, n= 0.4: mesh, iso-values of velocity norm and isobars.

5.3.2. Test case #2: Rc/Rs = 4Increasing the ratioRc/Rs has the obvious effect of re-

ducing significantly the value of the drag. The computationswere performed following the same lines as in the previouscase, starting from a coarse mesh and with a power-law in-dex of 1. The results are summarized in the following table:

Test case #2:Rc/Rs = 4

n Computed drag force Target valuesfrom Missirliset al. (2001)

1.0 (−0.0027,−0.0009,−9.28) −9.35880.8 (−0.0023,−0.0008,−7.62) −7.66980.6 (−0.0027,−0.0006,−6.18) −6.23590.4 (−0.0016,−0.0001,−4.94) −4.9977

Here again, the agreement is excellent and the computeddrag coefficients values are within 1% of those obtained inMissirlis et al. (2001)on the finest mesh. A typical mesh,iso-values of velocity norm and the isobars are illustratedin Fig. 16 for n = 0.4. The most important conclusion isthat what was done here in an axisymmetric geometry, canbe easily performed in an arbitrary fully three-dimensionalgeometry without modifications.

5.4. Flow in a static mixer

The last example is the flow in a low pressure drop (LPD)motionless mixer. The LPD mixer consists of five mixerelements each consisting of two semielliptical plates dis-tributed in a tubular housing. Incoming streamlines are split

Page 11: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2373

Fig. 17. Static mixer: initial and three adapted meshes.

and diverted by the plates thus provoking mixing; the largerthe number of plates, the better the mixing. This is an ex-ample of a complex geometry of industrial relevance.

In this example, the computations were started with aquasi uniform 25 828 element mesh. At the first attempt andsince the flow-field was rather complex, mesh adaptationproduced a mesh with more than 200 000 elements. It is thenclear that a much better solution was obtained. However, itwas decided to restrict the total number elements to fewerthan 50 000, i.e. twice the initial number of elements. Theidea was to show that mesh adaptation can produce a bettersolution by placing the elements where needed but withoutunduly increasing the computational burden. The value ofLin relation (13) was then reduced accordingly.

For this example, only the Newtonian case will be pre-sented. A parabolic profile was imposed at the inlet while noslip boundary conditions were imposed on the cylinder andthe mixer elements (the plates). Mesh adaptation producedthe three meshes presented inFig. 17.

The pressure drop:

/P = Pin − Pout,

in the static mixer will now be examined and compared withthe following empirical relation provided by the manufac-turer:

/P = 351η0Vne

D,

whereη0 is the viscosity,V the mean velocity,ne the num-ber of mixer elements andD is the diameter of the tubularhousing. The computations were performed in a non dimen-sional geometry and the following values were used:η0= 1,V= 0.5,ne= 5 andD = 1. With these values, the empiricalrelation gives/P = 877.5.

As can be seen in the following table, the pressure is notvery sensible to mesh adaptation in this problem. The fourmeshes yielded essentially the same pressure profile on thez-axis (in the flow direction) as illustrated inFig. 18. Thepressure drops by successive steps reveal the presence of theplates. The differences in the total pressure drop are indeedsmall and the computed values are close to 877.5, the errorbeing smaller than 2%.

Static mixer

Mesh Number of elements Pressure drop

1 25828 8802 31877 8553 47636 8634 49472 860

It is thus extremely difficult to discriminate between thecomputed solutions and the pressure field is of no help. Thevelocity field should reveal some differences but it is usu-ally hard to get valuable information from a picture ofthree-dimensional velocity vectors. Some interesting fea-tures of the flow can, however, be obtained by looking atthe velocity field on a plane passing through the axis of themixer. A close-up view is presented inFig. 19 where thevelocity vectors are shown between two element mixers forthree different meshes (#1, #2 and #4 in the above table). Onthe initial mesh (top) the presence of the plates is not welltaken into account since the mesh is too coarse to detectthe flow along the plates. On the other hand, as the mesh is

0 5-5-10-15-20 10 15 20

0

100

-100

200

300

400

500

600

700

800

900

z

P

Mesh 1Mesh 2Mesh 3Mesh 4

Fig. 18. Pressure on thez-axis for the four meshes.

Page 12: An adaptive remeshing strategy for shear-thinning fluid flow simulations

2374 A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375

Fig. 19. Velocity field between two mixing elements: meshes #1, #2 and #4.

progressively adapted (middle and bottom pictures inFig.19), details of the flow close to the mixer elements start toappear. An animation of the particle trajectories (not pre-sented) on the last adapted mesh has shown the presence ofrecirculation zones behind some of the plates, which havea strong influence on the residence time distribution. Theserecirculations were not observed on the initial (coarser)mesh since the velocity field was not rich enough close tothe mixer elements.

6. Conclusions

A complete mesh adaptation procedure allowing thecomputation of accurate solutions of shear-thinning fluidflow problems was presented. It combines a second-orderelement, an efficient iterative solver and a good adaptiveremeshing strategy. The variety of applications presentedin this work clearly demonstrates the full potential of theproposed methodology.

Page 13: An adaptive remeshing strategy for shear-thinning fluid flow simulations

A. Fortin et al. / Computers and Chemical Engineering 28 (2004) 2363–2375 2375

References

Ait Ali Yahia, D., Baruzzi, G., Habashi, W. G., Fortin, M., Dom-pierre, J., & Vallet, M.-G. (2002). Anisotropic mesh adaptation:Towards user-independent, mesh-independent and solver-independentCFD. Part II: Structured grids.Int. J. Numer. Methods Fluids, 39, 657–673.

Bertrand, F., & Tanguy, P. A. (2002). Krylov-based Uzawa algorithmsfor the solution of the Stokes equations using discontinuous-pressuretetrahedral finite elements.J. Comp. Phys., 181, 617–638.

Brezzi, F., & Fortin, M. (1991).Mixed and hybrid finite element methods.Springer-Verlag.

Dompierre, J., Vallet, M.-G., Bourgault, Y., Fortin, M., & Habashi, W. G.(2002). Anisotropic mesh adaptation: Towards user-independent andsolver-independent CFD. Part III: Unstructured meshes.Int. J. Numer.Methods Fluids, 39, 675–702.

Fortin, M., & Glowinski, R. (1983).Augmented Lagrangian methods:Applications to the numerical solution of boundary-value problems,volume 15 of Studies in mathematics and its applications. North-Holland.

Guénette, R., Fortin, A., Marcotte, J.-P., & Labbé, J. (2004). Iterativesolvers for quadratic discretizations of the Stokes problem.Int. J.Numer. Methods Fluids, 44(7), 695–720.

Habashi, W. G., Dompierre, J., Bourgault, Y., Ait Ali Yahia, D., Fortin,M., & Vallet, M.-G. (2000). Anisotropic mesh adaptation: Towardsuser-independent, mesh-independent and solver-independent CFD. PartI: General principles.Int. J. Numer. Methods Fluids, 32, 725–744.

Habashi, W. G., Dompierre, J., Bourgault, Y., Fortin, M., & Vallet, M.-G.(1998). Certifiable computational fluid dynamics through mesh opti-mization.AIAA J., 36, 703–711.

Hecht, F., & Mohammadi, B. (1997). Mesh adaptation by metric controlfor multi-scale phenomena and turbulence.AIAA, 97-0859.

Missirlis, K. A., Assimacopoulos, D., Mitsoulis, E., & Chhabra, R. P.(2001). Wall effects for motion of spheres in power-law fluids.J. NonNewtonian Fluid Mech., 96, 459–471.

Saad, Y. (1996).Iterative methods for sparse linear systems. PWS Pub-lishing Company.

van der Vorst, H. (1992). Bi-CGSTAB: A fast and smoothly convergingvariant of Bi-CG for the solution of non-symmetric linear system.SIAM J. Sci. Stat. Comput., 12, 631–644.