Numerical simulation of Weissenberg phenomena – the rod-climbing of viscoelastic fluids

Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

www.elsevier.com/locate/cma

Numerical simulation of Weissenberg phenomena ± the rod-climbing of viscoelastic ¯uids X.-L. Luo * CSIRO Mathematical and Information Sciences, Locked Bag 17, North Ryde, NSW 2113, Australia Received 30 October 1998

Abstract A ®nite element simulation of the Weissenberg e€ect, i.e., the rod-climbing of viscoelastic ¯uids, is presented. The ¯ow features axisymmetric swirling, free surface, gravity, surface tension, centrifugal force and all six viscoelastic stresses. An operator splitting algorithm is employed to solve the highly non-linear system of equations governing the upper-convected Maxwell (UCM) or the PhanThien±Tanner (PTT) ¯uid. Normal displacement of free surface is tracked explicitly in time and an elaborate orthogonal trajectory scheme applicable to unstructured as well as structured grids is developed to update the mesh. Comparison of numerical and experimental results shows good agreement. Ó 1999 Elsevier Science S.A. All rights reserved.

1. Introduction It is well-known that ¯uids with complex structure, such as macromolecular solutions and melts, behave in unexpected ways which cannot be described by Newton's viscous law. The Weissenberg e€ect [1±5] is one of the most interesting phenomena exhibited by non-Newtonian ¯uids: certain viscoelastic ¯uids in a cylindrical vessel will climb up a rotating rod against centrifugal force and gravity. The rod-climbing is one of the many second-order e€ects associated with the inequality of normal stresses in shear ¯ow [2,3]. It was suggested [2,4] that the normal stress acts like a hoop around the rod and forces the ¯uid inwards against the centrifugal force and upwards against the gravitational force. Numerical, analytical and experimental studies of this phenomena will improve the understanding of the rheological properties of non-Newtonian ¯uids and enable us to better predict their behaviours in numerous natural and industrial non-Newtonian ¯ows. Some discussion on the rheometrical applications of this ¯ow can be found in Refs. [2,4,5]. Several analytical studies on rod-climbing of viscoelastic ¯uids have been done using perturbation theory valid only for small angular speeds [6±11]. These studies have focused on the relationship between viscoelastic properties and rod-climbing quantities such as the climbing height and shape. One of the signi®cant ®ndings is that from the knowledge of the climbing one can determine the values of some important rheological parameters of simple ¯uids. Joseph et al. [7,12] have done some extensive experimental as well as analytical investigations into rod-climbing of several ¯uids, providing some valuable data for validating analytical and numerical results. Numerical simulation of rod-climbing of viscoelastic ¯uids is dicult in several aspects: (1) free surface with large deformation but without ®xed end points is present, some popular free surface tracking techniques, such as the ones used in extrusion simulations, are not e€ective here; (2) axisymmetric swirling ¯ow

*

Tel.: +61-02-932-53265; fax: +61-0293-25-3200. E-mail address: [email protected] (X.-L. Luo).

0045-7825/99/$ - see front matter Ó 1999 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 1 7 5 - 9

394

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

con®guration means all six stress components are non-zero, in addition to the extra angular velocity component, which makes the problem computationally more demanding than without swirling; (3) the general high Weissenberg number problem [13] may occur when the angular velocity is high and other factors such as surface tension, gravity and centrifugal force also add to the complexity of the problem; (4) since the main focus here is the secondary ¯ow, insucient accuracy will fail to capture correctly the nonlinear behaviour of this phenomenon and double precision is the preferred computing mode. Debbaut and Hocq [14] have made perhaps the ®rst attempt to numerically solve the full system of equations for the rod-climbing problem, using an algorithm developed by Marchal and Crochet [15] in which 4  4 sub-element interpolations for stress tensor were combined with a streamline-upwinding (SU) technique. Interesting solutions were obtained at relatively low values of Weissenberg (Wi) number. An implicit technique similar to that developed by Kistler and Scriven [16] was adopted for predicting the free surface location with structured quadrilateral mesh. They found at moderate Wi the free surface experiences large deformations and also shows a loss of smoothness in its ®nite element discretization. The present work explores new ways to overcome the diculties encountered in rod-climbing simulation and attempts to improve numerical accuracy without sacri®cing much on eciency. An elaborate mesh updating scheme is developed to minimize mesh distortion by making use of orthogonal or nearly orthogonal trajectory functions. Both the structured quadrilateral mesh and the unstructured triangular mesh will be used in the simulation to check the mesh-independence of the solutions. Section 2 presents governing equations for the axisymmetric swirling ¯ow of viscoelastic ¯uids, particularly of the upper-convected Maxwell (UCM) model and the Phan-Thien±Tanner (PTT) model and describes brie¯y the operator splitting algorithm. Section 3 is devoted to free surface treatment, with a detailed discussion on orthogonal trajectory approach for both structured quadrilateral and unstructured triangular meshes. Section 4 shows numerical results in comparison with perturbation theory and experimental measurements. Conclusions are given in Section 5. 2. Governing equations and operator splitting The motion of an incompressible viscoelastic ¯uid is governed by the following equations,   ou ‡ …u
r†u ÿ gDu ‡ rp ˆ qg ‡ r
s; q ot r
u ˆ 0;

…1† …2†

where g is the sum of solvent and molecular viscosities, u the velocity, p the pressure, g the gravity and s the extra viscoelastic stress tensor de®ned by s  r ‡ pI ÿ 2gd;

…3†

