- Email: [email protected]
An incremental difference formulation for viscoelastic flows and high resolution FEM solutions at high Weissenberg numbers
J. Non-Newtonian Fluid Mech., 79 (1998) 57 – 75
An incremental difference formulation for viscoelastic flows and high resolution FEM solutions at high Weissenberg numbers X.L. Luo * CSIRO Mathematical and Information Sciences, Locked Bag 17, North Ryde, NSW 2113, Australia Received 26 August 1997; received in revised form 14 January 1998
Abstract An incremental difference formulation (IDF) has been developed to further improve the efficiency and stability of the operator splitting algorithm for viscoelastic flows, taking advantage of the pseudo linearity of the governing equations. Using this new IDF/OS scheme and up to 200000 degrees of freedom, high resolution mesh-convergent solutions for the falling sphere problem are obtained up to Wi = 2.8, with a minimum value of K = 4.02 predicted at Wi=2.2. Beyond Wi =2.8, the solutions become significantly mesh-dependent, although steady convergence with individual meshes are still possible. Well structured quadrilateral elements lead to better accuracy than the entirely unstructured triangle meshes previously used. There is some local mesh-dependence of Tzz on the centreline starting at a finite distance downstream from the rear stagnation point, even for globally mesh-convergent solutions. This local mesh-dependence is most likely due to the presence of a thin boundary layer forming along the centreline in the stress wake, where Tzz exhibits a sharp radial gradient and its local maximum is not on the centreline but slightly off the centreline inside the boundary layer. Further mesh refinement in the radial direction, perhaps reducing the size by an order of magnitude, is needed to study the nature of this boundary layer. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Difference formulation; Operator splitting; Viscoelastic drag; Maxwell fluids; Sphere
1. Introduction In the last few years, several calculations of the sphere-in-tube problem for UCM fluid have reached Weissenberg numbers around or even beyond 2, with some disagreement among different groups on the predicted drag coefficients, or even among different meshes of the same group [1 – 9]. With very few exceptions, mesh-convergence for the sphere problem at Weis* Fax: + 61 0293 253200. 0377-0257/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S0377-0257(98)00063-9
58
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
senberg numbers around and beyond 2.0 has not been as convincingly demonstrated as lower Weissenberg numbers. The detailed and more meaningful comparison of stress distribution curves and contours among different calculations is scarce, owing to lack of published data. One of the purposes of this paper is to provide detailed data for future comparison of solutions at Weissenberg numbers around and beyond the value of 2.0. Jin et al [1] extended the explicitly elliptic momentum equation (EEME) formulation [10] to a class of constitutive equations of the Maxwell type, and obtained results for the sphere problem at high Weissenberg numbers with triangular meshes which are well-structured around the sphere surface but unstructured away from the sphere. Lunsmann et al. [2] were perhaps the first to convincingly demonstrate mesh-convergence of solutions up to Wi=1.6. Fan and Crochet [4] used a high-order finite element method to the sphere problem. The p-convergence up to Weissenberg number 2.0 is verified by means of a residual norm as well as the stress profiles in the stress boundary layer. They estimated an asymptotic value of 4.02 for the drag coefficient at Wi\ 2.0. Sun et al. [6] developed an adaptive viscoelastic stress splitting (AVSS) method and applied it to the sphere problem to obtain solutions up to Weissenberg number of 2.8. More recently, Baaijens et al [7] utilized a stabilised discontinuous Galerkin method, first introduced by Fortin and Fortin [11], for viscoelastic flows. They were able to reach Weissenberg number 2.7 with a moderate number of degrees of freedom, and predicted a minimum drag coefficient of 4.048 at Wi= 2.0, followed by a slow increase in the drag. Huang et al [8] used a new unstructured control volume method to get solutions at high Weissenberg numbers for the sphere problem. Unsmooth stress fields were reported at high Wi values [8]. Yang and Khlomami [9] employed an hp-adaptive finite element technique to the sphere problem to get mesh-convergent results beyond Weissenberg number 2.0. Their results showed that a finer mesh may give a much high Tzz peak than a coarser mesh in the stress wake at high Weissenberg numbers, indicating a possibly unbounded stress wake. The author developed an operator splitting algorithm with a decoupled formulation for transient as well as steady viscoelastic flows [5]. The high efficiency and rebustness of this algorithm enabled high resolution calculations of the sphere problem with over 105 degrees of freedom on a moderate workstation. Unfortunately, our previous calculations used completely unstructured triangle meshes (even on the sphere surface) which led to considerable inaccuracy and poor mesh-convergence at high Weissenberg numbers; This deficiency was due to the inability of unstructured triangle meshes to resolve the thin elastic stress boundary layers on the sphere properly. This inability of unstructured mesh with an uneven element height to resolve thin boundary layers has also been confirmed in the author’s work in turbulence modeling [12]. We point out that the unstructured meshes used by Jin et al [1], Sun et al. [6] and Huang et al [8] are well structured around the sphere surface. The main objective of the present work is to introduce an incremental difference formulation (IDF) for UCM fluids, taking advantage of the pseudo linearity of the momentum equation and UCM equations in a decoupled numerical procedure, and to present comprehensive high resolution results at high Weissenberg numbers for the sphere problem for future comparison with other work. In particular, apart from drag coefficient predictions, detailed stress distribution curves along the centreline and on the sphere surface and contour plots of velocity, pressure and stress fields will be provided, showing all the minimum and maximum values of the solutions.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
59
Much effort in this work involves getting mesh-convergent solutions, rather than trying to reach highest Wi number with a single mesh. Well structured quadrilateral elements with close to 200000 degrees of freedom have been used for high resolution calculations. We found it essential to check solution residuals of the constitutive equations in a decoupled scheme. Our present numerical work predicted a minimum drag coefficient value of 4.02 at Wi =2.2, after which the drag starts to increase slowly with Wi. Mesh-convergent solutions were obtained up to Wi =2.8, beyond which the solution became significantly mesh-dependent, although steady state results with individual meshes were still possible. New findings include the apparent local mesh-dependence of Tzz values on the centreline, which starts at a finite distance downstream from the rear stagnation point. A thin Tzz stress boundary layer forms along the centreline whose position coincides with the local mesh-dependence region. Due to its small thickness this thin boundary layer is likely to be missed by numerical solutions with relatively coarse meshes.
2. Governing equations and the incremental difference formulation Considering the following governing equations for the creeping motion of an incompressible viscoelastic fluid: du + 9p = 9 · T dt
(1)
9 · u = 0,
(2)
where u is velocity, p is pressure and T is the extra-stress tensor. For an upper-convected Maxwell fluid the extra-stress T satisfies the following constitutive equation: l
DT + T = 2h0d, Dt
(3)
where l is the relaxation time, h0 is the constant viscosity, d is the rate of deformation tensor, DT is the upper-convected time derivative: and Dt DT = u · 9T − LT − TLT, Dt
(4)
where the temporal term is omitted for steady flow and L is the velocity gradient tensor. Taking the stress splitting approach, the above system can be reformulated in terms of extra stress S and the Newtonian viscous stress N =2h0d: S T − N.
(5)
The momentum equation and the constitutive equation now become (u − h0Du + 9p =9 · S, (t
(6)
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
60
l
DS DN +S = − l . Dt Dt
(7)
2.1. Incremental difference formulation Eqs. (2), (6) and (7) form a non-linear system of a mixed type, the solution of which in general becomes harder as the Weissenberg number becomes larger. In the past, it is common for a decoupled scheme to start the solution for a higher Weissenberg number from a previously obtained solution at a lower Weissenberg number. In order words the solution at a lower Weissenberg number is often used as the initial guess for the Picard iteration procedure at a higher Wi. This so called incremental approach is usually not suitable for transient problems where one has to start the calculation with a true initial state. In our previous work [5], a transient operator splitting algorithm was developed which does not make use of the solution at a lower Weissenberg number—the time stepping for a given Weissenberg number always starts from zero initial values For steady flow, it is time consuming and unnecessary to do time stepping for all Weissenberg numbers starting from a true initial state. With our operator splitting code, one can of course utilize the incremental approach to save computer time without affect the final steady state solutions. However, in the present work we wish to go one step further: not only we make use of the solution at a lower Weissenberg number, but also we solve only the difference between two solutions at two Weissenberg numbers, taking advantage of the pseudo linearity of the equations in a decoupled procedure. This should improve both the efficiency and stability of the calculations for high Weissenberg numbers, as we will demonstrate in this paper. Let us denote two relaxation times as l1 and l2, assuming l2 \l1. The corresponding field variables at the two relaxation times will be distinguished by subscripts ‘1’ and ‘2’ respectively. The governing equations at l1 are (u1 (8) − h0Du1 + 9p1 = 9 · S1, (t (9) 9 · u1 = 0, DS1 DN1 +S1 = − l1 , (10) l1 Dt Dt and the equations at l2 are (u2 (11) − h0Du2 + 9p2 = 9 · S2, (t (12) 9 · u2 = 0 DS2 DN2 l2 +S2 = − l2 . (13) Dt Dt The momentum equation, continuity equation and the UCM equation are all pseudo linear in that when velocity, pressure and elastic stress are considered separately as in a segregated scheme, there is no quadratic or any higher order terms in the equations. This pseudo linearity is the basis of our new incremental difference scheme. By subtracting both sides of momentum and continuity equations at l1 from those at l2, one gets
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
61
(du − h0Ddu+ 9dp = 9 · dS, (t
(14)
9 · du =0
(15)
where d is used to denote the differences of variables due to an increment in l: du = u2 − u1
dp= p2 − p1
dS=S2 −S1.
(16)
The UCM equations need more than straight forward subtraction because of the cross terms involving tensor multiplication’s of stress and velocity gradient. Define the upper-convective D1 D2 Dd derivative operators , and (without the temporal term) corresponding to velocity Dt Dt Dt fields u1, u2 and du as D1( ) = u1 · 9( )− L1( )− ( ) −L1T, Dt
(17)
D2( ) = u2 · 9( )− L2( )− ( )L2T, Dt
(18)
Dd ( ) =du · 9( ) − dL( )− ( )dLT Dt
(19)
where ( ) represents an arbitrary stress tensor. It follows that D2( ) D1( ) Dd ( ) = + . Dt Dt Dt
(20)
Splitting S2 into S1 and dS and moving all S1 terms to the RHS, Eq. (13) becomes l2
D2dS D2S1 D2N2 + dS= − l2 − S1 − l2 Dt Dt Dt
(21)
The RHS of Eq. (21) can be expressed entirely in terms of increments of either dS=S2 −S1, or dN = N2 − N1 or dl =l2 − l1. Making use of Eq. (20), the RHS of Eq. (21) can be rewritten as RHS = − l2
n
D1S1 Dd S1 D(N1 +dN) Dd N2 − S1 − l2 −l2 + . Dt Dt Dt Dt
(22)
Now from Eq. (10), we get − l2
D1S1 D1N1 l2 −l2 = S1 Dt Dt l1
(23)
Substituting Eq. (23) into Eq. (22), we get RHS =
dl Dd S1 Dd N2 D1d N S1 − l2 . + + l1 Dt Dt Dt
Denoting S1 + N2 as G and substituting Eq. (24) into Eq. (21), one gets the final form
(24)
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
62
l2
D2dS dl DdG D1dN + dS= S1 − l2 + . Dt l1 Dt Dt
(25)
To sum up, Eqs. (14), (15) and (25) complete our formulation in terms of variable difference due to an increment in l. Note the operator Dd/Dt contains the difference in velocity and its gradient. A few necessary remarks include: Remark 1. All of the equations in the incremental difference formulation are exact, there is no approximation. Remark 2. All of the unknowns are in terms of difference between solutions at two Weissenberg numbers which significantly reduces the maximum values of the unknown stress components and their gradients in both the momentum and the UCM equations. For example, at l1 = 2.0 and l2 = 2.2, numerical results in the present work show that the maximum value in dSzz is only 19.9, while the maximum value of the full elastic stress component Szz at Weissenberg number 2.2 is 145.1, more than seven times larger than the difference. Consequently, the gradient 9dSzz is also much smaller than the full value 9Szz Remark 3. The incremental difference momentum and continuity equations (Eqs. (14) and (15)) are of exactly the same form as the original equations (Eqs. (1) and (2)), and the reformulated UCM equation takes the same form in LHS as the original, but with a few extra terms on the RHS. Therefore, the methods of solution for these equations can remain virtually unchanged from those for the full variable equations. In any decoupled scheme, the operators D1 D2 Dd and can all be evaluated readily from u1 and du and the same is true for N2 and dN, Dt’ Dt Dt given u1 and du. Remark 4. The Dirichlet boundary conditions for velocity difference du are simply zero everywhere in the new formulation. The boundary condition on dS can also be readily worked out, which is simply zero for the sphere problem.
3. IDF/OS scheme for viscoelastic flows As the governing equations in an IDF take the same form as the original full variable ones, IDF can be readily combined with our previously developed operator splitting (OS) algorithm [5] to produce a more efficient and more stable scheme-IDF/OS. Details regarding the OS algorithm can be found in [5]. Briefly, a transient algorithm is adopted to solve the nonlinear system of momentum, continuity and UCM equations in a segregated manner. The pre-conditioned conjugate gradient (PGG) algorithm of Glowinski [13] is used to solve the Quasi Stokes (QS) equation with viscoelastic stress contributions, which only involves solving a sequence of scalar Dirichlet problems associated with the Helmholtz operator and Neuman problems associated with the Laplacian operator, plus a few scalar intergrations. The matrix associated with the Helmholtz operator is always well conditioned (as the time step can be small), thus iterative solvers such as CG, SOR etc. will converge fast. For constant viscosity h0 and fixed time step Dt, all of the matrices in the PCG process remain the same for all time steps. As a result, the assembly and LU factorization for a direct solver or the incomplete factorization for an iterative solver need only be carried out once for all time steps, which saves on computing time considerably.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
63
In addition to fast convergence and high computing efficiency, the PCG scheme directly checks on the divergence of the velocity field as a convergence criterion. As a result, the incompressibility constraint is always well satisfied and the pressure solution obtained is usually more accurate than some other conventional methods such as the penalty type. A consistent SUPG method is used for the UCM equations, along with generalized minimal residual (GMRES) iterative solver for the asymmetric systems. Integration by parts is used to avoid the second derivatives in Eq. (25), and the extra surface integration term concerning the boundary flux of the deformation rate can be readily calculated.
4. High resolution solutions of the falling sphere problem After the IDF/OS scheme was implemented by the author using a high level FEM language Fasttalk [14], the new code was applied to the standard sphere-in-tube bench mark problem. As discussed in Section 1, our previous calculations with the OS scheme [5] used totally unstructured triangle meshes, which led to considerable inaccuracy at high Weissenberg numbers due to the inability of the unstructured mesh to resolve the thin elastic stress boundary layers. In the present work, we only used well-structured quadrilateral elements. A direct comparison of the old OS scheme with the new IDF/OS scheme was carried out with the same quadrilateral mesh, in order to demonstrate the higher efficiency and the stability of the latter. A well-structured triangle mesh proves to be equally as good as the quadrilateral ones. The improved efficiency of our new formulation allowed us to attack this difficult problem with very large degrees of freedom (up to 200000) on a moderate workstation. Our aim is to present detailed results on the stress as well as velocity and pressure fields to facilitate careful comparisons with other works in the future.
4.1. Some numerical details For the sphere-in-tube problem, the global quantity of interest is the non-dimensional drag coefficient: K=
F 6ph0Ua
(26)
as a function of the Weissenberg number Wi=
lU , a
(27)
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, l is the relaxation time of the UCM fluid and h0 is the viscosity. In the standard benchmark problem, the ratio of sphere radius over the cylinder radius is a/R = 0.5. As in previous OS schemes, a quadratic interpretation is used for velocity, the extra stresses, and linear for the pressure.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
64
4.1.1. Solution strategy We initially obtained solution at Wi= 0.6 with our previous conventional OS/SUPG formulation, i.e. the full velocity and stress component formulation. Consequently, we let l1 =0.6 and dl = 0.2 and used our new IDF/OS scheme to advance the solutions to higher Weissenberg numbers. The l1 was of course changed to l1 +dl after each convergent solution, ready for the next increment in Wi. The increment dl was fixed at 0.2 for all calculations. 4.1.2. Con6ergence criteria—checking equation residuals We believe all numerical work should report the convergence criteria used, which unfortunately is not always the case in the literature. In our previous work [5], the following convergence criteria for steady state solutions has been used: Convergence to steady state was considered to be reached if the following two conditions were satisfied simultaneously, u n + 1 − u n B eu u n+1
(28)
K n+1 − K n B eK, K n+1
(29)
and
where K n + 1 and K n are drag coefficients, and un + 1 and un are the velocity components at two consecutive time steps. The norms are the usual L2 norm. We have set eu =0.5×10 − 3 and eK = 10 − 5, respectively in our previous work. With the new formation, we found an even tighter convergence test can be applied to the momentum equation—we let eu = 10 − 4. Furthermore, due to the numerical solution of the UCM
Fig. 1. Partial view of mesh M1 and M3.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
65
Table 1 Mesh details
M1 M2 M3
No. nodes
No. ele.
df
Size (ds/dr)
17205 23653 31061
5580 7700 10140
108 810 149 618 196 506
0.0175/0.00202 0.0142/0.00166 0.0113/0.00128
The size ds is the arc length of elements on the sphere surface and dr is the radial length.
equation being in terms of incremental differences in stress components, it is necessary to check the residuals of the original UCM equations with the full stress components. Denoted by RS, the UCM equation residual vector produced by substituting the full stress solution (S2 =S1 + dS)into the original UCM equation at l2, the solution accuracy can be measured by the following residual norm, Res =
B Rs, S\ , BS, S\
(30)
where B, \ denotes a proper inner product, S is the solution vector of the discretized UCM equation. Note that Res in Eq. (30) vanishes at each time step in the actual discretized and linearised equations, i.e. At the n th time step, it vanishes when the RHS is calculated at the n−1 level and LHS at the n th level. When the same n th level solution is used on both the RHS and LHS, however, it may not vanish and its magnitude is a proper measure of convergence in the numerical iterations. The norm defined in Eq. (30) essentially checks solution errors relative to the scaling of the equation, as a result, it is dependent of the arbitrary scaling used in the equation. A fixed time step dt= 0.1 was used for all calculations.
4.1.3. Mesh refinement In addition to satisfying proper convergence criteria, mesh refinement provides another check on the accuracy of finite element results for complex viscoelastic flows. With very few exceptions, numerical results found in the literature [1–9] for the falling sphere problem often have unsatisfactory mesh-convergence or mesh-independence around or beyond Wi=2.0. In the present work, much effort has been devoted to get mesh convergent results for Wi\2.0 Table 2 Solution residuals of the UCM equation as a function of Wi for the three meshes Wi
Mesh M1
Mesh M2
Mesh M3
1.6 1.8 2.0 2.2 2.4 2.6 2.8
7.89×10−6 9.42×10−6 1.23×10−5 1.94×10−5 2.35×10−5 2.81×10−5 3.46×10−5
6.73×10−6 8.48×10−6 1.17×10−5 1.53×10−5 1.96×10−5 2.46×10−5 3.02×10−5
4.96×10−6 6.57×10−6 8.68×10−6 1.07×10−6 1.32×10−5 1.44×10−5 1.52×10−5
66
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Table 3 Predicted drag coefficient K as a function of Wi for the three meshes Wi
Mesh M1
Mesh M2
Mesh M3
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8
4.79525 4.52335 4.33170 4.20482 4.12188 4.06581 4.03604 4.02205 4.02089 4.02533 4.03844 4.05433
4.79602 4.52414 4.33201 4.20562 4.12271 4.06775 4.03726 4.02326 4.02158 4.02666 4.03904 4.05606
4.79668 4.52459 4.33340 4.20758 4.12473 4.06997 4.03957 4.02517 4.02375 4.02838 4.04171 4.05852
We have performed systematic calculations with three well structured quadrilateral meshes: mesh M1 with 108810 total df for pressure, velocity and stress variable, mesh M2 with 149618 total degrees of freedom, and mesh M3 with 196506 total df. Fig. 1 shows the key sections of mesh M1 and M3. The sphere radius is one, and cylinder radius is two. The upstream and downstream tube lengths are 20 and 30, respectively. There is a uniformly high concentration of nodes around the sphere surface and on the centreline near stagnation points. Table 1 gives details of the three meshes, including element sizes on the sphere surface. Note that the element sizes are uniform on the sphere surface, i.e. all elements on the sphere surface have the same size for a given mesh. A direct comparison between the new IDF/OS scheme and the old OS scheme was carried out with the same mesh M1 at Wi= 1.0 on a DEC 300 Model 600 workstation. The new scheme converged in 118 steps with a CPU time of 4 h 11 min, while the old scheme converged in 206 steps with a CPU time of 6 h 43 min.
Fig. 2. Comparison of drag coefficient K among seven different methods: — this work, M3; this work, M1; Lunsmann et al. [2]; + Jun et al. [1]; Fan and Crochet [4]; ×Khomami [3]; 2 Sun et al. [6]; Baaijens et al. [7].
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
67
Fig. 3. Trr and Tzz curves along the centreline and on the sphere surface, Wi =1.6. The range ( − p/2 BX B p/2) corresponds to the sphere surface.
4.2. Drag coefficient predictions Satisfactory mesh-convergent solutions were obtained up to Wi=2.8. The mesh-convergence was measured by the comparison of predictions between all three meshes on drag coefficient and stress distributions near stagnation points as well as on the sphere surface. Beyond Wi= 2.8, although steady state convergence could still be obtained with individual meshes, the mesh-convergence deteriorated rapidly, with mesh M1 predicting significantly lower drag coefficients than mesh M3. For example, at Wi=3.2, mesh M3 predicted K =4.081, while mesh M1 gave K= 3.942, the discrepancy was greater than 3%. Also, the contour plots of stress components are not smooth beyond Wi=2.8. Therefore, no results beyond Wi= 2.8 are reported here for lack of confidence in their accuracy. Table 2 gives the UCM equation residuals as a function of Wi values for all three meshes. In general, the stress equation residual becomes larger as Wi increases, and coarser mesh yields larger residuals. It is interesting to check the same residual of our previous calculations using unstructured meshes [5]. In general, it was found that the residuals were approximately two orders of magnitude larger in the old results than in the present ones. Not surprisingly, the predicted drag coefficients with the unstructured triangle meshes had poor mesh-convergence and were quite inaccurate [5]. The importance of checking stress equation residuals is evident here. Table 3 shows the predicted drag coefficient K for all three meshes, starting at Weissenberg number value 0.6. The predictions on drag coefficient K by the three meshes agree very well with each other. For example, at Wi=2.8, the relative difference in K between M1 and M3 is 0.1%. Our prediction shows that the drag coefficient reaches it minimum value of 4.02 at Wi= 2.2, after which it starts to slowly increase with Wi. It is interesting to note that the
68
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 4. Trr and Tzz curves along the centreline and on the sphere surface, Wi = 2.0.
minimum drag coefficient 4.02 predicted here is the same as the extrapolated asymptotic value reported by Fan and Crochet [4]. Fig. 2 shows a comparison of drag coefficient K among seven different methods, which include those of Jin et al. [1], Lunsmann et al. [2], Khomami [3], Fan and Crochet [4], Sun et al. [6] and Baaijents et al [7], in addition to the present results. Some of the drag data found in the literature do not have a satisfactory mesh-independence, in which case results from the finest mesh are used here for comparison. The latest drag data of Huang et al [8] and Yang and Khomami [9] are not readily available in the literature at the time of writing. As shown in Fig. 2, the agreement among most of the seven different methods is excellent up to Wi= 1.4. At higher Wi numbers, the agreement is less satisfactory, however, there is still a good agreement on the predicted trend.
Fig. 5. Trr and Trr curves along the centreline and on the sphere surface, Wi = 2.2.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
69
Fig. 6. Trr and Tzz curves along the centreline and on the sphere surface, Wi = 2.4.
4.3. Stress distribution cur6es In this section, we show extra-stress components Trr and Tzz near the front and rear stagnation points as well as on the sphere surface. Results at Wi=1.6, 2.0, 2.2, 2.4 and 2.8 from meshes M1 and M3 are compared in Figs. 3–7. The solid lines are M3 results and the dashed lines are M1 results. In general, mesh M2’s results are between those of M1 and M3 and are not shown here. Some useful observations made regarding the stress curves include: 1. The agreement on Trr between M1 and M3 is excellent. In fact, the two sets of curves are virtually indistinguishable from each other at all Weissenberg numbers.
Fig. 7. Trr and Trr curves along the centreline and on the sphere surface, Wi = 2.8.
70
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 8. Contour plots of streamline C, velocity components Ur, Uz, pressure P and the extra stress components Trr, Trz, Tzz and Tuu at Wi = 1.6. For C, Ur and Uz, 20 contour levels are uniformly distributed between the minimum and the maximum values. For P, Trr, Trz, Tzz and Tuu, 10 levels are uniformly distributed between the minimum and zero, and another 10 levels are uniformly distributed between zero and the maximum.
2. The agreement on Tzz between M1 and M3 is also excellent for the most part, the two sets of curves are indistinguishable on the sphere surface ( −p/2BXBp/2) and on the centreline except for a short distance downstream from the rear stagnation point. 3. The mesh-independence of stresses on the sphere surface explains the very good mesh-convergence of drag coefficient predictions as shown in Table 3. 4. The peak values of Trr and Tzz are higher on the upstream sphere surface than in the stress wake of the centreline for all Weissenberg numbers. As shown in Figs. 3–7, the local mesh dependence of Tzz starts at X=2.0, a finite distance away from the rear stagnation point (note: X =p/2 corresponds to the rear stagnation point), which implies no stagnation points singularity is involved here. The Tzz mesh-dependence region is also where the local peak value is in the stress wake. On the one hand, the fact that there is virtually no mesh dependence anywhere else and that there is an obvious local mesh dependence (even at moderate Wi values) after the rear stagnation point, may indicate a possibility of multiple solutions or bifurcations in that region. However, on the other hand, the fact that our finer mesh predicted a lower local peak and that
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
71
there is no evidence of unboundedness in our solutions may still point to a convergent nature of an unique solution. Also, the discrepancy did not increase significantly as Wi increases from 1.6 (Fig. 3) to 2.8 (Fig. 7). One likely explanation on the mesh dependence of Tzz curves in the stress wake is the presence of a thin boundary layer or sharp radial gradient of Tzz along the centreline downstream of the stagnation point, as will be revealed by the contour plots of the stress field, and discussed in Section 4.
4.4. Contour plots of the field solutions Contours of the solution fields at Wi = 1.6, 2.0, 2.2, 2.4 and 2.8 are shown in Figs. 8–12. All of the contour plots have 20 contour levels: for stream function C and velocity components Ur and Uz, the 20 levels are uniformly distributed between the minimum and the maximum values, and for pressure P, and stress components T66, Trz, Tzz and Tuu, we have 10 levels between the minimum (negative) value and zero, and another 10 levels between zero and the maximum (positive). In this way, more details of the solution can be revealed. All of the minimum and maximum values are given in the figures. Mesh M3 results have been used for the contour plots. Two general comments can be made regarding the contour plots:
Fig. 9. Contour plots of solutions at Wi = 2.0. Relative contour level are the same as in Fig. 8.
72
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 10. Contour plots of solutions at Wi =2.2. Relative contour levels are the same as in Fig. 8.
Firstly, the minimum value of Uz is zero at all Wi values, which means that Uz is always non-negative and therefore, there is no recirculation behind the sphere or anywhere else. Secondly, apart from quantitative changes such as the downstream shifts of C, Ur and Uz and the stress boundary layers becoming thinner, the basic patterns and structure of the contours remain unchanged for all Weissenberg numbers. In particular, there is no abrupt qualitative changes in the stress fields as Wi increases. The viscoelastic stress fields at all Wi values are characterised by large gradients and thin boundary layers adjacent to the sphere surface. The stress contours for the highest Wi value 2.8 shown in Fig. 12 appears to have some wiggles, although this did not seem to affect the numerical convergence for the computation at Wi=2.8. One very interesting feature is observed throughout the contours of Tzz. The contour pattern (at all Wi values) downstream of the rear stagnation point shows that the value of Tzz is lower on the centreline than slightly off the centreline, and a thin boundary layer forms along the centreline where Tzz increases then decreases sharply as one moves radially away from the centreline. In other words, the local maximum of Tzz is not on the centreline, but slightly off the centreline. As Wi increases ,this boundary layer becomes thinner. This sharp stress gradient in the thin boundary layer may well be responsible for the local mesh-dependence of Tzz as observed in the last section. The mesh size (dr0.01, see Fig. 1 and Table 1) along the centreline is approximately the same order of magnitude as the thickness of
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
73
the boundary layer, which is of course far from sufficient to adequately resolve the boundary layer, and the solutions are sensitive to the mesh in the high gradient stress awake. It is pointed out that the starting location of this boundary layer is approximately a half sphere radius downstream from the rear stagnation point, which coincides with the starting location of the Tzz mesh-dependence (X: 2.0). The formation of the boundary layer can be explained as follows: the higher value of Tzz on the side surface of the sphere (at X= 0) is sustained and convected downstream, which at downstream is slightly off the centreline and its value is higher than the Tzz value on the centreline. Two different mechanisms are involved here: off the centreline Tzz is convected from the side surface of the sphere where the flow is in mixed shear/elongation and is certainly a strong flow region [15], while on the centreline, Tzz is developed as a result of a pure elongational flow starting from the rear stagnation point where Tzz is zero. The ‘mismatch’ of the centreline with its neighbouring characteristic lines is simply a manifestation of the hyperbolic nature of the UCM equations. The association of the local mesh-dependence of Tzz with the formation of the boundary layer can also explain why coarser mesh predicts higher local peak for Tzz in the stress wake (on the centreline): a coarser mesh can not resolve the sharp radial gradient as well as a finer mesh can, and so it can only represent a smaller radial gradient, which means the Tzz value on the centreline has to be closer to the higher value off the centreline, therefore a coarser mesh predicts a higher Tzz value on the centreline.
Fig. 11. Contour plots of solutions at Wi =2.4. Relative contour levels are the same as in Fig. 8.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
74
Fig. 12. Contour plots of solutions at Wi =2.8. Relative contour levels are the same as in Fig. 8.
At Wi = 2.8, the thickness of this boundary layer appears to be d 0.01, the same as the mesh size in the radial direction in the sphere wake, which indicates the true boundary layer thickness is likely to be even smaller and the possibility of unboundedness of Tzz cannot be ruled out at high Wi values. The coarseness of mesh in the radial direction has a bearing both on accuracy and stability concerning the stress wake boundary layer, and the SUPG scheme does not help here because the radial direction is cross-stream. It is highly desirable to refine the mesh further to reduce the radial size by an order of magnitude to see if the stresses are indeed bounded, but unfortunately this is beyond the computing resources available to the author for the time being. Other components do not seem to form a sharp boundary layer in the wake, at least not at moderate Wi values. For example, a Trr component reduces its magnitude quickly from upstream to downstream, and only a small magnitude is sustained and convected downstream, which is close to the value at the centreline.
5. Conclusions The new algorithm based on incremental difference in variables due to an Wi increment further improves the efficiency and stability of the recently developed OS/SUPG scheme for UCM fluids. Using the new IDF/OS scheme, high resolution mesh-convergent solutions for the
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
75
falling sphere problem have been obtained up to Wi=2.8, beyond which the solutions become significantly mesh-dependent. The drag coefficient reaches its minimum value of 4.02 at approximately Wi= 2.2, after which it starts to increase slowly with Wi. Details of normal stress component distributions along the centreline as well as on the sphere surface have been given. The apparent local mesh-dependence of Tzz in the stress wake, even for globally mesh-independent solutions, gives rise to new questions regarding the cause of convergence loss at high Wi numbers. This local mesh-dependence may not necessarily point to the possibility of non-uniqueness of the solution, but rather it may simply be the result of a sharp boundary layer along the centreline in the wake which has not been adequately resolved in the numerical calculations. Furthermore, numerical work with finer mesh to resolve the stress wake is necessary to study the nature of the boundary layer and the exact cause of the local mesh-dependence of Tzz. The question of whether the stresses in the wake are bounded or not may best be investigated by a combination of boundary layer theory analysis and further numerical studies using meshes which are at least an order of magnitude finer.
Acknowledgements The author wishes to thank the referees for their helpful comments on a draft of this paper.
References [1] H. Jin, N. Phan-Thien, R.I. Tanner, Comput. Mech. 8 (1991) 409 – 422. [2] W.J. Lunsmann, L. Genieser, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 48 (1993) 63. [3] B. Khomami, paper presented at the 8th International Workshop on Numerical Methods in Non-Newtonian Flows, Cape Cod, October 1993. [4] Y. Fan, M.J. Crochet, J. Non-Newtonian Fluid Mech. 57 (1995) 283. [5] X.J. Luo, J. Non-Newtonian Fluid Mech. 63 (1996) 121 – 140. [6] J. Sun, N. Phan-Thien, R.I. Tanner, J. Non-Newtonian Fluid Mech. 95 (1996) 75. [7] F.P T. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, H.E.H. Mejjer, J. Non-Newtonian Fluid Mech. 68 (1997) 173–203. [8] X.-F. Huang, N. Phan-Thien, R.I. Tanner, Viscoelastic flow past a sphere in a tube: an unstructured control volume method, poster presentation at IUTAM 97, July, Sydney. [9] B. Yang, B. Khomami, Modeling of steady viscoelastic flows with hp-adaptive finite element techniques, poster presentation at IUTAM 97, Sydney. [10] R.C. King, M.R. Apelian, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 29 (1988) 147. [11] M. Fortin, A. Fortin, J. Non-Newtonian Fluid Mech. 32 (1989) 295 – 310. [12] X.L. Luo, Int. J. Num. Meth. Fluids 22 (1996) 1189 – 1205. [13] R. Glowinski, Finite element methods for the Numerical Simulation of Imcompressible Viscous Flow in Lectures in Applied Maths, vol. 28, 1991. [14] Fastflo Tutorial Guide, CSIRO Mathematical and information Sciences. More information is available from WWW site:http://www.nag.co.uk/simulation/Fastflo.html [15] R.I. Tanner, Engineering Rheol., Oxford Science Publications, 1985.
.
An incremental difference formulation for viscoelastic flows and high resolution FEM solutions at high Weissenberg numbers X.L. Luo * CSIRO Mathematical and Information Sciences, Locked Bag 17, North Ryde, NSW 2113, Australia Received 26 August 1997; received in revised form 14 January 1998
Abstract An incremental difference formulation (IDF) has been developed to further improve the efficiency and stability of the operator splitting algorithm for viscoelastic flows, taking advantage of the pseudo linearity of the governing equations. Using this new IDF/OS scheme and up to 200000 degrees of freedom, high resolution mesh-convergent solutions for the falling sphere problem are obtained up to Wi = 2.8, with a minimum value of K = 4.02 predicted at Wi=2.2. Beyond Wi =2.8, the solutions become significantly mesh-dependent, although steady convergence with individual meshes are still possible. Well structured quadrilateral elements lead to better accuracy than the entirely unstructured triangle meshes previously used. There is some local mesh-dependence of Tzz on the centreline starting at a finite distance downstream from the rear stagnation point, even for globally mesh-convergent solutions. This local mesh-dependence is most likely due to the presence of a thin boundary layer forming along the centreline in the stress wake, where Tzz exhibits a sharp radial gradient and its local maximum is not on the centreline but slightly off the centreline inside the boundary layer. Further mesh refinement in the radial direction, perhaps reducing the size by an order of magnitude, is needed to study the nature of this boundary layer. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Difference formulation; Operator splitting; Viscoelastic drag; Maxwell fluids; Sphere
1. Introduction In the last few years, several calculations of the sphere-in-tube problem for UCM fluid have reached Weissenberg numbers around or even beyond 2, with some disagreement among different groups on the predicted drag coefficients, or even among different meshes of the same group [1 – 9]. With very few exceptions, mesh-convergence for the sphere problem at Weis* Fax: + 61 0293 253200. 0377-0257/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S0377-0257(98)00063-9
58
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
senberg numbers around and beyond 2.0 has not been as convincingly demonstrated as lower Weissenberg numbers. The detailed and more meaningful comparison of stress distribution curves and contours among different calculations is scarce, owing to lack of published data. One of the purposes of this paper is to provide detailed data for future comparison of solutions at Weissenberg numbers around and beyond the value of 2.0. Jin et al [1] extended the explicitly elliptic momentum equation (EEME) formulation [10] to a class of constitutive equations of the Maxwell type, and obtained results for the sphere problem at high Weissenberg numbers with triangular meshes which are well-structured around the sphere surface but unstructured away from the sphere. Lunsmann et al. [2] were perhaps the first to convincingly demonstrate mesh-convergence of solutions up to Wi=1.6. Fan and Crochet [4] used a high-order finite element method to the sphere problem. The p-convergence up to Weissenberg number 2.0 is verified by means of a residual norm as well as the stress profiles in the stress boundary layer. They estimated an asymptotic value of 4.02 for the drag coefficient at Wi\ 2.0. Sun et al. [6] developed an adaptive viscoelastic stress splitting (AVSS) method and applied it to the sphere problem to obtain solutions up to Weissenberg number of 2.8. More recently, Baaijens et al [7] utilized a stabilised discontinuous Galerkin method, first introduced by Fortin and Fortin [11], for viscoelastic flows. They were able to reach Weissenberg number 2.7 with a moderate number of degrees of freedom, and predicted a minimum drag coefficient of 4.048 at Wi= 2.0, followed by a slow increase in the drag. Huang et al [8] used a new unstructured control volume method to get solutions at high Weissenberg numbers for the sphere problem. Unsmooth stress fields were reported at high Wi values [8]. Yang and Khlomami [9] employed an hp-adaptive finite element technique to the sphere problem to get mesh-convergent results beyond Weissenberg number 2.0. Their results showed that a finer mesh may give a much high Tzz peak than a coarser mesh in the stress wake at high Weissenberg numbers, indicating a possibly unbounded stress wake. The author developed an operator splitting algorithm with a decoupled formulation for transient as well as steady viscoelastic flows [5]. The high efficiency and rebustness of this algorithm enabled high resolution calculations of the sphere problem with over 105 degrees of freedom on a moderate workstation. Unfortunately, our previous calculations used completely unstructured triangle meshes (even on the sphere surface) which led to considerable inaccuracy and poor mesh-convergence at high Weissenberg numbers; This deficiency was due to the inability of unstructured triangle meshes to resolve the thin elastic stress boundary layers on the sphere properly. This inability of unstructured mesh with an uneven element height to resolve thin boundary layers has also been confirmed in the author’s work in turbulence modeling [12]. We point out that the unstructured meshes used by Jin et al [1], Sun et al. [6] and Huang et al [8] are well structured around the sphere surface. The main objective of the present work is to introduce an incremental difference formulation (IDF) for UCM fluids, taking advantage of the pseudo linearity of the momentum equation and UCM equations in a decoupled numerical procedure, and to present comprehensive high resolution results at high Weissenberg numbers for the sphere problem for future comparison with other work. In particular, apart from drag coefficient predictions, detailed stress distribution curves along the centreline and on the sphere surface and contour plots of velocity, pressure and stress fields will be provided, showing all the minimum and maximum values of the solutions.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
59
Much effort in this work involves getting mesh-convergent solutions, rather than trying to reach highest Wi number with a single mesh. Well structured quadrilateral elements with close to 200000 degrees of freedom have been used for high resolution calculations. We found it essential to check solution residuals of the constitutive equations in a decoupled scheme. Our present numerical work predicted a minimum drag coefficient value of 4.02 at Wi =2.2, after which the drag starts to increase slowly with Wi. Mesh-convergent solutions were obtained up to Wi =2.8, beyond which the solution became significantly mesh-dependent, although steady state results with individual meshes were still possible. New findings include the apparent local mesh-dependence of Tzz values on the centreline, which starts at a finite distance downstream from the rear stagnation point. A thin Tzz stress boundary layer forms along the centreline whose position coincides with the local mesh-dependence region. Due to its small thickness this thin boundary layer is likely to be missed by numerical solutions with relatively coarse meshes.
2. Governing equations and the incremental difference formulation Considering the following governing equations for the creeping motion of an incompressible viscoelastic fluid: du + 9p = 9 · T dt
(1)
9 · u = 0,
(2)
where u is velocity, p is pressure and T is the extra-stress tensor. For an upper-convected Maxwell fluid the extra-stress T satisfies the following constitutive equation: l
DT + T = 2h0d, Dt
(3)
where l is the relaxation time, h0 is the constant viscosity, d is the rate of deformation tensor, DT is the upper-convected time derivative: and Dt DT = u · 9T − LT − TLT, Dt
(4)
where the temporal term is omitted for steady flow and L is the velocity gradient tensor. Taking the stress splitting approach, the above system can be reformulated in terms of extra stress S and the Newtonian viscous stress N =2h0d: S T − N.
(5)
The momentum equation and the constitutive equation now become (u − h0Du + 9p =9 · S, (t
(6)
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
60
l
DS DN +S = − l . Dt Dt
(7)
2.1. Incremental difference formulation Eqs. (2), (6) and (7) form a non-linear system of a mixed type, the solution of which in general becomes harder as the Weissenberg number becomes larger. In the past, it is common for a decoupled scheme to start the solution for a higher Weissenberg number from a previously obtained solution at a lower Weissenberg number. In order words the solution at a lower Weissenberg number is often used as the initial guess for the Picard iteration procedure at a higher Wi. This so called incremental approach is usually not suitable for transient problems where one has to start the calculation with a true initial state. In our previous work [5], a transient operator splitting algorithm was developed which does not make use of the solution at a lower Weissenberg number—the time stepping for a given Weissenberg number always starts from zero initial values For steady flow, it is time consuming and unnecessary to do time stepping for all Weissenberg numbers starting from a true initial state. With our operator splitting code, one can of course utilize the incremental approach to save computer time without affect the final steady state solutions. However, in the present work we wish to go one step further: not only we make use of the solution at a lower Weissenberg number, but also we solve only the difference between two solutions at two Weissenberg numbers, taking advantage of the pseudo linearity of the equations in a decoupled procedure. This should improve both the efficiency and stability of the calculations for high Weissenberg numbers, as we will demonstrate in this paper. Let us denote two relaxation times as l1 and l2, assuming l2 \l1. The corresponding field variables at the two relaxation times will be distinguished by subscripts ‘1’ and ‘2’ respectively. The governing equations at l1 are (u1 (8) − h0Du1 + 9p1 = 9 · S1, (t (9) 9 · u1 = 0, DS1 DN1 +S1 = − l1 , (10) l1 Dt Dt and the equations at l2 are (u2 (11) − h0Du2 + 9p2 = 9 · S2, (t (12) 9 · u2 = 0 DS2 DN2 l2 +S2 = − l2 . (13) Dt Dt The momentum equation, continuity equation and the UCM equation are all pseudo linear in that when velocity, pressure and elastic stress are considered separately as in a segregated scheme, there is no quadratic or any higher order terms in the equations. This pseudo linearity is the basis of our new incremental difference scheme. By subtracting both sides of momentum and continuity equations at l1 from those at l2, one gets
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
61
(du − h0Ddu+ 9dp = 9 · dS, (t
(14)
9 · du =0
(15)
where d is used to denote the differences of variables due to an increment in l: du = u2 − u1
dp= p2 − p1
dS=S2 −S1.
(16)
The UCM equations need more than straight forward subtraction because of the cross terms involving tensor multiplication’s of stress and velocity gradient. Define the upper-convective D1 D2 Dd derivative operators , and (without the temporal term) corresponding to velocity Dt Dt Dt fields u1, u2 and du as D1( ) = u1 · 9( )− L1( )− ( ) −L1T, Dt
(17)
D2( ) = u2 · 9( )− L2( )− ( )L2T, Dt
(18)
Dd ( ) =du · 9( ) − dL( )− ( )dLT Dt
(19)
where ( ) represents an arbitrary stress tensor. It follows that D2( ) D1( ) Dd ( ) = + . Dt Dt Dt
(20)
Splitting S2 into S1 and dS and moving all S1 terms to the RHS, Eq. (13) becomes l2
D2dS D2S1 D2N2 + dS= − l2 − S1 − l2 Dt Dt Dt
(21)
The RHS of Eq. (21) can be expressed entirely in terms of increments of either dS=S2 −S1, or dN = N2 − N1 or dl =l2 − l1. Making use of Eq. (20), the RHS of Eq. (21) can be rewritten as RHS = − l2
n
D1S1 Dd S1 D(N1 +dN) Dd N2 − S1 − l2 −l2 + . Dt Dt Dt Dt
(22)
Now from Eq. (10), we get − l2
D1S1 D1N1 l2 −l2 = S1 Dt Dt l1
(23)
Substituting Eq. (23) into Eq. (22), we get RHS =
dl Dd S1 Dd N2 D1d N S1 − l2 . + + l1 Dt Dt Dt
Denoting S1 + N2 as G and substituting Eq. (24) into Eq. (21), one gets the final form
(24)
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
62
l2
D2dS dl DdG D1dN + dS= S1 − l2 + . Dt l1 Dt Dt
(25)
To sum up, Eqs. (14), (15) and (25) complete our formulation in terms of variable difference due to an increment in l. Note the operator Dd/Dt contains the difference in velocity and its gradient. A few necessary remarks include: Remark 1. All of the equations in the incremental difference formulation are exact, there is no approximation. Remark 2. All of the unknowns are in terms of difference between solutions at two Weissenberg numbers which significantly reduces the maximum values of the unknown stress components and their gradients in both the momentum and the UCM equations. For example, at l1 = 2.0 and l2 = 2.2, numerical results in the present work show that the maximum value in dSzz is only 19.9, while the maximum value of the full elastic stress component Szz at Weissenberg number 2.2 is 145.1, more than seven times larger than the difference. Consequently, the gradient 9dSzz is also much smaller than the full value 9Szz Remark 3. The incremental difference momentum and continuity equations (Eqs. (14) and (15)) are of exactly the same form as the original equations (Eqs. (1) and (2)), and the reformulated UCM equation takes the same form in LHS as the original, but with a few extra terms on the RHS. Therefore, the methods of solution for these equations can remain virtually unchanged from those for the full variable equations. In any decoupled scheme, the operators D1 D2 Dd and can all be evaluated readily from u1 and du and the same is true for N2 and dN, Dt’ Dt Dt given u1 and du. Remark 4. The Dirichlet boundary conditions for velocity difference du are simply zero everywhere in the new formulation. The boundary condition on dS can also be readily worked out, which is simply zero for the sphere problem.
3. IDF/OS scheme for viscoelastic flows As the governing equations in an IDF take the same form as the original full variable ones, IDF can be readily combined with our previously developed operator splitting (OS) algorithm [5] to produce a more efficient and more stable scheme-IDF/OS. Details regarding the OS algorithm can be found in [5]. Briefly, a transient algorithm is adopted to solve the nonlinear system of momentum, continuity and UCM equations in a segregated manner. The pre-conditioned conjugate gradient (PGG) algorithm of Glowinski [13] is used to solve the Quasi Stokes (QS) equation with viscoelastic stress contributions, which only involves solving a sequence of scalar Dirichlet problems associated with the Helmholtz operator and Neuman problems associated with the Laplacian operator, plus a few scalar intergrations. The matrix associated with the Helmholtz operator is always well conditioned (as the time step can be small), thus iterative solvers such as CG, SOR etc. will converge fast. For constant viscosity h0 and fixed time step Dt, all of the matrices in the PCG process remain the same for all time steps. As a result, the assembly and LU factorization for a direct solver or the incomplete factorization for an iterative solver need only be carried out once for all time steps, which saves on computing time considerably.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
63
In addition to fast convergence and high computing efficiency, the PCG scheme directly checks on the divergence of the velocity field as a convergence criterion. As a result, the incompressibility constraint is always well satisfied and the pressure solution obtained is usually more accurate than some other conventional methods such as the penalty type. A consistent SUPG method is used for the UCM equations, along with generalized minimal residual (GMRES) iterative solver for the asymmetric systems. Integration by parts is used to avoid the second derivatives in Eq. (25), and the extra surface integration term concerning the boundary flux of the deformation rate can be readily calculated.
4. High resolution solutions of the falling sphere problem After the IDF/OS scheme was implemented by the author using a high level FEM language Fasttalk [14], the new code was applied to the standard sphere-in-tube bench mark problem. As discussed in Section 1, our previous calculations with the OS scheme [5] used totally unstructured triangle meshes, which led to considerable inaccuracy at high Weissenberg numbers due to the inability of the unstructured mesh to resolve the thin elastic stress boundary layers. In the present work, we only used well-structured quadrilateral elements. A direct comparison of the old OS scheme with the new IDF/OS scheme was carried out with the same quadrilateral mesh, in order to demonstrate the higher efficiency and the stability of the latter. A well-structured triangle mesh proves to be equally as good as the quadrilateral ones. The improved efficiency of our new formulation allowed us to attack this difficult problem with very large degrees of freedom (up to 200000) on a moderate workstation. Our aim is to present detailed results on the stress as well as velocity and pressure fields to facilitate careful comparisons with other works in the future.
4.1. Some numerical details For the sphere-in-tube problem, the global quantity of interest is the non-dimensional drag coefficient: K=
F 6ph0Ua
(26)
as a function of the Weissenberg number Wi=
lU , a
(27)
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, l is the relaxation time of the UCM fluid and h0 is the viscosity. In the standard benchmark problem, the ratio of sphere radius over the cylinder radius is a/R = 0.5. As in previous OS schemes, a quadratic interpretation is used for velocity, the extra stresses, and linear for the pressure.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
64
4.1.1. Solution strategy We initially obtained solution at Wi= 0.6 with our previous conventional OS/SUPG formulation, i.e. the full velocity and stress component formulation. Consequently, we let l1 =0.6 and dl = 0.2 and used our new IDF/OS scheme to advance the solutions to higher Weissenberg numbers. The l1 was of course changed to l1 +dl after each convergent solution, ready for the next increment in Wi. The increment dl was fixed at 0.2 for all calculations. 4.1.2. Con6ergence criteria—checking equation residuals We believe all numerical work should report the convergence criteria used, which unfortunately is not always the case in the literature. In our previous work [5], the following convergence criteria for steady state solutions has been used: Convergence to steady state was considered to be reached if the following two conditions were satisfied simultaneously, u n + 1 − u n B eu u n+1
(28)
K n+1 − K n B eK, K n+1
(29)
and
where K n + 1 and K n are drag coefficients, and un + 1 and un are the velocity components at two consecutive time steps. The norms are the usual L2 norm. We have set eu =0.5×10 − 3 and eK = 10 − 5, respectively in our previous work. With the new formation, we found an even tighter convergence test can be applied to the momentum equation—we let eu = 10 − 4. Furthermore, due to the numerical solution of the UCM
Fig. 1. Partial view of mesh M1 and M3.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
65
Table 1 Mesh details
M1 M2 M3
No. nodes
No. ele.
df
Size (ds/dr)
17205 23653 31061
5580 7700 10140
108 810 149 618 196 506
0.0175/0.00202 0.0142/0.00166 0.0113/0.00128
The size ds is the arc length of elements on the sphere surface and dr is the radial length.
equation being in terms of incremental differences in stress components, it is necessary to check the residuals of the original UCM equations with the full stress components. Denoted by RS, the UCM equation residual vector produced by substituting the full stress solution (S2 =S1 + dS)into the original UCM equation at l2, the solution accuracy can be measured by the following residual norm, Res =
B Rs, S\ , BS, S\
(30)
where B, \ denotes a proper inner product, S is the solution vector of the discretized UCM equation. Note that Res in Eq. (30) vanishes at each time step in the actual discretized and linearised equations, i.e. At the n th time step, it vanishes when the RHS is calculated at the n−1 level and LHS at the n th level. When the same n th level solution is used on both the RHS and LHS, however, it may not vanish and its magnitude is a proper measure of convergence in the numerical iterations. The norm defined in Eq. (30) essentially checks solution errors relative to the scaling of the equation, as a result, it is dependent of the arbitrary scaling used in the equation. A fixed time step dt= 0.1 was used for all calculations.
4.1.3. Mesh refinement In addition to satisfying proper convergence criteria, mesh refinement provides another check on the accuracy of finite element results for complex viscoelastic flows. With very few exceptions, numerical results found in the literature [1–9] for the falling sphere problem often have unsatisfactory mesh-convergence or mesh-independence around or beyond Wi=2.0. In the present work, much effort has been devoted to get mesh convergent results for Wi\2.0 Table 2 Solution residuals of the UCM equation as a function of Wi for the three meshes Wi
Mesh M1
Mesh M2
Mesh M3
1.6 1.8 2.0 2.2 2.4 2.6 2.8
7.89×10−6 9.42×10−6 1.23×10−5 1.94×10−5 2.35×10−5 2.81×10−5 3.46×10−5
6.73×10−6 8.48×10−6 1.17×10−5 1.53×10−5 1.96×10−5 2.46×10−5 3.02×10−5
4.96×10−6 6.57×10−6 8.68×10−6 1.07×10−6 1.32×10−5 1.44×10−5 1.52×10−5
66
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Table 3 Predicted drag coefficient K as a function of Wi for the three meshes Wi
Mesh M1
Mesh M2
Mesh M3
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8
4.79525 4.52335 4.33170 4.20482 4.12188 4.06581 4.03604 4.02205 4.02089 4.02533 4.03844 4.05433
4.79602 4.52414 4.33201 4.20562 4.12271 4.06775 4.03726 4.02326 4.02158 4.02666 4.03904 4.05606
4.79668 4.52459 4.33340 4.20758 4.12473 4.06997 4.03957 4.02517 4.02375 4.02838 4.04171 4.05852
We have performed systematic calculations with three well structured quadrilateral meshes: mesh M1 with 108810 total df for pressure, velocity and stress variable, mesh M2 with 149618 total degrees of freedom, and mesh M3 with 196506 total df. Fig. 1 shows the key sections of mesh M1 and M3. The sphere radius is one, and cylinder radius is two. The upstream and downstream tube lengths are 20 and 30, respectively. There is a uniformly high concentration of nodes around the sphere surface and on the centreline near stagnation points. Table 1 gives details of the three meshes, including element sizes on the sphere surface. Note that the element sizes are uniform on the sphere surface, i.e. all elements on the sphere surface have the same size for a given mesh. A direct comparison between the new IDF/OS scheme and the old OS scheme was carried out with the same mesh M1 at Wi= 1.0 on a DEC 300 Model 600 workstation. The new scheme converged in 118 steps with a CPU time of 4 h 11 min, while the old scheme converged in 206 steps with a CPU time of 6 h 43 min.
Fig. 2. Comparison of drag coefficient K among seven different methods: — this work, M3; this work, M1; Lunsmann et al. [2]; + Jun et al. [1]; Fan and Crochet [4]; ×Khomami [3]; 2 Sun et al. [6]; Baaijens et al. [7].
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
67
Fig. 3. Trr and Tzz curves along the centreline and on the sphere surface, Wi =1.6. The range ( − p/2 BX B p/2) corresponds to the sphere surface.
4.2. Drag coefficient predictions Satisfactory mesh-convergent solutions were obtained up to Wi=2.8. The mesh-convergence was measured by the comparison of predictions between all three meshes on drag coefficient and stress distributions near stagnation points as well as on the sphere surface. Beyond Wi= 2.8, although steady state convergence could still be obtained with individual meshes, the mesh-convergence deteriorated rapidly, with mesh M1 predicting significantly lower drag coefficients than mesh M3. For example, at Wi=3.2, mesh M3 predicted K =4.081, while mesh M1 gave K= 3.942, the discrepancy was greater than 3%. Also, the contour plots of stress components are not smooth beyond Wi=2.8. Therefore, no results beyond Wi= 2.8 are reported here for lack of confidence in their accuracy. Table 2 gives the UCM equation residuals as a function of Wi values for all three meshes. In general, the stress equation residual becomes larger as Wi increases, and coarser mesh yields larger residuals. It is interesting to check the same residual of our previous calculations using unstructured meshes [5]. In general, it was found that the residuals were approximately two orders of magnitude larger in the old results than in the present ones. Not surprisingly, the predicted drag coefficients with the unstructured triangle meshes had poor mesh-convergence and were quite inaccurate [5]. The importance of checking stress equation residuals is evident here. Table 3 shows the predicted drag coefficient K for all three meshes, starting at Weissenberg number value 0.6. The predictions on drag coefficient K by the three meshes agree very well with each other. For example, at Wi=2.8, the relative difference in K between M1 and M3 is 0.1%. Our prediction shows that the drag coefficient reaches it minimum value of 4.02 at Wi= 2.2, after which it starts to slowly increase with Wi. It is interesting to note that the
68
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 4. Trr and Tzz curves along the centreline and on the sphere surface, Wi = 2.0.
minimum drag coefficient 4.02 predicted here is the same as the extrapolated asymptotic value reported by Fan and Crochet [4]. Fig. 2 shows a comparison of drag coefficient K among seven different methods, which include those of Jin et al. [1], Lunsmann et al. [2], Khomami [3], Fan and Crochet [4], Sun et al. [6] and Baaijents et al [7], in addition to the present results. Some of the drag data found in the literature do not have a satisfactory mesh-independence, in which case results from the finest mesh are used here for comparison. The latest drag data of Huang et al [8] and Yang and Khomami [9] are not readily available in the literature at the time of writing. As shown in Fig. 2, the agreement among most of the seven different methods is excellent up to Wi= 1.4. At higher Wi numbers, the agreement is less satisfactory, however, there is still a good agreement on the predicted trend.
Fig. 5. Trr and Trr curves along the centreline and on the sphere surface, Wi = 2.2.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
69
Fig. 6. Trr and Tzz curves along the centreline and on the sphere surface, Wi = 2.4.
4.3. Stress distribution cur6es In this section, we show extra-stress components Trr and Tzz near the front and rear stagnation points as well as on the sphere surface. Results at Wi=1.6, 2.0, 2.2, 2.4 and 2.8 from meshes M1 and M3 are compared in Figs. 3–7. The solid lines are M3 results and the dashed lines are M1 results. In general, mesh M2’s results are between those of M1 and M3 and are not shown here. Some useful observations made regarding the stress curves include: 1. The agreement on Trr between M1 and M3 is excellent. In fact, the two sets of curves are virtually indistinguishable from each other at all Weissenberg numbers.
Fig. 7. Trr and Trr curves along the centreline and on the sphere surface, Wi = 2.8.
70
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 8. Contour plots of streamline C, velocity components Ur, Uz, pressure P and the extra stress components Trr, Trz, Tzz and Tuu at Wi = 1.6. For C, Ur and Uz, 20 contour levels are uniformly distributed between the minimum and the maximum values. For P, Trr, Trz, Tzz and Tuu, 10 levels are uniformly distributed between the minimum and zero, and another 10 levels are uniformly distributed between zero and the maximum.
2. The agreement on Tzz between M1 and M3 is also excellent for the most part, the two sets of curves are indistinguishable on the sphere surface ( −p/2BXBp/2) and on the centreline except for a short distance downstream from the rear stagnation point. 3. The mesh-independence of stresses on the sphere surface explains the very good mesh-convergence of drag coefficient predictions as shown in Table 3. 4. The peak values of Trr and Tzz are higher on the upstream sphere surface than in the stress wake of the centreline for all Weissenberg numbers. As shown in Figs. 3–7, the local mesh dependence of Tzz starts at X=2.0, a finite distance away from the rear stagnation point (note: X =p/2 corresponds to the rear stagnation point), which implies no stagnation points singularity is involved here. The Tzz mesh-dependence region is also where the local peak value is in the stress wake. On the one hand, the fact that there is virtually no mesh dependence anywhere else and that there is an obvious local mesh dependence (even at moderate Wi values) after the rear stagnation point, may indicate a possibility of multiple solutions or bifurcations in that region. However, on the other hand, the fact that our finer mesh predicted a lower local peak and that
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
71
there is no evidence of unboundedness in our solutions may still point to a convergent nature of an unique solution. Also, the discrepancy did not increase significantly as Wi increases from 1.6 (Fig. 3) to 2.8 (Fig. 7). One likely explanation on the mesh dependence of Tzz curves in the stress wake is the presence of a thin boundary layer or sharp radial gradient of Tzz along the centreline downstream of the stagnation point, as will be revealed by the contour plots of the stress field, and discussed in Section 4.
4.4. Contour plots of the field solutions Contours of the solution fields at Wi = 1.6, 2.0, 2.2, 2.4 and 2.8 are shown in Figs. 8–12. All of the contour plots have 20 contour levels: for stream function C and velocity components Ur and Uz, the 20 levels are uniformly distributed between the minimum and the maximum values, and for pressure P, and stress components T66, Trz, Tzz and Tuu, we have 10 levels between the minimum (negative) value and zero, and another 10 levels between zero and the maximum (positive). In this way, more details of the solution can be revealed. All of the minimum and maximum values are given in the figures. Mesh M3 results have been used for the contour plots. Two general comments can be made regarding the contour plots:
Fig. 9. Contour plots of solutions at Wi = 2.0. Relative contour level are the same as in Fig. 8.
72
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
Fig. 10. Contour plots of solutions at Wi =2.2. Relative contour levels are the same as in Fig. 8.
Firstly, the minimum value of Uz is zero at all Wi values, which means that Uz is always non-negative and therefore, there is no recirculation behind the sphere or anywhere else. Secondly, apart from quantitative changes such as the downstream shifts of C, Ur and Uz and the stress boundary layers becoming thinner, the basic patterns and structure of the contours remain unchanged for all Weissenberg numbers. In particular, there is no abrupt qualitative changes in the stress fields as Wi increases. The viscoelastic stress fields at all Wi values are characterised by large gradients and thin boundary layers adjacent to the sphere surface. The stress contours for the highest Wi value 2.8 shown in Fig. 12 appears to have some wiggles, although this did not seem to affect the numerical convergence for the computation at Wi=2.8. One very interesting feature is observed throughout the contours of Tzz. The contour pattern (at all Wi values) downstream of the rear stagnation point shows that the value of Tzz is lower on the centreline than slightly off the centreline, and a thin boundary layer forms along the centreline where Tzz increases then decreases sharply as one moves radially away from the centreline. In other words, the local maximum of Tzz is not on the centreline, but slightly off the centreline. As Wi increases ,this boundary layer becomes thinner. This sharp stress gradient in the thin boundary layer may well be responsible for the local mesh-dependence of Tzz as observed in the last section. The mesh size (dr0.01, see Fig. 1 and Table 1) along the centreline is approximately the same order of magnitude as the thickness of
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
73
the boundary layer, which is of course far from sufficient to adequately resolve the boundary layer, and the solutions are sensitive to the mesh in the high gradient stress awake. It is pointed out that the starting location of this boundary layer is approximately a half sphere radius downstream from the rear stagnation point, which coincides with the starting location of the Tzz mesh-dependence (X: 2.0). The formation of the boundary layer can be explained as follows: the higher value of Tzz on the side surface of the sphere (at X= 0) is sustained and convected downstream, which at downstream is slightly off the centreline and its value is higher than the Tzz value on the centreline. Two different mechanisms are involved here: off the centreline Tzz is convected from the side surface of the sphere where the flow is in mixed shear/elongation and is certainly a strong flow region [15], while on the centreline, Tzz is developed as a result of a pure elongational flow starting from the rear stagnation point where Tzz is zero. The ‘mismatch’ of the centreline with its neighbouring characteristic lines is simply a manifestation of the hyperbolic nature of the UCM equations. The association of the local mesh-dependence of Tzz with the formation of the boundary layer can also explain why coarser mesh predicts higher local peak for Tzz in the stress wake (on the centreline): a coarser mesh can not resolve the sharp radial gradient as well as a finer mesh can, and so it can only represent a smaller radial gradient, which means the Tzz value on the centreline has to be closer to the higher value off the centreline, therefore a coarser mesh predicts a higher Tzz value on the centreline.
Fig. 11. Contour plots of solutions at Wi =2.4. Relative contour levels are the same as in Fig. 8.
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
74
Fig. 12. Contour plots of solutions at Wi =2.8. Relative contour levels are the same as in Fig. 8.
At Wi = 2.8, the thickness of this boundary layer appears to be d 0.01, the same as the mesh size in the radial direction in the sphere wake, which indicates the true boundary layer thickness is likely to be even smaller and the possibility of unboundedness of Tzz cannot be ruled out at high Wi values. The coarseness of mesh in the radial direction has a bearing both on accuracy and stability concerning the stress wake boundary layer, and the SUPG scheme does not help here because the radial direction is cross-stream. It is highly desirable to refine the mesh further to reduce the radial size by an order of magnitude to see if the stresses are indeed bounded, but unfortunately this is beyond the computing resources available to the author for the time being. Other components do not seem to form a sharp boundary layer in the wake, at least not at moderate Wi values. For example, a Trr component reduces its magnitude quickly from upstream to downstream, and only a small magnitude is sustained and convected downstream, which is close to the value at the centreline.
5. Conclusions The new algorithm based on incremental difference in variables due to an Wi increment further improves the efficiency and stability of the recently developed OS/SUPG scheme for UCM fluids. Using the new IDF/OS scheme, high resolution mesh-convergent solutions for the
X.L. Luo / J. Non-Newtonian Fluid Mech. 79 (1998) 57–75
75
falling sphere problem have been obtained up to Wi=2.8, beyond which the solutions become significantly mesh-dependent. The drag coefficient reaches its minimum value of 4.02 at approximately Wi= 2.2, after which it starts to increase slowly with Wi. Details of normal stress component distributions along the centreline as well as on the sphere surface have been given. The apparent local mesh-dependence of Tzz in the stress wake, even for globally mesh-independent solutions, gives rise to new questions regarding the cause of convergence loss at high Wi numbers. This local mesh-dependence may not necessarily point to the possibility of non-uniqueness of the solution, but rather it may simply be the result of a sharp boundary layer along the centreline in the wake which has not been adequately resolved in the numerical calculations. Furthermore, numerical work with finer mesh to resolve the stress wake is necessary to study the nature of the boundary layer and the exact cause of the local mesh-dependence of Tzz. The question of whether the stresses in the wake are bounded or not may best be investigated by a combination of boundary layer theory analysis and further numerical studies using meshes which are at least an order of magnitude finer.
Acknowledgements The author wishes to thank the referees for their helpful comments on a draft of this paper.
References [1] H. Jin, N. Phan-Thien, R.I. Tanner, Comput. Mech. 8 (1991) 409 – 422. [2] W.J. Lunsmann, L. Genieser, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 48 (1993) 63. [3] B. Khomami, paper presented at the 8th International Workshop on Numerical Methods in Non-Newtonian Flows, Cape Cod, October 1993. [4] Y. Fan, M.J. Crochet, J. Non-Newtonian Fluid Mech. 57 (1995) 283. [5] X.J. Luo, J. Non-Newtonian Fluid Mech. 63 (1996) 121 – 140. [6] J. Sun, N. Phan-Thien, R.I. Tanner, J. Non-Newtonian Fluid Mech. 95 (1996) 75. [7] F.P T. Baaijens, S.H.A. Selen, H.P.W. Baaijens, G.W.M. Peters, H.E.H. Mejjer, J. Non-Newtonian Fluid Mech. 68 (1997) 173–203. [8] X.-F. Huang, N. Phan-Thien, R.I. Tanner, Viscoelastic flow past a sphere in a tube: an unstructured control volume method, poster presentation at IUTAM 97, July, Sydney. [9] B. Yang, B. Khomami, Modeling of steady viscoelastic flows with hp-adaptive finite element techniques, poster presentation at IUTAM 97, Sydney. [10] R.C. King, M.R. Apelian, R.C. Armstrong, R.A. Brown, J. Non-Newtonian Fluid Mech. 29 (1988) 147. [11] M. Fortin, A. Fortin, J. Non-Newtonian Fluid Mech. 32 (1989) 295 – 310. [12] X.L. Luo, Int. J. Num. Meth. Fluids 22 (1996) 1189 – 1205. [13] R. Glowinski, Finite element methods for the Numerical Simulation of Imcompressible Viscous Flow in Lectures in Applied Maths, vol. 28, 1991. [14] Fastflo Tutorial Guide, CSIRO Mathematical and information Sciences. More information is available from WWW site:http://www.nag.co.uk/simulation/Fastflo.html [15] R.I. Tanner, Engineering Rheol., Oxford Science Publications, 1985.
.