Operator splitting algorithm for viscoelastic flow and numerical analysis for the flow around a sphere in a tube

Journalof Nou-Newtouian Fluid Mechanics ELSEVIER

J. Non-Newtonian Fluid Mech., 63 (1996) 121 140

,

Operator splitting algorithm for viscoelastic flow and numerical analysis for the flow around a sphere in a tube X.-L. Luo CSIRO Division o f Mathematics and Statistics, Locked Bag 17, North Ryde N S W 2113, Australia Received 17 July 1995: in revised form 16 October 1995

Abstract An operator splitting finite element algorithm for transient as well as steady viscoelastic flows has been implemented using a special high level language Fasttalk. The momentum equation is solved by operator splitting with a preconditioned conjugate gradient method, and the constitutive equation is integrated by an implicit scheme with a consistent streamline upwind Petrov-Galerkin (SUPG) method. The high efficiency of the present decoupled scheme has enabled us to do mesh refinement study using more than 105 degrees of freedom on a workstation computer for the viscoelastic flow around a sphere in a tube, and a speed of CPU zc N L3 has been delivered by the current code (N is the total number of unknowns). Steady state solutions were obtained for Deborah numbers De < 2.8 for the Upper-Convected Maxwell fluid. Comparison of the present results with those found in the literature (all due to coupled methods) shows excellent agreement on the predicted drag coefficient K for De _< 1.4. At higher De numbers (De _> 1.6) the present decoupled computation predicted no levelling off in K, while the other coupled computations predicted a levelling off in K followed by an earlier loss of convergence.

Keyword~': Decoupled method; Maxwell model; Operator splitting; Viscoelastic flows

1. Introduction

Numerical algorithms for viscoelastic flows are either of the decoupled type [1-6] or the mixed type [7-11]. In a decoupled method the momentum equation is solved separately from the constitutive equation, which is much less expensive than its coupled counterpart. A properly constructed decoupled scheme should in general have a larger radius of convergence than a coupled scheme which often relies on Newton type iterations. The decoupled method is also a natural choice for integral viscoelastic models. The numerical accuracy and computer efficiency 0377-0257/'96/$15.00 (© 1996- Elsevier Science B.V. All rights reserved SSDI 0377-0257(95)01411 -X

122

X.-L. Luo / J . Non-Newtonian Fluid Mech. 63 (1996) 121 140

of a decoupled algorithm depend on three aspects: (1) solving the constitutive equation given a velocity field; (2) solving the Navier Stokes (NS) equation given a viscoelastic stress field; (3) how the two equations are decoupled. High efficiency is particularly important for transient problems where the equations have to be solved many times, and for mesh refinement study and almost all three-dimensional problems where the problem size is very large. In some earlier work Tanner and co-workers [4,5] adopted a simple steady state Picard iteration scheme for the decoupling of constitutive equation and momentum equation, which combined with quite coarse grids and no upwinding of any kind, did not work very well by today's standard, reaching a Deborah (De) number value only about 0.5 for the standard sphere-in-tube problem. In the present work, a true time stepping scheme has been used which is not only more stable but also more robust than the steady type Picard iteration scheme, since no arbitrary relaxation factors are needed in the transient scheme. There are two difficulties traditionally associated with solving NS equations for incompressible flow--one is the incompressibility constraint or the coupling of pressure and velocity, the other is the nonlinear convection. For viscoelastic flows often the second difficulty is replaced by a material nonlinearity. The operator splitting (OS) method [12-14] treats nonlinearity and incompressibility in Navier-Stokes equations separately by 'splitting' operators with a fractional time stepping strategy. The method is suitable for both steady and transient problems and can readily be extended to include extra equations describing additional physical effects such as turbulence and viscoelasticity. For example, the author [15] has recently successfully applied an operator splitting algorithm to turbulent flow modelling using a nonlinear k - e model which, interestingly, has some resemblance to the Upper-Convected Maxwell (UCM) model [16,17]. The efficiency and accuracy of the OS algorithm mainly come from breaking the problem of solving NS equations into much simpler sub-problems and applying special methods to solve the sub-problems accurately and efficiently. The preconditioned conjugate gradient (PCG) algorithm of Glowinski [12] for quasi-Stokes (QS) problems is especially suitable for decoupled algorithms of viscoelastic flows, where often the Reynolds number is very small but viscoelastic effects dominate. The momentum equation in a decoupled scheme describes a quasi-Stokes problem with a complex non-Newtonian stress contribution, therefore it is essential for any decoupled scheme to have an accurate, efficient and robust QS problem solver. The PCG scheme proves to be faster and more robust than some other conventional algorithms such as the penalty method, and the former often gives better pressure solutions. The other essential part of a decoupled algorithm is solving the constitutive equation, the differential form of which is of hyperbolic type, often with large convection terms. The streamline upwind Petrov-Galerkin (SUPG) method [18] has been shown to be a numerically stable technique for solving linear hyperbolic equations. The effect of inconsistent streamline upwinding (SU) on numerical solutions has been a topic of much interest [7,19,20]. The SU method was first successfully applied to viscoelastic flow by Marchal and Crochet [7] in a mixed formulation. Subsequently Luo and Tanner [19] applied the same SU method in a decoupled formulation and achieved similar success, but they cautioned about the possible errors the SU may bring to the solution. Tanner and Jin [20] later on carried out a detailed comparative study on several numerical schemes including the SU method, and established quantitatively the inaccuracy of SU when applied to a one-dimensional case with equal or

X.-L. Luo / J . Non-Newtonian Fluid Mech. 63 (1996) 121-140

123