where r is the total stress tensor. For UCM ¯uid the extra viscoelastic stress satis®es the following constitutive equation,     os od T T ‡ u
rs ÿ Ls ÿ sL ‡ s ˆ ÿ2kgm ‡ u
rd ÿ Ld ÿ dL ; …4† k ot ot where k is the relaxation time, gm the molecular viscosity, L ˆ ruT the velocity gradient tensor, d ˆ …L ‡ LT †=2 is the rate of deformation tensor. The PTT model can also be expressed in terms of s,   os T ‡ u
rs ÿ …L ÿ nd†s ÿ s…L ÿ nd† ‡ ‰1 ‡ …k=gm † tr sŠs ‡ …2k tr s†d k ot   od ˆ ÿ2kgm …5† ‡ u
rd ÿ …L ÿ nd†d ÿ d…L ÿ nd†T ; ot where tr s is the trace of s, n is a material parameter controlling the ratio of second to ®rst normal stress di€erence, and  is related to elongational viscosity.

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

395

In the case of axisymmetric swirling ¯ow, we have altogether 10 unknown variables: three velocity components ur , uz and uh , pressure p and six extra viscoelastic stresses srr , srz , szz , shh , srh and szh . Unfortunately the full component form of the Maxwell model or the PTT model for axisymmetric swirling ¯ow cannot be readily found in text books and it is rather tedious to work them out from the tensor notations. The axisymmetric swirling condition means o=oh ˆ 0 holds for any scalar variables but …o^r=oh† ˆ h^ and ^ …oh=oh† ˆ ÿ^r, where ^r and h^ are the two unit directional vectors in radial and circumferential coordinates. From the above conditions one can work out the components for the following tensors involving gradients of velocity and stress, 0 1 u
rsrh ‡ uh …srrrÿshh † u
rsrz ÿ uhrshz u
rsrr ÿ 2uhrsrh B C uh …srr ÿshh † 2uh srh uh srz C; u
rs ˆ B …6† u
rs ‡ u
rs ‡ hh hz @ u
rshr ‡ r r r A u
rszr ÿ uhrszh 0 our or

B ou h L ˆ ruT ˆ B @ or ouz or

u
rszh ‡ uhrszr

ÿ urh

our oz

ur r

ouh oz

0

ouz oz

u
rszz

1 C C: A

…7†

After substituting Eqs. (6) and (7) into the Maxwell or PTT equation and some tedious algebraic operations one gets the full component form for each constitutive equation for axisymmetric swirling ¯ow. The equations are relatively lengthy compared with the no-swirling case. For example, the Maxwell equation for srh reads   osrh uh ouz ouh our ouh ‡ u
rsrh ‡ srr ‡ srh ÿ srr ÿ szh ÿ szr ‡ srh k ot r oz or oz oz         o ouh uh ouh uh uh our ouz ouh ouz ÿ ‡u
r ÿ ‡ 2 ÿ ‡ ˆ ÿkgm r or r r or oz or oz ot or   our ouh ouh our ouz ouh ‡2 ‡ : …8† ‡ kgm 2 or or oz oz or oz The component form for the momentum and continuity equations for axisymmetric swirling ¯ow are readily available from text books. In this work the surface tension rs has to be included rs ˆ

C ; Rp

where C is the surface tension coecient and Rp the principal radius of curvature in the free surface. Since inertia cannot be omitted here the full operator splitting algorithm as described in our early work [17] is adopted in this work. Brie¯y the constitutive equations are solved by a consistent Streamline-Upwind/Petrov-Galerkin (SUPG) scheme [21] in a segregated manner with an inner loop of iterations and the momentum and continuity equations are solved in three fractional time steps by the h-scheme of Glowinski [18±20]: First fractional step for the nth time step: Find un‡h 2 V, pn‡h 2 L2 …X†, 8U 2 V, 8q 2 L2 …X†,  Z Z Z  un‡h ÿ un dt n‡h n‡h dX ˆ qgU dX ‡ ‡ agrU
ru ÿ rUp Uq C U ds hDt ds oX Z h i ‡ ÿ bgrU
run ÿ grU
…run †T ÿ Uq…un
r†un ÿ rU
sn‡0:5 dX; Z

qr
un‡h ˆ 0:

…9† …10†

396

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Second fractional step: Find un‡1ÿh 2 V,8U 2 V,  Z  Z Z un‡1ÿh ÿ un‡h dt n‡1ÿh n‡h n‡1ÿh Uq dX ˆ qgU dX ‡ ‡ bgrU
ru ‡ U…u
r†u C Uds …1 ÿ 2h†Dt ds oX Z h i T ‡ ÿ agrU
run‡h ÿ grU
…run † ‡ rUpn‡h ÿ rU
sn‡0:5 dX: Third fractional step: Find un‡1 2 V, pn‡1 2 L2 …X†, 8U 2 V, 8q 2 L2 …X†,  Z  Z Z un‡1 ÿ un‡1ÿh dt n‡1 n‡1 Uq dX ˆ qgU dX ‡ ‡ agrU
ru ÿ rUp C U ds hDt ds oX Z h i n‡1ÿh n T n‡1ÿh n‡1ÿh ÿ grU
…ru † ÿ U…u
r†u ÿ rU
sn‡0:5 dX; ‡ ÿ bgrU
ru Z qr
un‡1 dX ˆ 0;

…11†

…12† …13†