varying mesh sizes. They also found that the SUBS (SU on Both Sides) is as effective as SU but without the errors associated with SU. The flow of U C M fluid past a sphere placed on the centerline of a cylindrical tube has received much attention in the last few years, partly due to the acceptance of this problem as a benchmark for numerical simulations at the 5th Workshop on Numerical Methods in Non-Newtonian Flows [21]. Controversy has arisen from numerical predictions as well as experiments. The existence of a steady solution with increasing Deborah number and the cause of numerical divergence at limiting De number are the important questions surrounding the calculations of U C M flow around a sphere. The calculations for the U C M fluid proved to be very difficult at the 5th Workshop [21] when the sphere-in-tube problem was first attempted, with all calculations failing to converge below D e = 1.0. Much progress has been made since then. Up to now the most successful calculations are all due to the mixed formulations [9,11,22 24]. Using the explicitly elliptic m o m e n t u m equation (EEME) formulation [8] or elastic viscous split stress (EVSS) formulation [10], Lunsmann et al. [22] first obtained steady state solutions for the UCM fluid for De values up to 1.6 with the standard sphere-in-tube geometry (ratio of sphere radius to cylinder radius is a/R = 0.5). Good convergence was demonstrated in their work by mesh refinement study. Jin et al. [9] extended the EEME formulation to a class of constitutive equations of the Maxwell type, and obtained convergent results for UCM and PTT models up to De values of 2.0 and 4.5 respectively. Crochet and Legat [23] and Khomami [11] used high-order elements in their mixed formulations and achieved similar success as the EEME formulation. More recently Fan and Crochet [24] developed a highorder EVSS method and obtained convergence up to Dc number of 2.1 for the standard sphere-in-tube problem. The calculations of the sphere-in-tube problem at high De numbers with any of the existing coupled schemes proved to be extremely expensive. So far the maximum number of degrees of freedom (DOF) the coupled schemes could afford is about 50 000. Although the decoupled schemes are much less expensive to use than the couled ones, tile ability of the decoupled schemes to get good results at high De numbers for the sphere-in-tube problem has not been demonstrated before now. The first purpose of this paper is to describe the implementation of OS/SUPG, a decoupied finite element scheme for viscoelastic flow incorporating operator splitting and SUPG. The second purpose is to apply the new scheme to the analysis of the standard benchmark sphere-in-tube problem, demonstrating the accuracy and efficiency of the OS/SUPG l/\~rmulation as a decoupled scheme. After a section on governing equations, the OS/SUPG formulation will be described in Section 3, and this will be followed by a section on some numerical details including solution strategy, convergence criteria and mesh refinement. Numerical results will be given in Section 5. A very fine unstructured triangle mesh with over l0 s D O F has been used to explore the limiting De number question. Some discussions will be given on the loss of convergence to steady state. Concluding remarks are given in the final section.

124

X.-L. Luo /J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

2. Governing equations

The m o t i o n of an incompressible viscoelastic fluid is governed by the following equations: Ou ~ - + ( u . V ) u + Vp = V. z,

(1)

V. u = 0,

(2)

where u is velocity, p is pressure and ~ is the extra-stress tensor. Note in the m o m e n t u m equation we have retained the convection term for generality. For the U C M fluid the extra-stress r satisfies the following constitutive equation: At + ~r = 2r/0d,

(3)

where 2 is the relaxation time, r/0 is the constant viscosity, d is the rate of deformation tensor, and Az/At is the upper-convective time derivative A~ & - ~ = ~ + u • Vz - Lz - zL T,

(4)

where L is the velocity gradient tensor. In Eq. (1) there is no explicit viscous stress term, which would prevent us from using most existing finite element algorithms designed for elliptic type N a v i e r - S t o k e s equations. However, we can reformulate the system in terms of extra viscoelastic stress s: (5)

s - • - 2q0d, and the m o m e n t u m equation and the constitutive equation now become 0u & + (u • V)u - ~/oAU + Vp = V" s, 2

+u'Vs-Ls-sLV

+s=_22qo

(6) +u'Vd-Ld-dLm

.

(7)

3. O S / S U P G formulation for viscoelastic flow