where t is the unit tangent vector, the tangent derivative of which gives the curvature for surface tension calculation and a, b and h are the operator splitting parameters satisfying a ‡ b ˆ 1 and 0 < h < 1. The sub-problems at ®rst and third fractional steps can be made identical and can be solved by the preconditioned conjugate gradient scheme of Glowinski [18], which involves solving a sequence of scalar Dirichlet problems associated with the Helmholtz operator and Neumann problems associated with the Laplacian operator. The sub-problem at second fractional step is linearized to become the classical linear di€usion-convection problem and can be solved directly. The values of a and b are chosen to yield identical Helmholtz operators in all three fractional steps, which contributes p considerably to the overall eciency of the operator splitting algorithm. A good choice for h is 1 ÿ 1= 2 for fast convergence [18±20]. The entire code is written in the Fasttalk language of Fast¯o, a ¯exible and powerful partial differential equation solving environment package developed at CSIRO [22]. 3. Free Surface Scheme Our free surface scheme mainly consists of three steps: (I) track explicitly the moving boundary and generate an intermediate mesh in the new domain; (II) ®nd a pair of orthogonal or nearly orthogonal trajectory functions in the new domain using the intermediate mesh; (III) adjust the mesh according to the orthogonal trajectories to minimize mesh distortion while maintaining the free surface shape. When the free surface movement and the corresponding mesh distortion are small, Procedure (I) is usually sucient. However, when the free surface movement is large, the accumulated mesh distortion will severely a€ect the quality of the mesh and consequently the quality of the non-linear computations if only Procedure (I) is carried out without further re®nement. Procedures (II) and (III) are developed here to prevent or reduce the deterioration of mesh qualities during free surface updating. Strictly speaking the orthogonal trajectory scheme developed here for mesh quality control are only applicable to a class of geometries where the domain can be represented by a pair of orthogonal trajectories, examples of which are many of the extrusion problems, some coating problems and the present rod-climbing problem. In many other cases where the domain does not have orthogonal boundaries, it is often possible to identify a sub-domain containing the free surface in which the orthogonal trajectory approach is applicable. 3.1. Tracking the free surface We will adopt Fast¯o's existing free surface module [22,23] to track the moving boundary, which is described below.

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

397

The solution domain is in general a time-varying domain X…t† which is to be generated by a continuous coordinate (mesh) deformation process X…t ‡ dt† ˆ X…t† ‡ dX…t†;

X…t† ! X…t ‡ dt†:

…14†

The ¯ow ®eld must satisfy the kinematic condition u
nˆ

dh dh ‡u
t dt ds

on dXf ;

…15†

where h is the normal displacement of the free surface, measured from current boundary dXf , t the unit tangential vector along the free surface. The kinematic condition relates the rate of displacement of the free surface to the current ¯uid velocity. The surface stress condition and the kinematic condition are normally suf®cient to de®ne a unique ¯ow ®eld and free boundary dXf . A kinematic iteration suitable for transient as well as steady free surface problems is carried out, at each time step we solve for a normal displacement function h…s† on the free surface satisfying the di€erential equation h dh d2 h ‡ un‡1
t ÿ a 2 ˆ un‡1
n Dt ds ds

on dXf :

…16†

With normal displacement of free surface known, the mesh displacement function dX ˆ …dX ; dY † on dX can be found accordingly and this function enables us to build a new mesh on the new domain X…t ‡ dt†. A Laplace equation is solved for each of the new mesh displacement coordinates dX …t ‡ dt† and dY …t ‡ dt† with Dirichlet boundary condition (BC) on free surface dX …t† ˆ hn
i;

dY …t† ˆ hn
j on dXf ;

…17†

where n is the unit normal vector of free surface at time t, and i and j are the unit norm of the X and Y coordinates, respectively. On other boundaries (non-moving) either natural or zero Dirichlet BC is imposed, depending on whether slip of a particular nodal coordinate is allowed. For example in the rodclimbing problem on the bottom of the container both dX and dY will have zero Dirichlet BC, but a natural BC for dY is imposed on the wall of the rod to allow nodal points to slip along the rod, although no-slip BC is still imposed for the momentum equation. The new mesh obtained will be treated only as an intermediate mesh on which further adjustment needs to be done to reduce mesh distortion. 3.2. Finding the orthogonal trajectory functions As shown in Fig. 1, we assume the computational domain X…t† is bounded by curves AB, DC, AD and BC, and they are orthogonal to each other at the corners (our present case of rod-climbing satisfys this restriction). Although many computational domains will not satisfy the above geometric restriction, often part of the domain which contains the free surface does satisfy the restriction, i.e., it is often possible to identify a local free surface region which are bounded arti®cially by orthogonal boundaries. We further assume the curves are suciently smooth so that a pair of orthogonal or nearly orthogonal functions W…t† and U…t† exist in the domain which satisfy rW
rU ˆ 0 (or rW
rU  0) everywhere. We seek such a pair of functions by solving Laplace equations with the following boundary conditions, 4 W ˆ 0;

W ˆ 0 on AB;

W ˆ 1 on DC;

4 U ˆ 0;

U ˆ 0 on AD;

U ˆ 1 on BC;

rW
n ˆ 0 on AD and BC; rU
n ˆ 0 on AB and DC:

…18† …19†

It is clear that if W and U are such a pair of orthogonal functions, then W0 ˆ aW ‡ b and U0 ˆ cU ‡ d are also orthogonal to each other, with a, b, c and d being arbitrary constants, therefore taking the values of zero and one as Dirichlet BC for W and U does not lose generality. At each time step after the new free surface is located by means of Procedure (I) as described in Section 3.1, W and U are solved in the new domain X…t ‡ dt† on the intermediate mesh, and they are used to modify the mesh to minimize mesh distortion, as discussed in the following section.