The task in our transient scheme is to find (s" + l, u n + ~, p , + ~) at the (n + 1)th time step from (s ~, u", p") at the n t h time step, given initial fields (s °, u °, pO) at t = 0. A segregated time stepping scheme incorprating operator splitting and S U P G is described in this section. The following functional spaces will be used: -

e [n'(n)]2},

(8)

"¢/'g ~ {(I~](l) C [ H i ( n ) ] ) 2, (1) -- g on Fo},

(9)

~ g = { ~ [ ~ ~ [H~(fl)]) 2, • = 0 on ro},

(lO)

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121 140

o~ - {qlq~L 2(f2)}

125

(ll)

where Hl(~) is the Sobelev space of functions which are Lebesgue-square integrable oll ~, together with their first derivatives.

3.1. Integration of UCM equations Eq. (7) is solved by an implicit time integration finite element scheme: find s" + ~~ ¢~, V@~ ~,

(

i (]) TS n' ~+- 1U =

S n~~I t d ~ n" VS n+l -- VTuns n+l --Sn+IVun-4-----U

0 ~ - ~ - 2qo~

At

+ f 2qoU"" VOd" dD-~2qo~Pd"u~ n dA,

(12)

where * is the finite element shape function in a Galerkin formulation. Integration by parts has been done to avoid the second derivative of velocity, which has produced a surface integration term concerning the boundary flux of deformation rate. In a consistent SUPG formulation [18], • is to be replaced by @+ W, in which W = (k/u.u)(u.V)O with k depending on local element size and velocity magnitude. When applying SUPG to the UCM in the form of Eq. (7), the second derivative of shape functions can not be avoided, with or without integration by parts. The inconsistent streamline upwinding (SU) [7,19] applied W to the convective term only, which adds artificial diffusivity to the constitutive equation and may effectively modify the problem to some extent. In general consistent SUPG should be applied. For steady flow, however, there is a more economical way of achieving consistency by using SUBS [20], which applys SU on both sides by the equation, i.e. ~ ~u" - Vs n + 1 and ~ ~u" • Vs" are added to left and right hand sides respectively. When steady state convergence is reached the effect of artificial upwinding vanishes in the SUBS formulation. Note in the above system, in an axisymmetric geometry, the three components Sr,., S:: and st: are coupled, and only soo can be solved separately. A full decoupling is readily made by treating the extra viscoelastic shear stress Sr: explicitly in the equations of Srr and s::, and solving for these normal stress components first before solving for st: with updated normal stresses. The asymmetric linear systems resulting from the UCM equatons are solved by a generalized minimal residual (GMRES) iterative solver.

3.2. Zero Reynolds number: the PCG scheme For a creeping flow the second term in Eq. (6) is omitted, and Eq. (6) together with the continuity equation (2) describe a quasi-Stokes (QS) problem, the variational form of which can be described as: Find u ' + J e ~ t ~. p ' + l ~ L 2 ( ~ ) ,

V(I)~7/, Vq~L2(~),

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

126

1 f ( ~ , u "+1 + v/0V~ " Vu,,+ - V~p n+ I) dO

=

;

(~c~,un-V~.s~+°5) df~+

f q V • u ~+ l df~

f(

)

rlo-~-~n-Pe+n's @dA

O,

(13) (14)

where s"+°5 = (s ~ + s n÷ 1)/2 and ~, = 1~At. Note the surface integration is not exactly on the usual total stress traction because we have used the incompressibility to drop the unsymmetric viscous stress term, which is necessary to maintain the symmetric property of the equation for the application of the preconditioned conjugate gradient (PCG) method. Details on the theory of the P C G scheme for QS problems have been given by Glowinski [12]. Here we only give a brief description of the P C G algorithm. We will use k for P C G iteration numbers, not to be confused with the time step n u m b e r n. The following P C G iteration is carried out for every time step: (1) Initially k = 0 and w ~ = p - ~ eL2(f~) is given. The right-hand side in Eq. (13) is denoted as J~, which is calculated from the k n o w n velocity and extra viscoelastic stress of the last time step. (2) Solve for ~ k - ~e {~/'g, ~(/o}, V~ ~ ~/'0,

where the notation {A, B} means A for the initial step, B subsequently. (3) If k > 0, c o m p u t e

p~- 1

f

v • u k - qv ~ - 1 df~

~ V • Z k - l w k - I dr2 qd

f~

If k = 0 , p o = 0 . (4) S e t p k = p k - i - - p w k-l, a n d u k={;C k l , u k I _ p k _ l ~ k - l } . (5) Solve for Qk ~ L2(f~), Vq ~ L2(f~),

f Qkq df~= f V " ukq dfL f~

f~

(6) Solve for yk the Poisson equation:

_ A yk = ~, Qk in fL with the b o u n d a r y conditions ~ Yk/On = 0 on Fo; yk= 0 on F~. (7) Set g~ = qoQ k + yk (8) If k > 0, test convergence:

X.-L. Luo ,' J. Non-Newtonian Fhdd Mech. 63 ~1996) 121- 140

127

If j'n V - ukg k df~ < 8, take p = pk, and u = u k, and stop PCG iteration. (9) Compute ;'k f V " ukg a df~ ~'

=

0

f2

f v • u k Igk

i df~

(10) Set w k = g k + ~,kWk(11) DO k : = k + 1, go to (2) and repeat the process. Although the PCG scheme appears to be much more sophisticated than a conventional velocity pressure decoupling method such as the penalty method, it only involves solving a sequence of scalar Dirichlet problems associated with the Helmholtz operator ~ , 1 - q 0 A and Neuman problems associated with the Laplacian operator A, plus a few scalar integrations. The matrix associated with operator c~,I-rh,A is always well conditioned (since :~, = l/At), thus iterative solvers such as CG, SOR etc. will converge fast. For constant viscosity r/o and fixed time step At, all the matrices in the PCG process remain the same for all time steps, thus the assembly and the LU factorization for a direct solver or the incomplete factorizati:on for an iterative solver need only to be done once for all time steps, which saves a lot of computing time. In addition to fast convergence and high computing efficiency, the PCG scheme directly checks on the divergence of velocity field as a convergence criterion, so the incompressibility constraint is always well satisfied and the pressure solution obtained is usually better than with some other conventional methods such as the penalty type.

3.3. Non-zero Reynolds number: fractional time stepping When the nonlinear convection term can not be omitted the full operator splitting algorithm [12 -14] is adopted. Eqs (6) and (2) for {u "+ I , p " + ' l are solved in three fractional time steps, the variational form of which can be described as follows. First fractional step: Find u ''~ ° ~ 1 " p " + ° E L2(~), W e Vqe L2(~),

@ =

OAt

~ - ~ q o V @ ' V u " + ° - V ~ p "+° dQ

[-[3qoV~'Vu"-@(u~.V)u"-V~.s"+°5]dfl+

J qVu"+°=0. Second fractional step: Find u ''+~ o ~ . Vqbe¢,

qoc~-pn+ns

~dA

(15)

(16)

128

X.-L. Luo /J. Non-Newtonian FluM Mech. 63 (1996) 121-140

q~

(1-20)At

4-flvloV~'Vu'+~-°+@(u~+°'V)u ~+~-° d~

= (-~vloV~'Vun+°+V~p'+°-V~'s'+°S)df~+ Third fractional step: Find un + l e U , p n + l ~ L 2(~'~), ~(I) @ ~]/',

~P

OAt

+

O~qOV(I) .

+ u" ld~ = O,

(17)

Vq ~ L 2(f~),

) ;

v//n+ 1 -- V@p,,+l df~ =

-~(un+l-°.v)un+l-°-V~.sn+°'5]d~+

fqV

Vlo~-~n-pn+n's ~dA.

[-fir/oVa"

qO-~n-Pn+

n

run+

• s]~ dA

/

1

0

(18) (19)

where ~, fl and 0 are the operator splitting parameters satisfying ~ + fl = 1 and 0 < 0 < 1. The sub-problems at the first and third fractional steps are identical and are of the type of steady quasi-Stokes problems, which can be solved by the PCG scheme described in the last section. The sub-problem at the second fractional step is linearized in the convection term to become the classical linear diffusion-convection problem and can easily be solved directly. The values of e and fl are chosen to yield identical Helmholtz operators in all three fractional steps, and this contributes considerably to the overall efficiency of the operator splitting algorithm. A good choice for 0 is 1 - 1 / x / ~ for fast convergence [12-14]. Overall the present OS/SUPG formulation has only first-order accuracy in time, which is sufficient for the present work since we are mainly concerned with the final steady state solution for the time being. For a true transient simulation, the present algorithm can easily be modified to have higher-order accuracy in time at the expense of computing time. For example, a R u n g e - K u t t a type or other predictor-corrector schemes can be applied here, as several transient algorithms for viscoelastic flow have done [25-28]. The SUBS is not valid for true transient problems because it brings in the SU effect during the time stepping process, affecting time dependent accuracy of the solution.

3.4. Finite element implementation The above OS/SUPG algorithm was implemented by the author using a special high level language called Fasttalk, which is a unique feature of Fastflo, a general purpose finite element partial differential equation solving package developed in the Division of Mathematics and Statistics, CSIRO, Australia. The Fasttalk language is a powerful tool for expressing partial differential equation solving algorithms in the finite element context• For example, although the present algorithm is quite complicated the entire Fasttalk code for this OS/SUPG formulation with non-zero Re, transient and three-dimensional capacities contains only about 1000 lines, while a Fortran version with similar capacities could be many times longer. Interested readers can consult Ref. [29], where a World-Wide-Web site is also given for more information•

X.-L. Luo /J. Non-Newtonian Fluid Mech. 63 (1996) 121 140

129

4. Some numerical details For the sphere-in-tube problem, the global quantity of interest is the non-dimensional drag coefficient F K- - 6~r/0 Ua

(20)

as a function of the Deborah n u m b e r 2U De -- - - , a

(21)

where F is the total drag force on the sphere, a is the radius of the sphere, U is the relative velocity of the sphere with respect to the cylinder, 2 is the relaxation time of the U C M fluid and ~1,~ is the viscosity. In the standard benchmark problem the ratio of sphere radius to cylinder radius is a / R = 0.5. 4. I. Mesh refinement and solution strategy

In addition to comparing with well established benchmark results, mesh refinement provides another check on the accuracy of finite element results for complex viscoelastic flows. Some of the previous calculations [2,4,5,21] which failed to converge for low to medium De numbers have used very coarse meshes in terms of current computing power. Only fine meshes have been used in the present study. We have performed systematic calculations with three unstructured triangle meshes: mesh M1 with 23 974 D O F for pressure, velocity and stress variables, mesh M2 with 46 484 DOF, and mesh M3 with 110 004 DOF. Fig. i shows the unstructured triangle mesh M2. The sphere radius is one, cylinder radius is two. The upstream tube length is 10 and the downstream length is 20. There is a high concentration of nodes around the sphere surface and on the centreline near stagnation points.

Fig. 1. Unstructured triangle mesh M2. The partial view is from : = - 3 to z = + 3.

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 1 2 1 - 1 4 0

130

The three meshes used here are not ideal in that the node distribution on the sphere surface is not uniform nor structured. Compared with the triangle meshes used by Jin et al. [9], the meshes used in the present study have a more gradual decrease in node density away from the sphere surface, especially in the stress wake, while the meshes used in Ref. [9] have a rather abrupt decrease in node density away from the sphere surface. On the other hand, the meshes in Ref. [9] have a much better node distribution on the sphere surface. For any Deborah number we started time stepping from uniform velocity u ° = (0, 1), zero pressure pO= 0 and zero extra viscoelastic stress s ° = 0, i.e. we did not use the solution at a lower De number as the initial guess, as is commonly done with steady type schemes. The above strategy proved not only more stable and more robust, but also more economic, because it allowed us to use a very fine mesh to explore the limiting De number question without having to go through small and medium De numbers first.

4.2. Convergence criteria There are three levels of convergence to be monitored during time stepping. At the first level is the convergence of iterative solvers for linear equations--conjugate gradient (CG) method for symmetric systems of Helmholtz and Laplace equations, and generalized minimal residual (GMRES) method for asymmetric systems of the U C M equations. The convergence criterion in the present study for both CG and G M R E S iterative solvers was to let the preconditioned residual be less than 1 0 - 9 . The next level is the convergence of P C G iterations for the QS problem. The criterion was to let the norm of the divergence projection (see Section 3.2).

f v ukg k df~ < 10 - 5 p

(22)

f~

Finally, convergence to steady state was considered to be reached if the following two conditions were satisfied simultaneously,

® -lun+~-unl <0.5 × 10 -3 I

(23)

+

and K n+l _ K n Kn+l

<10

5,

(24)

where K n + 1 and K" are drag coefficients, and un + i and u" are the velocity components at two consecutive time steps. The norm lu] is defined as

(25) in which ui denotes the ith node point value, and the summation is done over all the node points and for both components of the velocity.

X.-L. Luo / J . Non-Newtonian Fluid Mech. 63 (1996) 121-140

131

The above convergence criterion for the steady state checks relative changes in the whole field, not just at isolated points. The dimensionless time step length was typically set to 0.1 in our calculations, which means condition (23) requires the time derivative of the velocity norm to be less than 0.005. We found the extra condition (24) necessary, because often condition (23) was first satisfied without condition (24) being true. Many more time steps were required to satisfy the second criterion on drag force while the first criterion on velocity remained satisfied, with ® decreasing monotonically. On the other hand condition (24) alone is not sufficient either, since it could be fulfilled at some isolated time steps before condition (23) was satisfied. Alternatively, a criterion on the relative change in the norm of all the stres components can be imposed instead of condition (24). We believe condition (24) is also reliable because it effectively checks on the relative change of stresses on the sphere surface, the most sensitive place. It was found that condition (24) required the relative change in the norm of the extra viscoelastic stresses to be of the order of 10 3.

5. Results for U C M fluid

All computations were done on a DEC 300 Model 600 workstation. For the finest mesh M3, convergence for a high De number typically took less than 300 time steps and 7 CPU hours, starting from uniform velocity and zero extra viscoelastic stress field. For example:, to get a converged solution at De = 2.0, M3 with 110 004 D O F took 225 time steps and 6 h 24 min CPU time, M2 with 46 484 D O F took 268 time steps and 2 h 22 min CPU time, while M I with 23 974 D O F took 172 time steps and 52 min CPU time. Based on the CPU time and D O F of M I and M3, the current decoupled method has delivered a speed of CPU time in proportion to N ~3, where N is the total number of unknowns. The convergence of PCG iterations for the QS problem was indeed very fast, averaging about 3 iterations per time step, even for high De numbers. 5. 1. Lhniting Deborah number and the drag force

With mesh M1 steady state convergence was obtained for Deborah numbers up to De = 2.6, but with the finer mesh M2 the limiting Deborah number was De = 2.2. The finest mesh M3 has reached a De value of 2.8, only slightly higher than M l . With the finest mesh M3 we only did calculations for De > 1.6 and for the Newtonian case. Table 1 shows the predicted drag coefficient K for all the three meshes, in comparison with the EEME results of Lunsmann et al. [22]. The EEME results are those with the finest mesh of 28 558 D O F in Ref. [22]. The predictions of drag coefficient K by the three meshes agree very well with each other and with the EEME predictions for De < 1.4, and the maximum relative difference in K among different meshes is about 0.3%. The difference increased with De values. For De >_ 1.6 the maximum relative difference in K is about 2.4%, occurring between M1 and M3 at De = 2.6, the maximum De value reached by M1. The main reason for the discrepancy in K for different meshes is perhaps the huge difference in D O F between the meshes: M3 is about 4.5 times as fine as M1, or the difference in D O F between M3 and M1 is more than 76 000.

132

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

Table 1 OS/SUPG predicted drag coefficient K for the three meshes in comparison with the EEME results of Ref. [22] De

Mesh M1

Mesh M2

Mesh M3

EEME [22]

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8

5.94446 5.65667 5.18879 4.80047 4.51994 4.32807 4.19072 4.11039 3.99686 3.91849 3.82147 3.74469 3.67342 3.59848

5.94755 5.65774 5.18949 4.80869 4.52494 4.33268 4.19783 4.10757 3.97283 3.88511 3.78158 3.72772

5.94726

5.94739 5.65981 5.18640 4.79968 4.52391 4.33414 4.20572 4.12007 4.06348

3.94653 3.86386 3.74029 3.66375 3.59247 3.51423 3.43265

It is perhaps not very surprising that a larger difference in drag force predictions occurs at larger D e numbers among several fixed meshes, since higher stress values and gradients occur at higher D e numbers. The relative numerical accuracy m a y reduce as the D e value increases, especially for the coarser mesh M1. As seen from Table 1, at high D e numbers the finer mesh always predicted slightly lower values o f K (or a stronger viscoelastic effect) than the coarser mesh. It is believed that the predictions of the finer mesh M3 are more accurate than those of M1 and M2. Although the agreement a m o n g different meshes is already quite good, we believe that better structured meshes can deliver even better mesh independence in the drag force prediction. In our unstructured meshes the positioning of nodes on the sphere surface is quite arbitrary or completely unstructured, and the intervals of arc length between nodes are not uniform. The kind o f well structured node distribution on the sphere surface used by Jin et al. [9] is much preferred. Fig. 2 shows a comparison of drag coefficient K a m o n g seven different methods, which include those of L u n s m a n n et al. [22], Jin et al. [9], Crochet and Legat [23], Fan and Crochet [24]. Rasmussen and Hassager [30] and K h o m a m i [I1], in addition to the present results. While results o f all three meshes o f the present work are shown in Fig. 2, only results o f the finest mesh are shown for the other methods. Again, the agreement among different methods and meshes is very good up to about D e - - 1.4. At higher D e numbers (De >_ 1.6) there is a major difference between this work and the other computations based on mixed methods. The present decoupled calculation converged at higher D e values and predicted no levelling off in the drag coefficient K, while other coupled calculations encountered an earlier loss o f convergence and predicted an obvious trend of levelling off in K just before divergence.

X.-L. Luo / J . Non-Newtonian Fluid Mech. 63 (1996) 121 140

133

5.2. Velocity and stress fields"

The solution fields for a Newtonian fluid and for the UCM fluid at De = 1.0, 2.0 and 2.8 (the highest De value reached) are shown in Figs 3 and 4-6, respectively. All the contour plots have 16 contour levels equally spaced between the minimum and the maximum values. Since we have solved the U C M equation in terms of the extra viscoelastic stress tensor s, for the purpose of post-processing the total extra stress tensor • needs to be recovered by adding the viscous stress tensor to s. Explicitly computing the viscous stress tensor is subject to extra numerical inaccuracy because of the numerical differentiation of the velocity field. Of course this particular inaccuracy only appears in the post-processing, not in the actual simulation. One expects the minimum and maximum values of T shown in the figures to contain the post-processing inaccuracy just explained. The changes in the velocity and stress fields caused by viscoelasticity are apparent. The UCM velocity field has a downstream shift, which can be seen from both the Ur and the u_ contours. We have added an extra contour level very close to the centerline in the streamline plot in Fig. 6A for De = 2.8, which shows that there is no reverse flow or recirculation in the wake even at the highest De number. Consistent with the downstream shift of ur and u_ contours, the streamlines also shifted downstream. This p h e n o m e n o n can be qualitatively explained as a memory effect of the UCM fluid: take the cylinder as stationary, the sphere falls from top to bottom in a viscoel:astic fluid which has certain memory. Fluid in the sphere wake still 'remembers' the presence of the sphere a short while ago, thus the streamlines shift downstream as if the fluid still 'feels' the presence of the sphere downstream. The higher the De number, the stronger this memory effect. In a Newtonian flow around a sphere in a certain Reynolds number range, there is also a downstream shift of the main streamlines which is due to a recirculation behind the sphere. Fig. 6A shows there is no recirculation at all behind the sphere, and the downstream shift of streamlines observed here is certainly a unique feature of viscoelastic flow.

tO-

tl').

tO.

010

015

1'.0

115 De

210

215

310

Fig. 2. Comparison of drag coefficient K from seven different methods. - - - , this work, M 1; A~ this work, M2; , , this work, M3; V, Lunsmann et al. [22]; + , Jin et al. [9]: ?~, Crochet and Legat [23]: x , Khomami ( 1 1 ) : . , Rasmussen and Hassager [30]; Z~, Fan and Crochet [24].

134

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121 140

I mln.

= -2.00E00

rain.

= -3.76E

max.

= 0.00EO0

max.

= +3.76E

(a)

-

1

rain.

=

O.OOEO0

rain.

= -3.15E00

1

max.

= +1.66E00

max.

= +l.81E

rain + 1

max.

= =

503E00

rain

* 5 03EO0

max

299E00 : *8 01EO0

rain max

=

5.03E00

rain.

+5 03E00

max

=

1.02EO0

= +l,02E00

(b)

Fig. 3. A. Contour plots of streamline 0, velocity components ur, u: and pressure p for the Newtonian solution. Sixteen contour levels are equally spaced between minimum and maximum values. B. Contour plots of stress components v,, rr:, r:: and r0o for the Newtonian solution. The relative contour levels are the same as Fig. 3A. A n o t h e r interesting character of the U C M flow is revealed by the u_ contours. The u: contours become m o r e and more 'elongated' downstream as the De number increases, which means that at a higher D e number it takes a longer distance to accelerate from zero speed at the stagnation point to unity downstream in the centreline, and similarly it takes a longer distance to decelerate from the m a x i m u m speed to unity for fluid in the middle between the sphere and the cylinder wall. In other words the magnitude of the elongational rate immediately downstream of the sphere becomes smaller at a higher De number. This behaviour is not very surprising--it is essentially the resistance of the U C M fluid to stretching. The higher the De number is, the more rigid the U C M fluid becomes. Nevertheless, the elongational stress is still an increasing function of the D e number, despite the fact that the elongational rate decreases with De. The non-Newtonian stress fields are characterized by large gradients and thin boundary layers adjacent to the sphere surface. At high D e numbers only the axial stress component r:: has considerable magnitude away from the sphere surface, all the other components are only concentrated on the sphere surface, bearing in mind we have 16 contour levels equally distributed between the m i n i m u m and the maximum. Although the stress contours for the

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 ~1996) 121 140

135

highest De value 2.8 shown in Fig. 6B do not appear to be as smooth as those at lower De values, this did not seem to affect the numerical convergence for the computation at De = 2.8. Using 10 levels equally spaced between the minimum and zero and between zero and the maximum, Lunsmann et al. [22] showed the stress field for the Oldroyd-B model with [Y = 0.5 at De = 1.0. In order to compare with their contour plot in detail, Fig. 7 has used the same contour levels as in Ref. [22] to show our results for the same case of the Oldroyd-B model with [] = 0.5 at D e - - 1.0. All the details shown in Fig. 7 appear to be the same as in the corresponding plot in Ref. {22]. We note an opposite sign convention was used by Lunsmann et al. [22] for the stress components. 5.3. Loss ~[~ convergence The loss of convergence of numerical computations for the sphere-in-tube problem at high De numbers has been attributed to both the presence of the birefringent strand and the high stress gradients close to the sphere [31]. Some attention has also been paid to the elongational

t inilb ma,,

(a)

2[IOE00 : .O00EO0

rain = max=

340E 446E

1 1

rain = O.00EO0 max.=+l.70E00

nain max

= 7.00E00 ~157E+1

rain max

= 606E ~25flE~

1 1

nlln = 2 02E • I n~ax~ +312E~ 1

j ~

rain max

73~1"~ ~728E

I ~ I

mill nl~x

406E I +393E(1{}

(b)

Fig. 4. A. C o n t o u r plots o f streamline ~, velocity c o m p o n e n t s u~, u: and pressure p for the U C M solution at

De = 1.0. The relative c o n t o u r levels are the same as Fig. 3A. B. C o n t o u r plots of stress c o m p o n e n t s r,,, r,_. r_. and r~,, for the U C M solution at De = 0.1. The relative c o n t o u r levels are the same as Fig. 3A.

X.-L. Luo /J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

136

3

I'

rain. max.

-

2 00E00 -

0.00EOO

rain

-

3.25E

1

rain.

max

-

+4.79E

1

max..-

-

O 00E00

rain

+ 1.72E00

max.

-I.03E +1.68E

+ 1 + 1

(a)

rain

5 721'~

1

rain

2 6(IE

4 1

nlin

ma×

.3 29E

~ 1

max

~ 3 57~

~ 1

i]l,'ix.

:

5 52E +971E

I

rain

+ 1

max.

-

302E +9.06E00

(b)

Fig. 5. A. C o n t o u r plots of streamline ~,, velocity components ur, u_ and pressure p for the U C M solution at De = 2.0. The relative contour levels are the same as Fig. 3A. B. C o n t o u r plots of stress components rrr, rrz, r.: and ~oo for the U C M solution at De = 2.0. The relative contour levels are the same as Fig. 3A.

behaviour on the centreline [6,21,32]. It was suggested by Zheng et al. [6] that the aphysical elongational behaviour of a U C M fluid may be responsible for the loss of convergence to steady state. The present work found no evidence to support the above conjecture. The elongational stress and its gradient on the centreline near the stagnation points were found to be not as large as they were on the sphere surface at all the numbers covered so far by our computations. The sign of the elongational stress on the centreline is also physically correct: before the front stagnation point the flow upstream is purely compressed < 0) and r - z - rrr is negative, while after the rear stagnation point the flow downstream is purely elongated (du~/dz > 0) and r - - - rrr remains positive. A semi-analytic expression can be written for the elongational stress on the centreline. On the centreline the U C M equation for r____becomes a decoupled ordinary differential equation because ur = 0 and rr------0 on the centreline,

De

(duz/dz

( dr=

~u~r_.:) +

2 u_ d----~-- 2 ~

Ou_

r-: = 2qo c ~ '

(26)

1

X.-L. Luo /J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

137

which is linear in r:: and of first order, and it can be integrated exactly as r~z =

f e ~"d: q dz +

e -iWd:

c

e -.tw d:,

(27)

where w = (1 -2)t(~u:/O:))/(2u:), q = (2tlo/2u:)(Ou:/~z) and c is an integration constant. The value of w (but not q) becomes infinity at the stagnation point, and the value of r:= can be shown to be zero there. The elongational stress on the centreline is also well behaved in the vicinity of the stagnation point, for it can easily be shown from Eq. (26) that r=:(~s) =

2r/o~ 1 -

+ O(6s),

(28)

2~

where 6s is an arbitrarily small distance downstream from the stagnation point, and ~ = ~u:/?z is the varying elongational rate which approaches zero as 6s goes to zero. The above argument is valid for both stagnation points and also applys to the r~ component. Thus the U C M stresses on the centreline are Newtonian-like near the stagnation points.

j - -

) rain

-2.00E00

rain

=

max

= 0 00E00

max

=

(a)

3.00E ~497E

-

1

rain.

= 0,00E00

mln.

=

1

max.

= +1.73E00

max

= + L.99E

IlSE

+ 1 ~- 1

min max

=

- 1 60E00

= +405E

+ 1

rnin

= -3.05E

tnax.

= +368E

+ 1 + 1

rain max

=

2 76E00

rain

=

=

t 1 06[:] L 2

max

:

I 23E00 ~1.27E

~ 1

(b)

Fig. 6. A. C o n t o u r plots of streamline ~b, velocity c o m p o n e n t s u,., u: a n d pressure p for the U C M solution at De = 2.8. The relative c o n t o u r levels are the same as Fig. 3A except one extra level very close to the centreline. B. C o n t o u r plots of stress c o m p o n e n t s r~r, rr:, r:: a n d z0o for the U C M solution at De = 2.8. The relative c o n t o u r levels are the same as Fig. 3A.

138

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121 140

r

rnin,

=

max

=+1.56E+1

I 17EO0

rain max=

=

1.09E 1-2.24E+

# 1 1

rain max

~

1.40E00

=+590E

" I 1

max

=

6 09E

1

~1 8 4 E 0 0

Fig. 7. C o n t o u r plots of stress c o m p o n e n t s rrr, ~r~-, r.- a n d Zoo for the Oldroyd-B solution at fl = 0.5 a n d De = 1.0. T h e relative c o n t o u r levels are the same as those used in Ref. [22].

Upstream integration (27) can start from the inlet and finish at the front stagnation point, while downstream it can start from the rear stagnation point. In both cases the starting value for r._ is zero, which means the constant c in (27) is zero for both upstream and downstream segments. Thus the sign of r~_. is the same as q which is always negative upstream and positive downstream, therefore the aphysical behaviour of unbounded or negative Trouton viscosity which is possible in a uniform elongational flow just can not occur here on the centreline, unless q changes its sign on the upstream or downstream centreline, which has never been observed in any numerical prediction for the sphere-in-tube problem so far. The fact that M3 (DOF = 110 004) only slightly outperformed M1 ( D O F - - 2 3 974) in terms of De limit is interesting, and seems to suggest that after D O F exceeds a certain size (20 000, say), a dramatic increase in D O F does not necessarily improve performance in proportion, and the node distribution near the sphere may become a decisive factor. The reason that M1 performed better than M2 ( D O F - - 4 6 484) is perhaps due to mesh M1 having a better overall node distribution than M2. The present triangle mesh generator does not give the user much control on how the nodes are distributed along a specified curve. It may well be possible to get convergent steady state solutions at even higher De numbers by using fine meshes with carefully designed node distribution, preferably adaptive to the stress field. On the other hand, there is

X.-L. Luo /J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

139

always the possibility of the onset of physical instability inherent in the UCM fluid at high De numbers, or even aphysical stress fields occurring in the vicinity of the sphere.

6. Conclusions We have demonstrated the capacity of our decoupled operator splitting algorithm for solving highly viscoelastic flows. The high efficiency of the present decoupled OS/SUPG formulation has made it possible for us to use very fine meshes for complex viscoelastic flow computations on a modern workstation computer, obtaining a converged solution in a few hours with a DOF over 105. This computer efficiency (CPU speed and memory) is hard to achieve with a coupled method. The robustness and accuracy of the present OS/SUPG code have also enabled us to reach steady state convergence for De numbers up to De = 2.8 for the UCM fluid flow around a sphere in a tube, and the results agree very well with those of the more expensive coupled methods for De _< 1.4. For De >_ 1.6 the present decoupled scheme converged at De numbers higher than in other mixed methods and predicted no obvious levelling off in the drag coefficient K, while other coupled calculations encountered an earlier loss of convergence alter predicting an obvious trend of levelling off in K. This qualitative discrepancy needs further investigation. The following questions have now become even more pressing: should the drag coefficient K start to increase at a certain De value of UCM fluid?; does K approach an asymptotic value if it always decreases monotonically?; what is the magnitude of the asymptotic value for K if there is one? The present study suggests that the loss of convergence to steady state was m~t due to elongational behaviour on the centreline becoming aphysical. The loss of convergence due to high stress gradients around the sphere can be further investigated numerically by using an even larger number of DOF, which the present decoupled code can still afford. In future work the mesh can be even finer, but more effectively it needs to be better designed to have an optimal node distribution around the sphere and in the stress wake, preferably adaptive to the stress gradients. It is also highly desirable to design numerical experiments to study the existence and the nature of the birefringent tail [31]. The sphere problem solved here is two dimensional and steady; the full capacity of the present decoupled OS/SUPG algorithm will be better demonstrated in simulations of transient and three-dimensional non-Newtonian flow problems. The present decoupled formulation is also applicable to integral models, and it can be expected to perform well in simulations of non-Newtonian flow with high Reynolds numbers.

Acknowledgements This work is part of the Fastflo project which has been carried out in the CSIRO Division of Mathematics and Statistics, supported by the Department of Science, Industry and Technology of the Australian Government, BHP and Compumod. The author wishes to thank his colleague Dr. Nick Stokes and the three referees for their helpful comments on a draft of this paper.

140

X.-L. Luo / J. Non-Newtonian Fluid Mech. 63 (1996) 121-140

References [1] M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, 1984, [2] O. Hassager and C. Bisgaard, J. Non-Newtonian Fluid Mech., 12 (1983) 153. [3] B. Caswell and M.Viriyayuthakorn, J. Non-Newtonian Fluid Mech., 12 (1983) 13. [4] X.-L. Luo and R.I. Tanner, J. Non-Newtonian Fluid Mech., 21 (1986) 179. [5] F. Sugeng and R,I. Tanner, J. Non-Newtonian Fluid Mech., 20 (1986) 281-291. [6] R. Zheng, N. Phan-Thien and R.I. Tanner, J. Non-Newtonian Fluid Mech., 36 (1990) 27. [7] J.M. Marchal and M.J. Crochet, J. Non-Newtonian Fluid Mech., 26 (1987) 77. [8] R.C. King, M.R. Apelian, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 29 (1988) 147. [9] H. Jin, N. Phan-Thien and R.I. Tanner, Comput. Mech., 8 (1991) 409-422. [10] D. Rajagopalan, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 36 (1990) 159. [11] B. Khomami, paper presented at the Eighth International Workshop on Numerical Methods in Non-Newtonian Flows, Cape Cod, October 1993. [12] R. Glowinski, Finite element methods for the numerical simulation of incompressible viscous flow, Lectures in Applied Mathematics, Vol. 28, AMS, Providence, RI, 1991, pp. 219-301 [13] M.O. Bristeau, R. Glowinski and J. Periaux, Comput. Phys. Rep., 6 (1987) 73-87. [14] E. Dean, R. Glowinski and C.H. Li, Comp. Phys. Comm., 53 (1989) 401-39. [15] X.-L. Luo, Operator splitting computation of turbulent flow in an axisymmetric 180° narrowing bend using several k-~ models and wall functions, Int. J. Num. Methods Fluids, (1996), in press. [16] C.G. Speziale, On nonlinear k - - I and k - - ~ models of turbulence, J. Fluid Mech., 178 (1986) 459-475. [17] R.S. Rivlin, The relation between the flow of non-Newtonian fluids and turbulent Newtonian fluids, Q. Appl. Math., 15 (1957) 212. [18] A.N. Brooks and T.J.R. Hughes, Comput. Methods Appl. Mech. Eng., 32 (1982) 199. [19] X.-L. Luo and R.I. Tanner, J. Non-Newtonian Fluid Mech., 31 (1989) 143-162. [20] R.I. Tanner and H. Jin, J. Non-Newtonian Fluid Mech., 41 (1991) 171. [21] 5th Workshop on Numerical Methods in Non-Newtonian Flows, J. Non-Newtonian Fluid Mech., 29 (1988). [22] W.J. Lunsmann, L. Genieser, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 48 (1993) 63. [23] M.J. Crochet and V. Legat, J. Non-Newtonian Fluid Mech., 42 (1992) 283. [24] Y. Fan and M.J. Crochet, J. Non-Newtonian Fluid Mech., 57 (1995) 283. [25] P.J. Northey, R.C. Armstrong and R.A. Brown, J. Non-Newtonian Fluid Mech., 36 (1990) 109. [26] P. Townsend and M.F. Webster, Proc. NUMETA 87, 2 (1987) T12/1-11. [27] F. Olsson, J. Non-Newtonian Fluid Mech., 51 (1994) 309-340. [28] T. Sato and S.M. Richardson, J. Non-Newtonian Fluid Mech., 51 (1994) 249-275. [29] N. Stokes, Fasttalk Manual for Fastflo Version 2.4, CSIRO Division of Mathematics and Statistics, April 1995. More information is available from WWW site: http://www.mel.dms.CSIRO.AU/nick/world/CC05.html. [30] H.K. Rasmussen and O. Hassager, J. Non-Newtonian Fluid Mech., 46 (1993) 289. [31] O.G. Harlen, J. Non-Newtonian Fluid Mech., 37 (1990) 157. [32] M.D. Chilott and J.N. Rallison, J. Non-Newtonian Fluid Mech., 29 (1988) 381.