398

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Fig. 1. An illustration of orthogonal trajectory mesh.

3.3. Adjusting mesh according to the orthogonal trajectories Assuming a structured quadrilateral mesh is initially built at time t on orthogonal trajectories as illustrated in Fig. 1, so all the element lines are either on W ˆ const or U ˆ const lines. As shown in Fig. 1 the element E with corner points Q and P are on four trajectory lines passing the corner points: U ˆ Up , U ˆ Uq , W ˆ Wq and W ˆ Wq . In general the intermediate mesh generated through Procedure (I) does not automatically satisfy the orthogonality condition because the domain has changed its shape and nodal points have been moved. The task now is to adjust the position of each nodal point so that its W and U values, as solved in the new domain, recover the original constants and therefore the orthogonality of the mesh is restored. Referring to Fig. 1, the nodal point P has initial trajectory values of Wp and Up at location Xp ˆ(Xp , Yp ). For a convenient discussion we denote the nodal point P in any one of the three ways: P …Up ; Wp †, P …Xp ; Yp † or P …Xp †. At time t ‡ dt on the intermediate mesh the solution of W and U is denoted by W0 and U0 . In the intermediate mesh the nodal point P becomes P 0 which is denoted by P 0 …W0p ; U0p † or P 0 …Xp0 ; Yp0 † or P 0 …X0p †. One could use a searching process similar to particle tracking to ®nd the position for P where its orthogonal function coordinates (W; U) recover its original values. Instead we propose here a more robust and ef®cient algorithm, making use of the gradient ®elds of W0 and U0 and the fact that the intermediate location is not far from the target location. The assumption of the intermediate location being not far from the target location is always valid one, since the overall non-linear iterations of the viscoelastic ¯ow system restrict the time step magnitude and the free surface movement at any one time step. Denoting the target location as X00p ˆ …Xp00 , Yp00 ), Taylor series expansion gives Wp ÿ W0p  …rW†0p
…X00p ÿ X0p †;

…20†

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Up ÿ U0p  …rU†0p
…X00p ÿ X0p †; 0

399

…21† 0

where the gradients …rW†p and …rU†p are evaluated at location X0p . From the two Eqs. (20) and (21) we can solve for the two unknown components Xp00 and Yp00 . Unfortunately Eqns. (20) and (21) are only approximations, the new values of W00 and U00 at the new location X00p are still not equal to Wp and Up , but closer to them than W0p and U0p were. Therefore it is necessary to carry out an iterative process until the ®nal target position is located where W00p ˆ Wp and U00p ˆ Up . One diculty of the iterative process is that at X00p the values W00 and U00 are not readily known, since X00p is not a nodal point of the intermediate mesh any more after being moved from X0p . In fact X00p may not necessarily be within the same mesh. Again we chose to avoid a searching process to ®nd (W00p ; U00p ), instead we form a new intermediate mesh based on new locations of X00 and solve for W and U ®elds on the new intermediate mesh (but in the same domain with the same free surface). Using the new intermediate mesh and the new solutions of W and U the above procedure of adjustment can be repeated until a preset tolerance of Wp ÿ W00p and Up ÿ U00p is met. Fortunately the above iterative process typically takes only two or three iterations before the norm of Wp ÿ W00p and Up ÿ U00p (for all the nodes) goes to below 10ÿ6 . We have also found it is essential to use double precision for the mesh adjustment process, and this is especially true in the case of ®ne mesh (or high local mesh concentration) where the adjustment is only a very small fraction of the tiny mesh size. Although the above procedures appear to be computer-time consuming, it only takes about 15% of total CPU time for solving the entire viscoelastic problem of rod-climbing. For unstructured triangular meshes the concept of orthogonal grid does not naturally apply, but the idea of using orthogonal trajectories to minimize mesh distortion is still applicable, if we keep the W and U values of each nodal point the same as the starting values of the original mesh at all time steps, the shape of each triangle is nearly preserved, thus the distortion of the mesh is reduced, as will be illustrated below by an example. 3.4. Example 1: Structured quadrilateral mesh Fig. 2(a) shows an initial quadrilateral mesh on a rectangular domain. Assuming the top boundary is a free surface and it is subjected to an arti®cial normal displacement described by a sine function at every time step, we proceed to update the mesh using the above orthogonal trajectory scheme. Note the sine function is imposed on the normal displacement, not on the free surface shape itself, so the shape of the free surface is not necessarily sine wave after a few time steps. Fig. 2(b) shows the mesh updated using our orthogonal trajectory scheme after 5 time steps and Fig. 2(c) shows an enlarged view of a section near the free surface. It is clear the orthogonality of the mesh is preserved despite the large deformation of the domain. As a comparison Fig. 3(a) shows the mesh updated only by Procedure (I) as described in Section 3.1, i.e., without orthogonal trajectory adjustment. Obviously the mesh updated without the orthogonal trajectory adjustment has lost orthogonality and it is very much distorted, especially near the free surface, as the enlarged view in Fig. 3(b) illustrates more clearly. In this example the movement of free surface at each time step is unrealistically large. In a real viscoelastic ¯ow simulation the free surface movement at each time step is made very small by limiting the magnitude of the time step, since the non-linear iteration requires small change of domain in each time step to ensure stability and convergence. 3.5. Example 2: Unstructured triangular mesh Fig. 4(a) is the original triangular mesh on a rectangular domain. This time we assume both the top and the bottom boundaries are free surfaces subjecting to continuous normal displacement described by a sine function.

400

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Fig. 2. Structured quadrilateral mesh: (a) initial mesh; (b) mesh updated using orthogonal trajectory scheme; (c) enlarged view of the updated mesh near the centre of the free surface.

Fig. 4(b) shows the mesh updated using our orthogonal trajectory scheme after 8 time steps, with Fig. 4(c) showing an enlarged view of a free surface section near the top-right corner. One can see the mesh distortion is very small and the overall quality of the mesh after the large deformation is still good, compatible with the original mesh. As a comparison, Fig. 5(a) shows the mesh updated without the procedure of orthogonal trajectory adjustment, and an enlarged view of a part of the free surface (the same section as in Fig. 4(c)) is shown in Fig.

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

401

Fig. 3. Structured quadrilateral mesh: (a) mesh updated without using orthogonal trajectory scheme; (b) enlarged view of the updated mesh near the centre of the free surface.

5(b). Obviously the mesh is distorted signi®cantly, especially near the free surface, as the enlarged view in Fig. 5(b) illustrates more clearly. Some of the angles in the triangle mesh are either less than 20 or larger than 160 . The di€erence in mesh quality shown in the comparison of Figs. 4(c) and 5(b) convincingly demonstrates the advantages of using our orthogonal trajectory scheme, this di€erence often means a di€erence between good and poor convergence or even between convergence and divergence in viscoelastic ¯ow simulations involving large deformation of the domain due to free surface motion.

4. Numerical results The ¯ow description and boundary conditions are shown in Fig. 6. A rod of radius R rotates at a constant angular speed x about the axis of a stationary cylindrical vessel containing a viscoelastic ¯uid. In the secondary ¯ow in the r±z plane the viscoelastic stress will force the ¯uid to ¯ow inward against centrifugal force and climb upwards against gravity.

402

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Fig. 4. Unstructured triangular mesh: (a) initial mesh; (b) mesh updated using orthogonal trajectory scheme; (c) enlarged view of the updated mesh near the top-right corner of free surface.

Four-dimensionless numbers are needed to characterize the rod-climbing ¯ow, the Weissenberg number Wi, Reynolds number Re, Stokes number St and the capillary number Ca, which are de®ned as follows: Wi ˆ kx;

Re ˆ

qxR2 ; g

St ˆ

gx ; qgR

Ca ˆ

gxR ; C

where xR is taken as the characteristic velocity in the ¯ow.

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

403

Fig. 5. Unstructured triangular mesh: (a) mesh updated without using orthogonal trajectory scheme; (b) enlarged view of the updated mesh near the top-right corner of free surface.

Our numerical study mainly consists of two parts: ®rstly the Maxwell model is used to solve the rodclimbing problem under the conditions of creeping ¯ow (Re ˆ 0), zero surface tension (Ca ˆ 1) and low angular speed, and the results are compared with perturbation theory which is valid only for small angular speed. In the second part of our numerical study the PTT model is used to simulate some real experiments of rod-climbing ¯ow by Beavers and Joseph [7], taking into consideration all the important factors a€ecting the ¯ow: surface tension, gravity, centrifugal force and second normal stress di€erence. For the ®rst part only a triangular mesh was used, but for the second part a quadrilateral mesh as well as a triangular mesh was used to check the mesh-independence of the numerical solutions. 4.1. Maxwell model results in comparison with theory At a low angular velocity and ignoring surface tension and inertia, the second-order theory (Joseph et al. [7,9]) gives the free surface shape for the UCM ¯uid as

404

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Fig. 6. Geometry, ¯ow description and boundary conditions for the rod-climbing problem.

z…r† 2Wi2 ˆ ; R D…r=R†4

…22†

where D is the ratio of the Wi to the St. This expression gives a 1=r4 curve for free surface shape and a linear relation between climb height and the square of Wi. Fig. 7 shows the starting mesh and geometry. The depth of the ¯uid is 6R and the radius of the vessel is 10R and the bottom of the rod is one R distance above the bottom of the vessel. Assuming D ˆ 0:01 (same order of magnitude as in real experiments), three small angular velocity cases of Wi ˆ 0:01; 0:015; 0:02 were computed. Since the climbing height in these cases is expected to be small, the resulting mesh deformation is also small, so it is suf®cient to just use Procedure (I) as described in Section 3.1 without using the elaborate mesh adjustment Procedures (II) and (III). Fig. 8 compares the numerical and analytical predictions on the free surface shape, which shows the two agree very well, perhaps except near the contact point (r=R ˆ 1) on the rod-climbing. Indeed at low Wi the free surface shows a 1=r4 dependence, except near the contact point where the numerical simulation predicts a ¯at angle (90 ) while analytical one gives a sharper angle. The di€erence in contact angle between numerical and perturbation analysis is an interesting one and worth some comments, experimental observations [7] reveal that if there is no static wetting due to surface tension, i.e., no initial climbing when the ¯uid is stationary (x ˆ 0), the contact angle at the rod remains ¯at (90 ) after the climbing and the free surface shape near the rod is convex. Since both numerical and analytical results in Fig. 4 are based on zero static climbing (due to zero surface tension), the above experimental observation means our numerical prediction is more realistic than the

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

405

Fig. 7. Geometry and starting mesh for the rod-climbing ¯ow of UCM ¯uid.

Fig. 8. Predictions of the free surface shapes. N, this work using UCM; ± , perturbation theory [7,9].

perturbation analysis. From computational point of view, the absence of surface tension means the contact angle is not a boundary condition, but it is part of the solution. In other words, one cannot specify the contact angle in the absence of surface tension, otherwise the problem would be incorrectly over-speci®ed. At higher Wi values or higher angular rod speeds the perturbation theory over-predicts the climbing height as compared with numerical predictions. As shown in Fig. 8 the difference between numerical and perturbation theory predictions in climbing height at W i ˆ 0:02 is already more than what can be expected from the difference in contact angle. The higher the angular velocity, the larger the difference is. For example, at W i ˆ 0:05 our numerical prediction on climbing height is h=R ˆ 0:32, while the perturbation theory (Eq. 22))) gives h=R ˆ 0:5 which is about 50% higher than the numerical prediction. The numerical work of Debbaut and Hocq [14] also con®rmed this trend.

406

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

4.2. PTT model results in comparison with experiments To realistically simulate the experiments of Beavers and Joseph [7], other factors such as surface tension, centrifugal force and second normal stress di€erence have to be included. The PTT model was used to cope with the non-zero second normal stress di€erence. In the experiments reported by Beavers and Joseph [7], STP (a 26.6% polyisbutylene solution in petroleum oil) was used, for which the following material constants were given at temperature 25:5 C, q ˆ 0:89 g cmÿ3 ;

g ˆ gs ‡ gm ˆ 146 poise; ÿ2

C ˆ 30:9 dyn cm s ;

  gs k ˆ 0:0162 s; 1ÿ gm

where gs is the newtonian solvent viscosity. The cylindrical vessel has a radius of RV ˆ 15:25 cm, the height of the ¯uid at rest is 7.7 cm. A few experiments were carried out with three rods (R ˆ 0:635; 0:476; 0:317) rotating at di€erent speeds, and we have selected four cases to simulate: three cases with R ˆ 0:635 cm and one case with R ˆ 0:476. The operating parameters and dimensionless numbers for the four cases are listed in Table 1. The case with R ˆ 0:476 and x ˆ 5 rev: sÿ1 had one of the highest rotating speeds and climbing heights among all the steady state experiments done with STP and coated rods at room temperature [7]. The experiments found [7] that at higher rotating speeds the steady drop-like con®guration loses its stability to a time-periodic motion in which the climbing ¯uid may collapse downwards to the main body of ¯uid. The present study only deals with steady state solutions. Four material parameters in the PTT model need to be determined or estimated before starting the simulation: gm , k, n and . In the absence of direct experimental data, we have let  ˆ 0 and taken the commonly used ratio gs =gm ˆ 1=9 to determine gm from g ˆ gs ‡ gm ˆ 146 poise and the relaxation time k from …1 ÿ …gs =gm ††k ˆ 0:016 s. The second normal stress parameter n can be estimated from the climbing constant given by the perturbation theory,    g …23† b^ ˆ 1 ÿ s kg …1 ÿ 2n†; gm where b^ is the climbing constant which can be measured by carrying out low angular speed experiments [7]. For the case of STP at temperature 26 C the experimental value of b^ was found to be b^ ˆ 1:02 [7], from which and Eq. (23) we get an estimate of n ˆ 0:284. Experiments [7] also revealed that the climbing constant is quite sensitive to temperature. For example, at 27:9 C the value of b^ was found to be 0:86, comparing with b^ ˆ 1:02 at 26 C. From Eq. (23) a b^ value of 0:86 gives n ˆ 0:318. We have done calculations with both n ˆ 0:284 and n ˆ 0:318. The corresponding ratios of second to ®rst normal stress di€erence given by PTT model are (N 2=N 1 ˆ ÿn=2) ÿ0:142 and ÿ0:159, respectively, for the two n values. In order to study the rod-climbing due to viscoelastic force without the e€ect of the static climbing, the rods in the experiments were coated with Scotchgard, so that when the ¯uid was stationary there was no static climbing due to surface tension, and the contact angle of ¯uid with the rod was ¯at (90 ) at all rotational speeds. In this way the static climbing was eliminated but at the same time the dynamic e€ect of surface tension on rod-climbing was maintained. Table 1 The operating parameters and dimensionless numbers for the four experiments simulated in this work R (cm)

x …rev: sÿ1 )

Re

Wi

St

Ca

0.635 0.635 0.635 0.476

1.70 2.10 2.60 5.00

0.0263 0.0325 0.0402 0.0435

0.1944 0.2401 0.2973 0.5718

2.813 3.475 4.302 11.037

32.05 39.59 49.02 70.66

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

407

The two starting ®nite element meshes are shown in Fig. 9, the height (H0 ) and the width (RV ÿ R) of the domain are 12:126R and 23:016R, respectively, the same as in the experiment. The triangular mesh has a total of 3873 nodes which gives a total of 36731 degrees of freedom (DOF) for the rod-climbing problem, and the quadrilateral mesh has a total of 3516 nodes with a DOF of 32769. Both meshes have a high concentration of nodes near the contact point of free surface on the rod. Because relatively large free surface movements occured in the experiments (an order of magnitude larger than the low angular speed cases considered in the previous section), the full orthogonal trajectory free surface scheme was employed to reduce mesh distortion and therefore to obtain better convergence at higher Wi values. The convergence criterion for steady state was set as H

j un‡1 ÿ un j < 10ÿ4 ; j un‡1 j dt

i.e., the pre-set tolerance has to be satis®ed by the norm (computed over the entire domain) of the time derivative of the relative velocity di€erence between two consecutive time steps. It was found this convergence criterion would make the relative di€erence in free surface height between two time steps to be less than 0.02%, i.e.,

Fig. 9. Geometry and starting meshes for the simulation of experiments.

408

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

Fig. 10. Predicted free surface pro®les in comparison with experimental measurements for the three cases with R ˆ 0.635 at different rod speeds. Ð, this work using triangular mesh with n ˆ 0:318; ÿ ÿ ÿ, this work using triangular mesh with n ˆ 0:284; ‡, this work using quadrilateral mesh with n ˆ 0:318 (for x ˆ 2:60 rev: sÿ1 only); }, measurements by Beavers and Joseph [7]; , Debbaut and Hocq [14] with n ˆ 0:275 (for x ˆ 1:70 rev: sÿ1 only).

hn‡1 ÿ hn < 2  10ÿ4 ; hn where h is de®ned as the rise above the initial stationary level H0 , as illustrated in Fig. 6. We found typically about 400 time steps with dimensionless time step dt ˆ 4 are required to reach steady state convergence. We

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

409

also checked the total mass conservation and found the relative error to be less than 10ÿ5 after a few hundred time steps, so no attempt was necessary to explicitly balance the total mass in any of the computations. Fig. 10 compares the predicted free surface pro®les with experimental measurements for the three cases with R ˆ 0:635 at di€erent rod speeds, showing computational results of both n ˆ 0:284 and n ˆ 0:318. Prediction of Debbaut and Hocq [14] for rod speed of 1.70 rev. sÿ1 is also shown. Free surface curves predicted for higher rod speeds are not available from Ref. [14]. There is little di€erence in predicted results between triangular mesh and quadrilateral mesh, as shown in Fig. 10 for the case of x ˆ 2:60 rev: sÿ1 , the case with the largest rotation speed in Fig. 10. The agreement between predictions and experiments is in general good, given that the maximum absolute di€erence (with dimension) in climbing height between our prediction with n ˆ 0:284 and the measurement is less than 0.05 cm, i.e., less than half a millimeter, and the maximum di€erence between our prediction with n ˆ 0:318 and the measurement is less than 0.02 cm. Calculations with n ˆ 0:284 over-predicted the climbing height at all rod speeds, while results with n ˆ 0:318 are closer to experimental measurements. The present numerical results con®rmed that in general the second normal stress reduces rod-climbing. It is worth emphasising here that the climbing constant b^ is a relation derived from perturbation theory valid only for small angular velocities, so the value of n ˆ 0:284 (or n ˆ 0:318 at slightly di€erent temperature) computed from the climbing constant is only an estimation. A more valid comparison of numerical predictions with experimental measurements can only be done with more reliable data from direct measurements of second normal stress di€erence. Fig. 11 shows the ®nal quadrilateral mesh, the ®nal triangular mesh and the velocity vectors (in the r±z plane) in the region near the contact point for the case of x ˆ 2:60 …rev: sÿ1 † with n ˆ 0:318. One can see the quality of the triangular mesh is still good after the relatively large deformation of the domain due to free surface movement, and the orthogonality of the quadrilateral mesh is preserved. The corner

Fig. 11. Final quadrilateral mesh (a) triangular mesh (b) and the velocity vectors (c) in the region near the contact point for the case of x ˆ 2:60 …rev: sÿ1 † with n ˆ 0:318.

410

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

recirculation is clockwise while the main recirculation is counterclockwise. The velocity arrows on the free surface show clearly that the kinematic condition of Eq. (15) is well satis®ed by the numerical solution. Fig. 12 plots contours of the ®rst normal stress N 1 ˆ shh ÿ srr and the second normal stress N 2 ˆ srr ÿ szz for the same case, taking r±h plane as the principle shearing plane. Both normal stresses take their absolute maximum values on the surface of the rod towards its lower part. Note the maximum N 1 value of 3:8876 and the maximum (absolute) N 2 value of ÿ0:61813 gives a ratio of N 2=N 1 ˆ ÿ0:159, consistent with PTT model prediction of N 2=N 1 ˆ ÿn=2 for simple shear ¯ow. The positive ®rst normal stress indeed forms `hoops' around the rod, it creates a positive tension in h direction, forcing the ¯uid inwards against the centrifugal force and upwards against the gravitational force. The negative second normal stress acts vertically as a tension and since it is zero on the free surface and it has higher (absolute) values towards the lower part of the rod, it pulls the ¯uid downwards which cancels out some of the ®rst normal stress e€ect and reduces the climbing height. Contours of the azimuthal velocity component uh and pressure p are shown in Fig. 13.

Fig. 12. Contours of (a) ®rst normal stress N 1 ˆ shh ÿ srr and (b) second normal stress N 2 ˆ srr ÿ szz for the case of x ˆ 2:60 …rev: sÿ1 † with n ˆ 0:318. Nine countour levels are uniformly distributed between the minimum and the maximum values.

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

411

Fig. 13. Contours of (a) azimuthal velocity uh and (b) pressure p. Relative contour levels for uh are evenly distributed between the minimum (0) and the maximum (…uh †max ˆ 2:60  2pR), with relative increment of duh =…uh †max ˆ 0:05; Relative contour levels for pressure are evenly distributed from zero (on the free surface) to maximum with relative increment of dp=pmax ˆ 0:02.

Fig. 14. Three-dimensional free surface shapes for the case of R ˆ 0:476 cm, x ˆ 5 …rev: sÿ1 †. Left: this work with n ˆ 0:318; right: photo from Ref. [7].

For the case of R ˆ 0:476 and x ˆ 5 rev: sÿ1 , no free surface curve was given in ref. [7], but a photograph was presented instead. In order to compare our numerical prediction with the experiment, we have created three-dimensional graphics from the computed free surface pro®le. Fig. 14 compares the numerically predicted three-dimensional free surface with the photograph presented in ref. [7]. The similarity of the two is satisfactory, apart from the apparent re¯ection in the lower half of the photo and some other optical e€ects in the photo. According to ref. [7], the apparent wavelike shape of the free surface in the photo was just an optical e€ect due to a light re¯ection from the rod. Both the numerical graphics and the experimental photo show smooth bell-shaped free surface with similar heights. The predicted climbing height in this case is h=R ˆ 1:093, while the estimated climbing height in the experiment by measuring the photo is about h=R ˆ 1:1, although accurate estimation is dicult due to the three-dimensional aspect and the optical e€ects in the photo. The agreement on the width of the bell-shaped free surface between prediction and experiment appears to be not as good as that of the height. The bell-shape predicted is noticeably narrower than that of the experiment shown in the photo.

5. Conclusions Full ®nite element simulations of steady state rod-climbing ¯ow have been performed using UCM and PTT models. The predictions are in good agreement with second-order theory and experimental data. The orthogonal trajectory scheme developed in this work is applicable to an unstructured triangular mesh as

412

X.-L. Luo / Comput. Methods Appl. Mech. Engrg. 180 (1999) 393±412

well as a structured quadrilateral mesh, and it provides ¯exibility and good quality in mesh updating for successful numerical simulation of free surface ¯ow involving relatively large deformation in the domain. Without using the orthogonal trajectory scheme the simulation failed to converge even for very moderate climbing height. Although the current numerical scheme is already a time stepping one, time accuracy achieved here is only ®rst-order because of the way the mesh and solution ®elds are updated at each time step. A very challenging further work is to develop a more elaborate transient scheme for viscoelastic free surface ¯ow which is second-order accurate in time and use it to simulate the periodic behaviour of rod-climbing ¯ow at higher angular speeds. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

K. Weissenberg, Nature 159 (1947) 310±311. R.B. Bird, R.C. Armstrong O. Hassager, Dynamics of Polymeric Liquids, 2nd ed., vol. 1, Wiley, New York, 1987. R.I. Tanner, Engineering Rheology, Clarendon Press, Oxford, 1985. H.A. Barnes, J.F. Hutton, K. Walters, An introduction to Rheology, Elsevier, Amsterdam, 1989. K. Walters, Rheometry, Chapman and Hall, London, 1975. H. Giesekus, Rheol. Acta 1 (1961) 404±413. G.S. Beavers, D.D. Joseph, The rotating rod viscometer, J. Fluid Mech. 69 (1975) 475±511. D.D. Joseph, G.S. Beavers, Free surface induced by the motion of viscoelastic ¯uids, in: R.S. Rivlin (Ed.), The Mechanics of Viscoelastic Liquids, AMD 22, ASME. Noe York, 1977, pp. 59±99. Y. Yoo, D.D. Joseph, G.S. Beavers, Higher-order theory of the Weissenberg e€ect, J. Fluid Mech. 92 (1979) 529±590. Beavers, Y. Yoo, D.D. Joseph, The free surface on a liquid between cylinders rotating at di€erent speeds. Part III, Rheol. Acta 19 (1980) 19±31. D.D. Joseph, Fluid Dynamics of Polymeric Liquids, Springer, New York, 1990. D.D. Joseph, G.S. Beavers, A. Cers, C. Dewald, A. Hoger, P.T. Phan, Climbing constants for various liquids, J. Rheology 28 (1984) 325±345. M.J. Crochet, A.R. Davies, K. Walters, Numer. Simulation of Non-Newtonian Flow, Elsevier, New York, 1984. B. Debbaut, B. Hocq, J. Non-Newtonian Fluid Mech. 43 (1992) 103±126. J.M. Marchal, M.J. Crochet, Non-Newtonian Fluid Mech. 26 (1987) 77±114. S.F. Kistler, L.E. Scriven, Coating ¯ows, in: J.R.A. Pearson, S.M. Richardson (Eds.), Computational Analysis of Polymer Processing, Applied Sciences Publishers, London, 1983, pp. 243±299. X.L. Luo, J. Non-Newtonian Fluid Mech. 63 (1996) 121±140. R. Glowinski, Finite element methods for the Numerical Simulation of Incompressible Viscous Flow. in: Lectures in Applied Maths, vol. 28 1991. M.O. Bristeau, R. Glowinski, J. Periaux, Comput. Phy. Rep. 6 (1987) 73±187. E. Dean, R. Glowinski, C.H. Li, in: Goldstein, J. Rosecrans (Eds.), Mathematics Applied to Science, S. & Sod, G. 1988. A.N. Brooks, T.J.R. Hughes, Comput. Meth. Appl. Mech. Eng. 32 (1982) 199. Fast¯o Tutorial Guide (1997), CSIRO Mathematical and Information Sciences. More information is available from WWW site:http://www.nag.co.uk/simulation/Fast¯o/fast¯o.html. J.R. Mooney, F.R. DeHoog, Modelling steady state free surface ¯ows using Fast¯o, in: Proceedings of the 12th Australasian Fluid Mechanics Conference, University of Sydney, 1995, pp. 407±410.