- Email: [email protected]
On the numerical simulation of axisymmetric swirling flows of differential viscoelastic liquids: the rod climbing effect and the Quelleffekt
Journal of Non-Newtonian Fluid Mechanics, 43 (1992) 103-126
103
Elsevier Science Publishers B.V., Amsterdam
On the numerical simulation of axisymmetric swirling flows of differential viscoelastic liquids: the rod climbing effect and the Quelleffekt Benoit Debbaut 1,2 and Bernard Hocq 2 ’ Unite’de Mtkanique Appliquke, Universite’ Catholique de Louvain, Place du Levant 2, B-1348 Louvain-la-Neuve (Belgium) 2 Polyflow S.A., Place de I’Universite 16, B-1348 Louvain-la-Neuve (Belgium) (Received May 23, 1991; in revised form November 22, 1991)
Abstract
A finite-element investigation is presented of some Weissenberg effects occurring in swirling flows of viscoelastic liquids in geometries of revolution. Numerical results are shown for the rod-climbing flow and the so-called Quelleffekt of viscoelastic liquids of the differential type. For both types of flow, comparisons with earlier numerical or analytical predictions and some experimental measurements are also shown. Keywords: rod climbing; Weissenberg
effect; viscoelastic liquids; swirling flows; numerical
simulation.
1. Introduction
Early experimental observations by Garner and Nissan [l] with a solution of rubber in benzene showed a rod-climbing effect. This led Weissenberg to introduce the hypothesis of pull, or normal stress, acting along streamlines [2]. Such a pull appears to be responsible for the so-called Weissenberg effect (or rod climbing). Indeed, as a physical explanation of this effect, it was suggested in Ref. 3 that the normal stress “acts like a hoop stress around the rod” and in Ref. 4 “strangulates the fluid, and forces it Correspondence to: B. Debbaut, Unite de Mecanique AppliquCe, Universite Louvain, Place du Levant 2, B-1348 Louvain-la-Neuve, Belgium.
0377-0257/92/$05.00
Catholique
0 1992 - Elsevier Science Publishers B.V. All rights reserved
de
104
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103-126
inwards against the centrifugal force and upwards against the gravitational force”. Such a flow belongs to the class of axisymmetric swirling flows which, from a kinematic viewpoint, are characterized by three velocity components depending upon two space variables (r, z) in a cylindrical reference frame. The rheometrical interest of this class of flows has already been discussed [3-51. In the literature, several papers have been devoted to analytical developments using perturbation theory, for rod climbing [6-121, indicating a relationship between rod-climbing measurements and the viscoelastic properties of the liquid. On the basis of analytical developments for Couette flow with a free surface, Serrin [6] and Giesekus [7] discussed the origin of the fluid climbing on the inner cylinder. In Ref. 7, Giesekus also introduced the concept of positive and negative Weissenberg effects. In Refs. 8-12, the authors showed that a knowledge of the climbing was sufficient for determining an important viscoelastic parameter of a simple fluid. Experimental measurements carried out by Beavers and Joseph [S] on STP (polyisobutylene in petroleum oil) showed good quantitative agreement with analytical predictions. Joseph et al. [13] have tabulated experimental rod-climbing measurements for several fluids. Yoo et al. [lo] also analysed the development of secondary motions as a function of the rod angular velocity. Hu et al. [14] and Joseph [15] have discussed the control of secondary motions by the extensional viscosity, which was already known to be important in vortex development in contraction flows [16]. Among experimental works for various fluids in a disk and cylinder system, Hill [17] noted a reversal in the direction of the secondary motion as the fluid becomes “more conductive to viscoelastic secondary flows”. A phenomenon related to the Weissenberg effect is the Quelleffekt investigated by Bijhme et al. [18]. Here, the viscoelastic liquid is contained in a cylindrical vessel. The fluid is set into motion by a rotating disk located at the bottom of the vessel. It flows upwards along the axis of symmetry, as a source (QueU), and produces a bulge in the free surface around the axis. This effect can be understood as a “rod-climbing without a rod” [4]. In Ref. 18, the authors analysed the effects of inertia and normal stress differences upon bulge shape. Also, the predicted results showed good agreement with experimental measurements for both silicone oil and a 2.5% polyacrylamide (PAA) aqueous solution. PAA, which is much more elastic than silicone oil, is also characterized by a non-constant shear viscosity with a small plateau and by a non-constant normal stress coefficient with a larger plateau. However, the authors in Ref. 18 displayed results obtained in the range of constant normal stress coefficient and
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
105
showed that the bulge shape is not affected by a non-constant viscosity. An important feature, seen both theoretically and experimentally, is that climbing and axial bulge deformation are quadratic in the angular velocity of the rotating device (rod or disk), for low angular velocity [g-13,18]. In view of the air-fluid interface in the experiments, surface tension effects have also been investigated [8,12,18,19]. In particular, in experiments with STP, Beavers and Joseph [8] removed the static contact angle by coating the rod with Scotchgard; with coated and uncoated rods, the corresponding theoretical predictions still matched the experimental measurements. Delgado et al. [19] have analysed the limiting case as the gravitational force vanishes. In their experiments with PAA, Bijhme et al. [18] found that the influence of surface tension is of the order of a few percent. Finally, we should mention some studies on the stability of rod-climbing flow. As the angular velocity increases, the flow can show a time-periodic behaviour, the breathing instability, which is discussed in Refs. 8 and 20. Our purpose in the present paper is to extend the capabilities of numerical algorithms widely used for 2-D flows of differential viscoelastic liquids [21-251 to axisymmetric swirling flows. Numerical simulations of such flows have already been done for viscoelastic fluids of integral type [26]. Although a large variety of axisymmetric swirling flows of viscoelastic liquids are discussed in Ref. 12, we shall investigate here the finite-element simulation of two flows with a free surface: the rod climbing and the Quelleffekt of both Oldroyd-B and Johnson-Segalman liquids [27,28]. Such flow problems have been selected, not only in view of their rheometrical interest, but also because analytical and experimental results already exist [8,18]. Furthermore, exp”erimenta1 measurements were obtained in a range of shear rates for which normal stress coefficients are constant. In particular, we shall focus on the free surface shape and on the development of secondary motions. The present work follows an early finite-element investigation done by Georges and Hocq [29] on the rod-climbing effect for both Newtonian and second-order fluids. The equations and the finite-element discretization are discussed in Section 2. In Section 3, we analyse the rod-climbing flow of an upper-convected Maxwell liquid and identify some interesting features. Sections 4 and 5 are devoted to rod climbing and Quelleffekt of Oldroyd-B [27] and Johnson-Segalman [28] liquids; comparisons with earlier predictions and experimental results are discussed. 2. Equations
and discretization
Let us denote the rate of deformation d = $7~
+ bT),
tensor d as follows:
(1)
106
B. Debbaut and B. Hocq/J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
where v is the velocity field. For the models considered stress tensor T can be decomposed as follows: T= Tl + T2,
here, the extra(2)
where Tl and T2 are, respectively, the viscoelastic and the Newtonian components to the extra-stress tensor T. These components obey the following constitutive relationships:
where q is the total shear viscosity, A is the relaxation time, and the symbols v and A denote the upper-convected and 1ower;convected derivative operators, respectively. In (3), 5 is a material parameter controlling the ratio of the second to the first normal stress differences, while C#Jis the ratio of the Newtonian solvent viscosity to that of the solution. As 5 vanishes, one recovers the Oldroyd-B model [27]; furthermore, as 4 and 5 vanish, one recovers the upper-convected Maxwell liquid. The constitutive relationships (3) are coupled with the momentum and incompressibility equations, i.e.: -Vp+VT+pg=pvVv,
(4)
v*v=o,
(5) where p is the pressure field, p is the fluid density and g is the gravity vector. For the class of axisymmetric swirling flows that we consider here, the unknown fields T,, v and p depend only upon two space variables in a cylindrical system of axes (e,, eZ), and do not vary in the azimuthal direction e,. However, we must point out here that the fields Tl and v now have 6 and 3 components respectively, instead of 4 and 2 as in the case for 2-D axisymmetric flows [21]. In terms of numerical simulation, which will be discussed below, swirling flow problems can be handled with the same computational algorithms that were used for 2-D flow situations. The computational cost increases due to the additional amount of components for Tl and v, but remains much lower than that required for a full 3-D simulation [30]. In addition to eqns. (l)-(5), additional relationships must be satisfied. Indeed our flow problems involve free surfaces, the location of which is a priori unknown. The kinematic condition on free surfaces can be written as follows: v*n=o,
(6)
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (192)
103- 126
107
where it is the normal vector to the free surface. In some simulations, surface tension a, will he taken into account. This obeys the relationship: 1 -+’
q=r i
R1
I n
R,
’
(7)
where I’ is the surface tension coefficient, while R, and R, are the principal radii of curvature of the free surface. Our flow problems occur in closed volumes, so that any solution must satisfy the global form of volume conservation: /V
dV-
V, = 0,
where V,, is the initial volume at rest. 3 characterize the flow, we use non-dimensional numbers. Therefore, in addition to the given material and physical data q, A, p, I and g, we select a characteristic dimension R and angular velocity R of the rotating device. We define the Reynolds, Weissenberg, Stokes and capillary numbers, respectively, as follows: pflR2 Re = rl ’ We=hlR, 770 St = PgR ’ 77fiR Ca = r *
Although Re is generally very small, acceleration terms in (4) will be taken into account in some simulations. It has been shown [7,18,26] that inertia forces generate a secondary motion opposed to that of viscoelasticity. Moreover, although the Stokes number St may be large, gravity forces cannot be neglected. Indeed, both rod climbing and the Quelleffekt result from viscoelastic normal stresses acting against gravitational forces. A further dimensionless number, 9, can be obtained as the ratio of the Weissenberg number to the Stokes number: $+_
QgR 77 *
(10)
This number is independent of the angular velocity CI, so that, for a given fluid with its material properties A, p and q, in a given environment specified by g, 9 characterizes the flow geometry.
108
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) IO3- 126
The numerical simulation is performed by means of two finite-element algorithms. (i) The algorithm developed by Marchal and Crochet [22], which combines 4 X 4 sub-element interpolations for the viscoelastic extrastress tensor Z’i with a streaml@e upwgding (SU) technique applied to the convective term appearing in ?‘r and T1. Mesh-convergence properties for the velocity field and streamlines were already proved in Refs. 22 and 23. (ii> The earlier MIX1 algorithm [21] with biquadratic interpolation for the viscoelastic extra-stress is also used for several reasons. As will be shown, interesting solutions are obtained at relatively low values of We. Moreover, our flow problems are smooth enough in the regions of interest, so that the technique MIX1 can be successfully used, as suggested in Ref. 31. This technique is also cheaper from a computational viewpoint for refined meshes. The use of two different algorithms will also allow us to assess the convergence of the solutions, although solutions at higher We will be obtained only with the 4 X 4 SU technique. In this paper we do not intend to track possible bifurcation points; thus, calculations will be pursued as long as interesting results can be obtained. For both techniques, biquadratic and bilinear interpolants are used for the velocity and pressure fields, respectively. The finite-element mesh, and the free surface in particular, are interpolated by means of biquadratic shape functions when surface tension is taken into account, while a bilinear interpolation is used otherwise. For predicting the free surface location, we use an implicit technique similar to that developed by Kistler and Striven in Ref. 32; the kinematic equation (8) is integrated either by means of the Galerkin technique or a streamline-upwind/Petrov-Galerkin formulation [33]. The latter will be used when surface tension is neglected. The set of nodal unknowns for T,, Y, p and for the free surface location is computed by means of a Newton iterative algorithm, which has the major advantage of quadratic convergence. 3. The rod climbing of an upper-convected Maxwell liquid As a first example of axisymmetric swirling flows, we focus on the rod climbing or Weissenberg effect of an upper-convected Maxwell liquid. In Fig. 1, we display a view of the flow problem in a cross-sectional plane containing the axis of symmetry. A fluid contained in a cylindrical vessel is set into motion by means of a rotating rod, the axis of which coincides with that of the vessel. In the present computations, the radius of the fluid container is ten times that of the rod. As far as boundary conditions are concerned, the fluid sticks to the wall and to the rod which rotates at an
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
I
109
wall
Fig. 1. Rod climbing: flow description
and boundary conditions.
angular velocity a. Along the axis, we impose the classical symmetry boundary conditions, while the kinematic equation (6) is applied along the stress-free upper surface. The flow domain is closed, so that no boundary conditions are required for the viscoelastic stress tensor Ti. In the present section, we focus on the interaction between gravity and normal stresses. Therefore, we take Re = 0 and Cu + 00. For the finite-element discretization, we use the meshes displayed in Fig. 2. They are characterized by a fine definition in the vicinity of the free surface, and near the rod. This is especially true for WEISS~, where a large number of small quadrilateral elements are located in that region, in order to capture both the small secondary motions which are expected there [lo] and the free surface shape. Our first objective in this section is to assess mesh convergence of the numerical results. Both MIX1 [21] and 4 x 4-SU [22] algorithms are used with WEISS& while the first is used with WEISS~ in view of its lower computational cost. Table 1 displays the number of degrees of freedom for the meshes of Fig. 2 with the corresponding algorithms. Beside the free surface shape, we focus attention on the evolution of two quantities as a function of We. (i) The first is the climbing h, as displayed in Fig. 1; for the sake of generality, we use a non-dimensional climbing h/R. (ii) Numerical results will show the development of two secondary vortices, one of them being located under the bulge. The ratios of intensity of the secondary vortices to the azimuthal flow rate are also of interest. In Figs. 3(a)-(d), we plot both the climbing h and the vortex intensities as a function of We for 99 = 0.01. These curves are obtained by means of the three strategies described above and indicate a fairly good mesh
B. Debbaut and B. Hocq /.J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
Fig. 2. Rod climbing: finite-element
meshes.
convergence, even in a range of We numbers where h differs from a quadratic dependence with respect to We. The graphs in Fig. 3 also show that both h and vortex intensities are increasing functions of We. Furthermore, one sees in Fig. 3(b) that for low We, the climbing h is a quadratic function of We as already noted [8-131 and matches the analytically predicted values. This emphasizes the reliability of our numerical predictions. In Refs. 8-13 we find an analytical expression for the free surface shape, which is valid for low angular velocity R. In the particular case of the
TABLE 1 Degrees of freedom for the three numerical strategies Field
WEISSI,
MIX1
Tl u P
10146 5073 446
39486 5073 446
18126 9063 786
Total
15665
45005
27975
WEISSI,
4 x 4-su
WEISS2,
MIX1
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
0.12
h
x
:
111
a
0.3
ox
8
q
0.12
@O
* . 0.4
. a
We0
W&
0.0 0
5
10
15
a
20 x 10-3
x 10”
0
b
0.2
0.4
0.6
0.8
1
x 10”
x 10-4
5 vortex
intensity
4
x B d
3
dp on
2
a”
q
. 1
9 ox
1
.
n+ We
0 ,.m= 0.0
0 0.04
0.08
C
0.12
0.0
0.04
0.06
0.12
d
Fig. 3. (a), (b) Climbing h as a function of We*; (c) intensity of the large vortex as a function of We, for 9 = 0.01; (d) intensity of the small vortex as a function of We, for 6 = 0.01. (+), MIX1 algorithm with WEISSI; (X 1,4X4-W algorithm with WEISSI; (o), MIX1 algorithm with w&s2; ( -1, analytical prediction.
creeping flow of an upper-convected 8-13 reduces to h(r) -= R
Maxwell fluid, the expression in Refs.
2We2 ‘B( r/R)4 ’
(11)
in non-dimensional form. In Fig. 4 we have plotted the numerical prediction of the free surface shape for We = 0.005, 0.01 and 0.015, together with the analytical counterpart given by (11). Excellent agreement is obtained; for low We, the free surface indeed shows a shape following l/r4. In Fig. 5 we display the free surface shape and streamlines of the secondary motions in the neighbourhood of the bulge obtained with WEIW, at various values of We for 9 = 0.01. The growth of a small vortex under
112
.05
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 1
h(r) R t,
*04 .-a We = .015
.a3
.’
‘,
r*
.02
.Ol
.OO 1.0
1.8 r 2.0 R Fig. 4. Free surface shape obtained at various values of We for 6 = 0.01. (+I, MIX1 algorithm with WEISSI; (X), 4 X4-SU algorithm with WEISSI; (01, MIX1 algorithm with ), analytical prediction. WEISSZ; (---
We = 0.025
We = 0.0665
1.2
1.4
1.6
We = 0.0644
We = 0.099
Fig. 5. Free surface shape and streamlines obtained at various values of We for 9 = 0.01. A fluid particle moves counterclockwise along a positive streamline, and clockwise along a negative streamline.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
113
the bulge is indeed observed. As no inertia is involved in the present simulation, the two vortices are generated and driven solely by viscoelasticity. In the large vortex, fluid particles move counterclockwise, while they move clockwise in the small vortex under the bulge and thus climb the rod as already sketched by Yoo et al. [lo]. The separating line between the two opposed secondary motions starts near to a point where the free surface changes curvature. Predicting the location of this point with accuracy seems to be difficult, as it would require an extremely refined mesh. In particular, it is the re-entrant point of the opposite boundary streamlines of both secondary vortices. In the vicinity of this point, at high We, rapid slope changes in the free surface can occur within an extremely short distance as also sketched by Yoo et al. [lo] and experimentally shown in Refs. 8, 14 and 20. At high We, the free surface experiences large deformations and also shows a loss of smoothness in its finite-element discretization, which pollutes the numerical solutions and leads to unrealistic numerical predictions. If a Galerkin technique is used for integrating eqn. (61, spurious oscillations in the free surface appear much earlier. Despite this unrealistic loss of smoothness in the free surface, numerical results can still be obtained. In view of their lack of interest, they are not displayed here. Also, we were not interested in capturing bifurcation points that we believe are numerical artefacts in view of the spurious oscillations. Several possible causes can explain this loss of smoothness which has been observed with all our strategies. Surface tension has been neglected in the present computations, which would have led to a reduction of the changes in curvature. Thus we have investigated here a more critical situation than that with surface tension. Next, qualitatively speaking, pictures in Refs. 8 and 20 showed similar experimental bulge shapes characterized by a rapid slope change over a short distance at high rod angular velocity, while a time-periodic instability was evidenced at higher angular velocity. A transient algorithm combined with an adaptive mesh technique (in the region of rapid change in the free surface) should be used. We must point out here that the streamlines displayed in Fig. 5 are not fluid trajectories; rather they are traces of trajectories in a cross-sectional plane containing the axis of symmetry. Along those traces, a fluid particle moves extremely slowly, and the major flow remains the swirling one. The graphs of Figs. 3(c) and Cd) showing the vortex intensities emphasize this statement. In Fig. 6, we have plotted the contour lines of the azimuthal velocity component ug and the pressure field for We = 0.0644 and B = 0.01. We observe in particular that the azimuthal velocity component is slightly higher in the bulge than along the rod. We believe that it originates from a
114
B. Debbaut and B. Hocq/J.
@ E
Non-Newtonian Fluid Mech. 43 (1992) 103-126
rl
,005
u, ==
I
&i
1.
.Ol
I ,i
b
Fig. 6. Contour lines of the azimuthal velocity component (b), for We = 0.0644 and S? = 0.01.
ug (a>, and of the pressure field p
quasi-rigid motion experienced by the fluid in a region where the shear stress is virtually zero. As far as the pressure field is concerned, its contour lines reveal the expected results. The free surface is an isobar. Far from the rod, the pressure contour lines are controlled by gravity. Close to the rod, they shift from their hydrostatic counterpart and move upwards. This arises from the azimuthal extra-stress component (or normal stress) which rapidly increases towards the rod, as stated in Refs. 2, 3 and 17. If one follows a horizontal line towards the rod, the pressure increases. This is illustrated in Fig. 7, where we give the pressure profile along the radial line passing through the free surface at a long distance from the rod, at We = 0.0644 and 9 = 0.01.
x 1o-3
12
8
5 4 2 3 1 Fig. 7. Pressure profile along the radial line passing through the free surface at long distance from the rod, at We = 0.0644, for 22 = 0.01.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
4. The rod climbing of a Johnson-Segalman
115
liquid
Beavers and Joseph [8] reported theoretical predictions on rod climbing for a fourth-order Rivlin-Ericksen fluid. For low values of a, their predictions showed a good quantitative agreement with experimental measurements done with STP (a 26.6% polyisobutylene solution in petroleum oil). Our purpose here is to check the validity of finite-element predictions for such a flow. In Ref. 8, the diameter of the cylindrical vessel is 30.5 cm; the height of fluid at rest is 7.7 cm. For the rotating rod, we have selected a radius of 0.635 cm. Material data for STP are estimated or given (at a temperature of 25.5”C): p =
0.89 g cmp3
q = 146 P,
’
(1 - 4)h = 0.0162 s,
5 = 0.275, I = 30.9 dyn cm-‘, g = 981 cm ss2.
(12) In view of the non-vanishing second normal stress difference shown by the STP, a Johnson-Segalman model is selected; for 4 in (3) and (12), we have taken the value l/9. From (121, it follows that B = 0.0693. In their experiments, Beavers and Joseph [S] coated the rod with Scotchgard, so that at the contact point the free surface is normal to the rotating rod. In the present simulation, we impose that the free surface has a vanishing slope on the rod. For low angular velocity a, when normal stress differences, inertia forces and surface tension are taken into account, the climbing h can be predicted analytically as follows [g-13], in a dimensional form:
(13)
P = (1 - 4)rlA - 2(1- 4b7~5
(14)
is the so-called climbing constant and is a linear combination of the first and second normal stress coefficients. Equation (13) predicts a quadratic dependence of the climbing h upon R, and thus upon We.
116
B. Debbaut and B. Hocq /J. Non-Newtonian Fluid Mech. 43 (1992) 103 - 126
Fig. 8. Rod climbing: finite-element
mesh.
For the numerical simulation of the present flow as experimentally investigated [8], we use the finite-element mesh shown in Fig. 8 together with the 4 x 4-SU finite-element technique [22]. This mesh is characterized by 30438,3915 and 345 degrees of freedom, respectively for T,, v and p. A fine discretization is specified in the vicinity of the contact point of the rotating rod and the free surface. Such a fine mesh has been selected in order to obtain a good description of the flow patterns. The flow region located below the rod (still considered in Fig. 2) has been removed as in the experimental device described in Ref. 8. To avoid a discontinuity in the azimuthal velocity component at the left bottom point in Fig. 8, vanishing tangential forces are imposed over the last element. Numerical simulations have been performed up to We = 0.360. Figure 9 displays the climbing h/R as a function of We2. For low We, we find that our numerical predictions compare well with the experimental measure-
.6
h
We2
.oo
.04
a3
.12
.16
Fig. 9. Climbing h as a function of We*. (predictions; (01, experimental data.
), analytical prediction;
( + 1, numerical
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
117
y .15
.lO
.05
1.0
1.4
1.8
2.2
Fig. 10. Free surface shape obtained numerical prediction.
R at We = 0.195. (o), experimental
data [8]; (-),
ments [S], and with the analytical predictions from eqn. (13). In particular, we note that h is quadratic in We, as was shown experimentally [S]. At higher We, however, a slight discrepancy arises among these results. The uncertainty on some material parameters in (12), and the selection of a Johnson-Segalman model with one single relaxation time and a purely viscous stress component, are possible causes for such a discrepancy. In Fig. 10, we have plotted the free surface profile obtained at We = 0.195. A comparison with the experimental measurements carried out by Beavers and Joseph [8] is also shown. Again, a slight discrepancy is observed, similar to that mentioned above. A complex development of flow patterns is observed as a function of We. This is illustrated in Fig. 11, where we have displayed the streamlines in a cross-sectional plane containing the axis of symmetry for various values of We. We believe that such a complex development of secondary motions results from the opposing interaction of normal stress differences and inertia. At low We, we essentially observe a large vortex driven by inertia with two centres, which is more intense than the small vortices located on the bottom of the vessel and under the free surface, which are driven by the normal stresses. The vortex intensities have been normalized with respect to the azimuthal flow rate; one finds that the secondary motions are extremely weak as compared to the primary flow. At low We, the direction of the fluid particles clearly indicates the physical origin of these vortices: the fluid particles move clockwise in a vortex driven by inertia and are identified by negative values of the stream function, while they move
118
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
5.0 IO*
l:o 105
We=.045
We= .109
-i.oIO-5
11010-S
We= .142
! We= i
4.0
10-J
.195
-l.o 1o-5
5.0 lD_5
A.5 10”
We= 240
We = .298
-0.25 10-4
.A.5 104
i
/ I
! We= .322
We = .360
Fig. 11. Streamlines obtained at various values of We. A fluid particle moves counterclockwise along a positive streamline and clockwise along a negative streamline.
counterclockwise in a vortex driven by viscoelasticity, where they are identified by positive values of the stream function. These viscoelastic vortices grow with We; this is especially true for the one located under the free surface. This last vortex eventually becomes
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
119
more intense than the original inertial vortex. Meanwhile, a third small viscoelastic vortex is generated under the bulge, the intensity of which is very low. We note here that the two viscoelastic vortices under the free surface show a behaviour similar to that seen in the previous section. At higher We, the original inertial vortex is split into two parts. Both these parts are progressively squeezed by the larger viscoelastic vortex, which results from the assembly of the first two viscoelastic vortices and which still grows. Independent of this phenomenon, the third viscoelastic vortex located under the bulge grows slowly. Calculations have not been pursued beyond We = 0.360, as convergence becomes difficult to achieve. Moreover, as previously shown in Fig. 9, a discrepancy occurs between numerical results and experimental measurements of the climbing at high We. These problems are not yet fully understood. A possible cause can be found in the free surface, which shows at least two stagnation points for high We. 5. The Quel1efik.t
Another axisymmetric swirling flow of interest is the so-called QueZZeffekt in Refs. 4 and 18. Figure 12 illustrates this flow. A cylindrical vessel, the radius of which is R, is filled with a fluid up to filling height H. As a disk, located at the bottom of the container, rotates around its axis, the fluid is set into motion and flows upwards along this axis, like a source. A bulge is then generated at the upper surface of the liquid. This bulge can be characterized by its vertical displacement h along the axis of symmetry. For simulating this flow, we have selected the boundary conditions as
wall
:
Fig. 12. Quelleffekt flow description
and boundary conditions.
120
B. Debbaut and B. Hocq/
1
I
I
I
Fig. 13. Quelleffekt: finite-element
I
J. Non-Newtonian
Fluid Mech. 43 (1992) 103-126
I lttl
mesh.
follows. The fluid sticks to the walls of the container and to the rotating disk; symmetry boundary conditions are imposed along the axis, while eqn. (6) is applied on the upper surface which is stress free. To avoid possible difficulties arising from a discontinuity in the azimuthal velocity boundary condition at the lower right-hand side corner of the domain, vanishing tangential forces are imposed on the bottom within a distance of O.OlR from the wall. As for the rod-climbing simulations of the two previous sections, no boundary conditions are specified for the viscoelastic stress tensor T,, since the flow domain is closed. Since our investigations in Section 3 assess the reliability of the results, we have performed our computations with the finite-element mesh of Fig. 13. In particular, it is characterized by 15006, 1953 and 176 degrees of freedom for T,, v and p respectively. It also contains relatively small elements along the boundary. In this section, we again use the finite-element technique developed by Marchal and Crochet [22]. Biihme et al. [18] reported numerical and experimental results for the Quelleffekt of a 2.5% aqueous polyacrylamide solution. In their experimental device, the radius of the vessel was 4.5 cm. Fluid material data are also found in Ref. 18: p =
1.005 g cme3
~=76P (1 -@A
= 0.776 s,
%$ = 0.22, g = 981 cm se2.
(15)
In view of the non-vanishing second normal stress difference, we have selected a Johnson-Segalman model; for 4, we have taken the value l/9. From (15), it follows that 9 = 5.09. Also, in order to observe the effects of the sole first normal stress difference, we also performed computations with an Oldroyd-B model for which we selected 4 = 0.1; therefore % = 5.03.
B. Debbaut and B. Hocq / J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
121
model, at We = 0.323, and H = 1.
Fig. 14. Streamlines obtained with the Johnson-Segalman
For the experiments discussed in Ref. 18, the height of filling H of the fluid in the cylindrical vessel was an independent parameter. For H ranging from 0.5 to 1.25, numerical predictions and experimental measurements carried out by Biihme et al. 1181quantitatively coincide. Figure 14 shows the contour lines obtained with the Johnson-Segalman model at We = 0.323, for H = 1. Our predictions show a very good quantitative agreement with that obtained in Ref. 18. In Fig. 15, we plot the displacement h of the upper surface along the axis of symmetry as a function of H for both Johnson-Segalman and -3
-3
x 10
x 10
16
16
8
8
4
4
0
0 0.5
a
0.9
1.3
0.5
1.7
b
0.9
1.3
#
1.7
Fig. 1.5. Displacement h of the upper surface along the axis of symmetry as a function of H. predictions of Ref. 18; (+), (a) Oldroyd-B model; (b) Johnson-Segalman model. ( -_), present numerical predictions.
122
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103-126 x
10
-3
2.5
8
2.0
6
1.5 4 1.0 2
0.5 0.0
0 0.0
a
0.2
0.4
0.6
0.0
0.04
0.08
0.12
b
Fig. 16. Graph of h as a function of We *, for H = 0.5. (a) Johnson-Segalman model; (b) Oldroyd-B model. ( + ), Numerical predictions; (), predictions of Ref. 18.
Oldroyd-B models, at We = 0.323 and 0.319, respectively. We also display the predictions of Ref. 18, which showed good agreement with experimental measurements since both first and second normal stress differences were taken into account. Our predictions match well those of Ref. 18. One observes in Fig. 15 that h decreases as H increases, and thus as the effects of the rotating disk move away. Moreover, we find that selecting the Oldroyd-B model leads to a larger displacement of the upper surface along the axis of symmetry, as was predicted [18]. This results from the vanishing second normal stress difference of this model. Figure 16 displays the evolution of h as a function of We2 for both Johnson-Segalman and Oldroyd-B models for H = 0.5. Comparison with results of Ref. 18 is also displayed, and again shows excellent agreement. One observes here that for low values of We, h is a quadratic function of We. Again one finds that selecting the Oldroyd-B model, and thus removing the second normal stress difference, leads to a larger deformation of the upper free surface. At high We, h is no longer quadratic in We, as would be expected, since the predictions of Ref. 18 are only valid for low We. In Fig. 17, we display the free surface profiles for various values of H, obtained with both Johnson-Segalman and Oldroyd-B models, at We = 0.323 and 0.319 respectively. These profiles show good agreement with the predictions of Ref. 18. However, for low H, they are characterized by a slight loss of smoothness near the wall of the vessel. The finite-element discretization is probably not sufficient in this region to capture the abrupt deformation experienced by the free surface. In Fig. 18, we plot the intensity of the secondary motion as a function of We, which has been normalized with respect to the azimuthal flow rate.
B. Debbaut and B. Hocq/ .I. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 0.02
0.02
7
y
y 0.01
0.01
0.00
0.00
-0.01
-0.01
-0.02
a
123
-0.02 0.0
0.2
0.4
0.6
0.6
1.0
b
o.0
0.2
0.4
0.6
0.6
1.0
Fig. 17. Free surface shape for various values of H. (a) Johnson-Segalman We = 0.323; (b) Oldroyd-B model, at We = 0.319.
model, at
Again, we find that the secondary motion is extremely weak as compared to the primary one, so that a fluid particle must do a large amount of rotation around the axis of the vessel before recovering its initial location in the (r, z) plane. 6. Conclusions A finite-element technique has been used for simulating axisymmetric swirling flows of viscoelastic liquids of the differential type. Convergence of the results has been proved on the basis of different meshes and numerical algorithms.
0.12 Vortex
intensity
a
0.08
0.04
We
0.00 0.0
0.2
0.4
0.6
0.8
Fig. 18. Vortex intensity as a function of We. (a) Johnson-Segalman model.
model; (b) Oldroyd-B
124
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
Our technique has been successfully applied to both rod-climbing flow and the Quelleffekt. Our present results on rod climbing and surface displacement agree well with previous predictions. Moreover, we have been also able to quantify the secondary motion intensity for both these flows and to show that the secondary motions are extremely weak as compared to the main azimuthal flow. For the rod-climbing flow, a full simulation of a realistic experiment has been performed, including both inertia terms and surface tension. The development of complex flow patterns has also been shown, as the rod angular velocity increases. In particular, we have successively found that two viscoelastic vortices grow, split the large original inertial vortex into two parts, and are eventually combined into a larger viscoelastic vortex which progressively squeezes the remaining parts of the split inertial vortex. Numerical simulations performed for the Quelleffekt have shown the opposed effects of the first and second normal stress differences on the surface displacement. We have shown that, as far as bulge displacement is concerned, the second normal stress difference acts against the first one. Further flow problems will be investigated in later work. In particular, it is meaningful to perform transient computations which could possibly predict flow instabilities at high Weissenberg numbers. Extension to differential viscoelastic models with multiple relaxation times will also be attempted. Acknowledgement The authors wish to acknowledge the collaboration of Dr. J. Rosenberg from Polyflow S.A. who improved the final version of the English text. References 1 F.H. Garner and A.H. Nissan, Rheological properties of high viscosity solutions of long molecules, Nature, 158 (1946) 634-635. 2 K. Weissenberg, A continuum theory of rheological phenomena, Nature, 159 (1947) 310-311. 3 H.A. Barnes, J.F. Hutton and K. Walters, An Introduction to Rheology, Elsevier, Amsterdam, 1989, Chapter 5. 4 R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, John Wiley, New York, 2nd edn., 1987, pp. 62-65. 5 K. Walters, Rheometry, Chapman and Hall, London, 1975. 6 J. Serrin, Poiseuille and Couette flow of non-Newtonian fluids, Z. Angew. Math. Mech., 39 (1959) 295-299. 7 H. Giesekus, Einige Bemerkungen zum FlieBverhalten elasto-viskoser Fliissigkeiten in stationaren Schichtstromungen, Rheol. Acta, 1 (1961) 404-413.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 8 G.S. Beavers and D.D. Joseph, The rotating rod viscometer,
9
10 11 12 13 14
15 16 17 18 19
20 21 22 23 24 25
26 27 28 29
125
J. Fluid Mech., 69 (1975) 475-511. D.D. Joseph and G.S. Beavers, Free surfaces induced by the motion of viscoelastic fluids, in R.S. Rivlin (Ed.), The Mechanics of Viscoelastic Liquids, AMD 22, ASME, New York, 1977, pp. 59-99. Y. Yoo, D.D. Joseph and G.S. Beavers, Higher-order theory of the Weissenberg effect, J. Fluid Mech., 92 (1979) 529-590. G.S. Beavers, J.Y. Yoo and D.D. Joseph, The free surface on a liquid between cylinders rotating at different speeds. Part III, Rheol. Acta, 19 (1980) 19-31. D.D. Joseph, Fluid Dynamics of Polymeric Liquids, Springer-Verlag, New York, 1990, Chapter 17. D.D. Joseph, G.S. Beavers, A. Cers, C. Dewald, A. Hoger and P.T. Than, Climbing constants for various liquids. J. Rheology, 28 (1984) 325-345. H.H. Hu, 0. Riccius, K.P. Chen, M. Arney and D.D. Joseph, Climbing constant, second-order correction of Trouton’s viscosity, wave speed and delayed die swell for Ml, J. Non-Newtonian Fluid Mech., 35 (1990) 287-307. D.D. Joseph, Remarks on inertial radii, persistent normal stresses, secondary motions and non-elastic extensional viscosities, J. Non-Newtonian Fluid Mech., 32 (1989) 107-114. B. Debbaut and M.J. Crochet, Extensional effects in complex flows, J. Non-Newtonian Fluid Mech., 30 (1988) 169-184. C.T. Hill, Nearly viscometric flow of viscoelastic fluids in the disk and cylinder system. II. Experimental, Trans. Sot. Rheol., 16 (1972) 213-245. G. BBhme, R. VoR and W. Warnecke, Die freie Oberflache einer Fliissigkeit iiber einer rotierende Scheibe, Rheol. Acta, 24 (1985) 22-33. A. Delgado, S. Berg, R. Kroger and H.J. Rath, Theoretical considerations on the steady rod climbing at vanishing influence of gravity, in D.R. Oliver (Ed.), Proc. of the Golden Jubilee of the British Society of Rheology and 3rd European Rheology Congress, Edinburgh, Elsevier, 1990, pp. 129-131. D.D. Joseph, Stability of Fluid Motion II, Springer tracts in Natural Philosophy, Berlin, 1976, pp. 232-237. M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of non-Newtonian Flow, Elsevier, Amsterdam, 1984. J.M. Marchal and M.J. Crochet, A new mixed finite element for calculating viscoelastic flow, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. B. Debbaut, J.M. Marchal and M.J. Crochet, Numerical simulation of highly viscoelastic flows through an abrupt contraction, J. Non-Newtonian Fluid Mech., 29 (1988) 119-146. M.J. Crochet, Numerical simulation of viscoelastic flow: a review, Rubber Chem. Technol., 62(3) (1989) 425-455. R. Keunings, Simulation of viscoelastic fluid flow, in C.L. Tucker III (Ed.), Fundamentals of Computer Modelling for Polymer Processing, Carl Hanser Verlag, Munchen, 1989, pp. 403-469. S. DuPont and M.J. Crochet, Swirling flows of viscoelastic fluids of the integral type in rheogoniometers, Chem. Eng. Commun., 53 (1987) 199-221. J.G. Oldroyd, On the formulation of rheological equation of state, Proc. R. Sot., London, Ser. A, 200 (1950) 523-541. M.W. Johnson, Jr., and D. Segalman, A model for viscoelastic fluid behavior which allows non-affine deformation, J. Non-Newtonian Fluid Mech., 2 (1977) 255-270. X. Georges et B. Hocq, Simulation numerique de l’effect Weissenberg, Memoire, Unite de Mecanique Appliquee, Louvain-la-Neuve, 1986.
126
B. Debbaut and B. Hocq/J.
30 0. Wambersie,
Non-Newtonian Fluid Mech. 43 (1992) 103-126
Simulation numerique de i’extrusion tridimensionnelle par une methode pseudo-transitoire basCe sur les gradients conjugues, Ph. D. Thesis, Unite de Mecanique Appliquee, Louvain-la-Neuve, 1990. 31 R. Keunings, Progress and challenges in computational rheology, Rheol. Acta, 29 (1990) 556-570. 32 S.F. Kistler and L.E. Striven, Coating flows, in J.R.A. Pearson and SM. Richardson (Eds.), Computational Analysis of Polymer Processing, Applied Science Publishers, London, 1983, pp. 243-299. 33 A.N. Brooks and T.J.R. Hughes, Streamline-upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible NavierStokes equations, Comput. Methods Appl. Mech. Eng., 32 (1982) 199-2.59.
103
Elsevier Science Publishers B.V., Amsterdam
On the numerical simulation of axisymmetric swirling flows of differential viscoelastic liquids: the rod climbing effect and the Quelleffekt Benoit Debbaut 1,2 and Bernard Hocq 2 ’ Unite’de Mtkanique Appliquke, Universite’ Catholique de Louvain, Place du Levant 2, B-1348 Louvain-la-Neuve (Belgium) 2 Polyflow S.A., Place de I’Universite 16, B-1348 Louvain-la-Neuve (Belgium) (Received May 23, 1991; in revised form November 22, 1991)
Abstract
A finite-element investigation is presented of some Weissenberg effects occurring in swirling flows of viscoelastic liquids in geometries of revolution. Numerical results are shown for the rod-climbing flow and the so-called Quelleffekt of viscoelastic liquids of the differential type. For both types of flow, comparisons with earlier numerical or analytical predictions and some experimental measurements are also shown. Keywords: rod climbing; Weissenberg
effect; viscoelastic liquids; swirling flows; numerical
simulation.
1. Introduction
Early experimental observations by Garner and Nissan [l] with a solution of rubber in benzene showed a rod-climbing effect. This led Weissenberg to introduce the hypothesis of pull, or normal stress, acting along streamlines [2]. Such a pull appears to be responsible for the so-called Weissenberg effect (or rod climbing). Indeed, as a physical explanation of this effect, it was suggested in Ref. 3 that the normal stress “acts like a hoop stress around the rod” and in Ref. 4 “strangulates the fluid, and forces it Correspondence to: B. Debbaut, Unite de Mecanique AppliquCe, Universite Louvain, Place du Levant 2, B-1348 Louvain-la-Neuve, Belgium.
0377-0257/92/$05.00
Catholique
0 1992 - Elsevier Science Publishers B.V. All rights reserved
de
104
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103-126
inwards against the centrifugal force and upwards against the gravitational force”. Such a flow belongs to the class of axisymmetric swirling flows which, from a kinematic viewpoint, are characterized by three velocity components depending upon two space variables (r, z) in a cylindrical reference frame. The rheometrical interest of this class of flows has already been discussed [3-51. In the literature, several papers have been devoted to analytical developments using perturbation theory, for rod climbing [6-121, indicating a relationship between rod-climbing measurements and the viscoelastic properties of the liquid. On the basis of analytical developments for Couette flow with a free surface, Serrin [6] and Giesekus [7] discussed the origin of the fluid climbing on the inner cylinder. In Ref. 7, Giesekus also introduced the concept of positive and negative Weissenberg effects. In Refs. 8-12, the authors showed that a knowledge of the climbing was sufficient for determining an important viscoelastic parameter of a simple fluid. Experimental measurements carried out by Beavers and Joseph [S] on STP (polyisobutylene in petroleum oil) showed good quantitative agreement with analytical predictions. Joseph et al. [13] have tabulated experimental rod-climbing measurements for several fluids. Yoo et al. [lo] also analysed the development of secondary motions as a function of the rod angular velocity. Hu et al. [14] and Joseph [15] have discussed the control of secondary motions by the extensional viscosity, which was already known to be important in vortex development in contraction flows [16]. Among experimental works for various fluids in a disk and cylinder system, Hill [17] noted a reversal in the direction of the secondary motion as the fluid becomes “more conductive to viscoelastic secondary flows”. A phenomenon related to the Weissenberg effect is the Quelleffekt investigated by Bijhme et al. [18]. Here, the viscoelastic liquid is contained in a cylindrical vessel. The fluid is set into motion by a rotating disk located at the bottom of the vessel. It flows upwards along the axis of symmetry, as a source (QueU), and produces a bulge in the free surface around the axis. This effect can be understood as a “rod-climbing without a rod” [4]. In Ref. 18, the authors analysed the effects of inertia and normal stress differences upon bulge shape. Also, the predicted results showed good agreement with experimental measurements for both silicone oil and a 2.5% polyacrylamide (PAA) aqueous solution. PAA, which is much more elastic than silicone oil, is also characterized by a non-constant shear viscosity with a small plateau and by a non-constant normal stress coefficient with a larger plateau. However, the authors in Ref. 18 displayed results obtained in the range of constant normal stress coefficient and
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
105
showed that the bulge shape is not affected by a non-constant viscosity. An important feature, seen both theoretically and experimentally, is that climbing and axial bulge deformation are quadratic in the angular velocity of the rotating device (rod or disk), for low angular velocity [g-13,18]. In view of the air-fluid interface in the experiments, surface tension effects have also been investigated [8,12,18,19]. In particular, in experiments with STP, Beavers and Joseph [8] removed the static contact angle by coating the rod with Scotchgard; with coated and uncoated rods, the corresponding theoretical predictions still matched the experimental measurements. Delgado et al. [19] have analysed the limiting case as the gravitational force vanishes. In their experiments with PAA, Bijhme et al. [18] found that the influence of surface tension is of the order of a few percent. Finally, we should mention some studies on the stability of rod-climbing flow. As the angular velocity increases, the flow can show a time-periodic behaviour, the breathing instability, which is discussed in Refs. 8 and 20. Our purpose in the present paper is to extend the capabilities of numerical algorithms widely used for 2-D flows of differential viscoelastic liquids [21-251 to axisymmetric swirling flows. Numerical simulations of such flows have already been done for viscoelastic fluids of integral type [26]. Although a large variety of axisymmetric swirling flows of viscoelastic liquids are discussed in Ref. 12, we shall investigate here the finite-element simulation of two flows with a free surface: the rod climbing and the Quelleffekt of both Oldroyd-B and Johnson-Segalman liquids [27,28]. Such flow problems have been selected, not only in view of their rheometrical interest, but also because analytical and experimental results already exist [8,18]. Furthermore, exp”erimenta1 measurements were obtained in a range of shear rates for which normal stress coefficients are constant. In particular, we shall focus on the free surface shape and on the development of secondary motions. The present work follows an early finite-element investigation done by Georges and Hocq [29] on the rod-climbing effect for both Newtonian and second-order fluids. The equations and the finite-element discretization are discussed in Section 2. In Section 3, we analyse the rod-climbing flow of an upper-convected Maxwell liquid and identify some interesting features. Sections 4 and 5 are devoted to rod climbing and Quelleffekt of Oldroyd-B [27] and Johnson-Segalman [28] liquids; comparisons with earlier predictions and experimental results are discussed. 2. Equations
and discretization
Let us denote the rate of deformation d = $7~
+ bT),
tensor d as follows:
(1)
106
B. Debbaut and B. Hocq/J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
where v is the velocity field. For the models considered stress tensor T can be decomposed as follows: T= Tl + T2,
here, the extra(2)
where Tl and T2 are, respectively, the viscoelastic and the Newtonian components to the extra-stress tensor T. These components obey the following constitutive relationships:
where q is the total shear viscosity, A is the relaxation time, and the symbols v and A denote the upper-convected and 1ower;convected derivative operators, respectively. In (3), 5 is a material parameter controlling the ratio of the second to the first normal stress differences, while C#Jis the ratio of the Newtonian solvent viscosity to that of the solution. As 5 vanishes, one recovers the Oldroyd-B model [27]; furthermore, as 4 and 5 vanish, one recovers the upper-convected Maxwell liquid. The constitutive relationships (3) are coupled with the momentum and incompressibility equations, i.e.: -Vp+VT+pg=pvVv,
(4)
v*v=o,
(5) where p is the pressure field, p is the fluid density and g is the gravity vector. For the class of axisymmetric swirling flows that we consider here, the unknown fields T,, v and p depend only upon two space variables in a cylindrical system of axes (e,, eZ), and do not vary in the azimuthal direction e,. However, we must point out here that the fields Tl and v now have 6 and 3 components respectively, instead of 4 and 2 as in the case for 2-D axisymmetric flows [21]. In terms of numerical simulation, which will be discussed below, swirling flow problems can be handled with the same computational algorithms that were used for 2-D flow situations. The computational cost increases due to the additional amount of components for Tl and v, but remains much lower than that required for a full 3-D simulation [30]. In addition to eqns. (l)-(5), additional relationships must be satisfied. Indeed our flow problems involve free surfaces, the location of which is a priori unknown. The kinematic condition on free surfaces can be written as follows: v*n=o,
(6)
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (192)
103- 126
107
where it is the normal vector to the free surface. In some simulations, surface tension a, will he taken into account. This obeys the relationship: 1 -+’
q=r i
R1
I n
R,
’
(7)
where I’ is the surface tension coefficient, while R, and R, are the principal radii of curvature of the free surface. Our flow problems occur in closed volumes, so that any solution must satisfy the global form of volume conservation: /V
dV-
V, = 0,
where V,, is the initial volume at rest. 3 characterize the flow, we use non-dimensional numbers. Therefore, in addition to the given material and physical data q, A, p, I and g, we select a characteristic dimension R and angular velocity R of the rotating device. We define the Reynolds, Weissenberg, Stokes and capillary numbers, respectively, as follows: pflR2 Re = rl ’ We=hlR, 770 St = PgR ’ 77fiR Ca = r *
Although Re is generally very small, acceleration terms in (4) will be taken into account in some simulations. It has been shown [7,18,26] that inertia forces generate a secondary motion opposed to that of viscoelasticity. Moreover, although the Stokes number St may be large, gravity forces cannot be neglected. Indeed, both rod climbing and the Quelleffekt result from viscoelastic normal stresses acting against gravitational forces. A further dimensionless number, 9, can be obtained as the ratio of the Weissenberg number to the Stokes number: $+_
QgR 77 *
(10)
This number is independent of the angular velocity CI, so that, for a given fluid with its material properties A, p and q, in a given environment specified by g, 9 characterizes the flow geometry.
108
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) IO3- 126
The numerical simulation is performed by means of two finite-element algorithms. (i) The algorithm developed by Marchal and Crochet [22], which combines 4 X 4 sub-element interpolations for the viscoelastic extrastress tensor Z’i with a streaml@e upwgding (SU) technique applied to the convective term appearing in ?‘r and T1. Mesh-convergence properties for the velocity field and streamlines were already proved in Refs. 22 and 23. (ii> The earlier MIX1 algorithm [21] with biquadratic interpolation for the viscoelastic extra-stress is also used for several reasons. As will be shown, interesting solutions are obtained at relatively low values of We. Moreover, our flow problems are smooth enough in the regions of interest, so that the technique MIX1 can be successfully used, as suggested in Ref. 31. This technique is also cheaper from a computational viewpoint for refined meshes. The use of two different algorithms will also allow us to assess the convergence of the solutions, although solutions at higher We will be obtained only with the 4 X 4 SU technique. In this paper we do not intend to track possible bifurcation points; thus, calculations will be pursued as long as interesting results can be obtained. For both techniques, biquadratic and bilinear interpolants are used for the velocity and pressure fields, respectively. The finite-element mesh, and the free surface in particular, are interpolated by means of biquadratic shape functions when surface tension is taken into account, while a bilinear interpolation is used otherwise. For predicting the free surface location, we use an implicit technique similar to that developed by Kistler and Striven in Ref. 32; the kinematic equation (8) is integrated either by means of the Galerkin technique or a streamline-upwind/Petrov-Galerkin formulation [33]. The latter will be used when surface tension is neglected. The set of nodal unknowns for T,, Y, p and for the free surface location is computed by means of a Newton iterative algorithm, which has the major advantage of quadratic convergence. 3. The rod climbing of an upper-convected Maxwell liquid As a first example of axisymmetric swirling flows, we focus on the rod climbing or Weissenberg effect of an upper-convected Maxwell liquid. In Fig. 1, we display a view of the flow problem in a cross-sectional plane containing the axis of symmetry. A fluid contained in a cylindrical vessel is set into motion by means of a rotating rod, the axis of which coincides with that of the vessel. In the present computations, the radius of the fluid container is ten times that of the rod. As far as boundary conditions are concerned, the fluid sticks to the wall and to the rod which rotates at an
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
I
109
wall
Fig. 1. Rod climbing: flow description
and boundary conditions.
angular velocity a. Along the axis, we impose the classical symmetry boundary conditions, while the kinematic equation (6) is applied along the stress-free upper surface. The flow domain is closed, so that no boundary conditions are required for the viscoelastic stress tensor Ti. In the present section, we focus on the interaction between gravity and normal stresses. Therefore, we take Re = 0 and Cu + 00. For the finite-element discretization, we use the meshes displayed in Fig. 2. They are characterized by a fine definition in the vicinity of the free surface, and near the rod. This is especially true for WEISS~, where a large number of small quadrilateral elements are located in that region, in order to capture both the small secondary motions which are expected there [lo] and the free surface shape. Our first objective in this section is to assess mesh convergence of the numerical results. Both MIX1 [21] and 4 x 4-SU [22] algorithms are used with WEISS& while the first is used with WEISS~ in view of its lower computational cost. Table 1 displays the number of degrees of freedom for the meshes of Fig. 2 with the corresponding algorithms. Beside the free surface shape, we focus attention on the evolution of two quantities as a function of We. (i) The first is the climbing h, as displayed in Fig. 1; for the sake of generality, we use a non-dimensional climbing h/R. (ii) Numerical results will show the development of two secondary vortices, one of them being located under the bulge. The ratios of intensity of the secondary vortices to the azimuthal flow rate are also of interest. In Figs. 3(a)-(d), we plot both the climbing h and the vortex intensities as a function of We for 99 = 0.01. These curves are obtained by means of the three strategies described above and indicate a fairly good mesh
B. Debbaut and B. Hocq /.J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
Fig. 2. Rod climbing: finite-element
meshes.
convergence, even in a range of We numbers where h differs from a quadratic dependence with respect to We. The graphs in Fig. 3 also show that both h and vortex intensities are increasing functions of We. Furthermore, one sees in Fig. 3(b) that for low We, the climbing h is a quadratic function of We as already noted [8-131 and matches the analytically predicted values. This emphasizes the reliability of our numerical predictions. In Refs. 8-13 we find an analytical expression for the free surface shape, which is valid for low angular velocity R. In the particular case of the
TABLE 1 Degrees of freedom for the three numerical strategies Field
WEISSI,
MIX1
Tl u P
10146 5073 446
39486 5073 446
18126 9063 786
Total
15665
45005
27975
WEISSI,
4 x 4-su
WEISS2,
MIX1
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
0.12
h
x
:
111
a
0.3
ox
8
q
0.12
@O
* . 0.4
. a
We0
W&
0.0 0
5
10
15
a
20 x 10-3
x 10”
0
b
0.2
0.4
0.6
0.8
1
x 10”
x 10-4
5 vortex
intensity
4
x B d
3
dp on
2
a”
q
. 1
9 ox
1
.
n+ We
0 ,.m= 0.0
0 0.04
0.08
C
0.12
0.0
0.04
0.06
0.12
d
Fig. 3. (a), (b) Climbing h as a function of We*; (c) intensity of the large vortex as a function of We, for 9 = 0.01; (d) intensity of the small vortex as a function of We, for 6 = 0.01. (+), MIX1 algorithm with WEISSI; (X 1,4X4-W algorithm with WEISSI; (o), MIX1 algorithm with w&s2; ( -1, analytical prediction.
creeping flow of an upper-convected 8-13 reduces to h(r) -= R
Maxwell fluid, the expression in Refs.
2We2 ‘B( r/R)4 ’
(11)
in non-dimensional form. In Fig. 4 we have plotted the numerical prediction of the free surface shape for We = 0.005, 0.01 and 0.015, together with the analytical counterpart given by (11). Excellent agreement is obtained; for low We, the free surface indeed shows a shape following l/r4. In Fig. 5 we display the free surface shape and streamlines of the secondary motions in the neighbourhood of the bulge obtained with WEIW, at various values of We for 9 = 0.01. The growth of a small vortex under
112
.05
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 1
h(r) R t,
*04 .-a We = .015
.a3
.’
‘,
r*
.02
.Ol
.OO 1.0
1.8 r 2.0 R Fig. 4. Free surface shape obtained at various values of We for 6 = 0.01. (+I, MIX1 algorithm with WEISSI; (X), 4 X4-SU algorithm with WEISSI; (01, MIX1 algorithm with ), analytical prediction. WEISSZ; (---
We = 0.025
We = 0.0665
1.2
1.4
1.6
We = 0.0644
We = 0.099
Fig. 5. Free surface shape and streamlines obtained at various values of We for 9 = 0.01. A fluid particle moves counterclockwise along a positive streamline, and clockwise along a negative streamline.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
113
the bulge is indeed observed. As no inertia is involved in the present simulation, the two vortices are generated and driven solely by viscoelasticity. In the large vortex, fluid particles move counterclockwise, while they move clockwise in the small vortex under the bulge and thus climb the rod as already sketched by Yoo et al. [lo]. The separating line between the two opposed secondary motions starts near to a point where the free surface changes curvature. Predicting the location of this point with accuracy seems to be difficult, as it would require an extremely refined mesh. In particular, it is the re-entrant point of the opposite boundary streamlines of both secondary vortices. In the vicinity of this point, at high We, rapid slope changes in the free surface can occur within an extremely short distance as also sketched by Yoo et al. [lo] and experimentally shown in Refs. 8, 14 and 20. At high We, the free surface experiences large deformations and also shows a loss of smoothness in its finite-element discretization, which pollutes the numerical solutions and leads to unrealistic numerical predictions. If a Galerkin technique is used for integrating eqn. (61, spurious oscillations in the free surface appear much earlier. Despite this unrealistic loss of smoothness in the free surface, numerical results can still be obtained. In view of their lack of interest, they are not displayed here. Also, we were not interested in capturing bifurcation points that we believe are numerical artefacts in view of the spurious oscillations. Several possible causes can explain this loss of smoothness which has been observed with all our strategies. Surface tension has been neglected in the present computations, which would have led to a reduction of the changes in curvature. Thus we have investigated here a more critical situation than that with surface tension. Next, qualitatively speaking, pictures in Refs. 8 and 20 showed similar experimental bulge shapes characterized by a rapid slope change over a short distance at high rod angular velocity, while a time-periodic instability was evidenced at higher angular velocity. A transient algorithm combined with an adaptive mesh technique (in the region of rapid change in the free surface) should be used. We must point out here that the streamlines displayed in Fig. 5 are not fluid trajectories; rather they are traces of trajectories in a cross-sectional plane containing the axis of symmetry. Along those traces, a fluid particle moves extremely slowly, and the major flow remains the swirling one. The graphs of Figs. 3(c) and Cd) showing the vortex intensities emphasize this statement. In Fig. 6, we have plotted the contour lines of the azimuthal velocity component ug and the pressure field for We = 0.0644 and B = 0.01. We observe in particular that the azimuthal velocity component is slightly higher in the bulge than along the rod. We believe that it originates from a
114
B. Debbaut and B. Hocq/J.
@ E
Non-Newtonian Fluid Mech. 43 (1992) 103-126
rl
,005
u, ==
I
&i
1.
.Ol
I ,i
b
Fig. 6. Contour lines of the azimuthal velocity component (b), for We = 0.0644 and S? = 0.01.
ug (a>, and of the pressure field p
quasi-rigid motion experienced by the fluid in a region where the shear stress is virtually zero. As far as the pressure field is concerned, its contour lines reveal the expected results. The free surface is an isobar. Far from the rod, the pressure contour lines are controlled by gravity. Close to the rod, they shift from their hydrostatic counterpart and move upwards. This arises from the azimuthal extra-stress component (or normal stress) which rapidly increases towards the rod, as stated in Refs. 2, 3 and 17. If one follows a horizontal line towards the rod, the pressure increases. This is illustrated in Fig. 7, where we give the pressure profile along the radial line passing through the free surface at a long distance from the rod, at We = 0.0644 and 9 = 0.01.
x 1o-3
12
8
5 4 2 3 1 Fig. 7. Pressure profile along the radial line passing through the free surface at long distance from the rod, at We = 0.0644, for 22 = 0.01.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
4. The rod climbing of a Johnson-Segalman
115
liquid
Beavers and Joseph [8] reported theoretical predictions on rod climbing for a fourth-order Rivlin-Ericksen fluid. For low values of a, their predictions showed a good quantitative agreement with experimental measurements done with STP (a 26.6% polyisobutylene solution in petroleum oil). Our purpose here is to check the validity of finite-element predictions for such a flow. In Ref. 8, the diameter of the cylindrical vessel is 30.5 cm; the height of fluid at rest is 7.7 cm. For the rotating rod, we have selected a radius of 0.635 cm. Material data for STP are estimated or given (at a temperature of 25.5”C): p =
0.89 g cmp3
q = 146 P,
’
(1 - 4)h = 0.0162 s,
5 = 0.275, I = 30.9 dyn cm-‘, g = 981 cm ss2.
(12) In view of the non-vanishing second normal stress difference shown by the STP, a Johnson-Segalman model is selected; for 4 in (3) and (12), we have taken the value l/9. From (121, it follows that B = 0.0693. In their experiments, Beavers and Joseph [S] coated the rod with Scotchgard, so that at the contact point the free surface is normal to the rotating rod. In the present simulation, we impose that the free surface has a vanishing slope on the rod. For low angular velocity a, when normal stress differences, inertia forces and surface tension are taken into account, the climbing h can be predicted analytically as follows [g-13], in a dimensional form:
(13)
P = (1 - 4)rlA - 2(1- 4b7~5
(14)
is the so-called climbing constant and is a linear combination of the first and second normal stress coefficients. Equation (13) predicts a quadratic dependence of the climbing h upon R, and thus upon We.
116
B. Debbaut and B. Hocq /J. Non-Newtonian Fluid Mech. 43 (1992) 103 - 126
Fig. 8. Rod climbing: finite-element
mesh.
For the numerical simulation of the present flow as experimentally investigated [8], we use the finite-element mesh shown in Fig. 8 together with the 4 x 4-SU finite-element technique [22]. This mesh is characterized by 30438,3915 and 345 degrees of freedom, respectively for T,, v and p. A fine discretization is specified in the vicinity of the contact point of the rotating rod and the free surface. Such a fine mesh has been selected in order to obtain a good description of the flow patterns. The flow region located below the rod (still considered in Fig. 2) has been removed as in the experimental device described in Ref. 8. To avoid a discontinuity in the azimuthal velocity component at the left bottom point in Fig. 8, vanishing tangential forces are imposed over the last element. Numerical simulations have been performed up to We = 0.360. Figure 9 displays the climbing h/R as a function of We2. For low We, we find that our numerical predictions compare well with the experimental measure-
.6
h
We2
.oo
.04
a3
.12
.16
Fig. 9. Climbing h as a function of We*. (predictions; (01, experimental data.
), analytical prediction;
( + 1, numerical
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
117
y .15
.lO
.05
1.0
1.4
1.8
2.2
Fig. 10. Free surface shape obtained numerical prediction.
R at We = 0.195. (o), experimental
data [8]; (-),
ments [S], and with the analytical predictions from eqn. (13). In particular, we note that h is quadratic in We, as was shown experimentally [S]. At higher We, however, a slight discrepancy arises among these results. The uncertainty on some material parameters in (12), and the selection of a Johnson-Segalman model with one single relaxation time and a purely viscous stress component, are possible causes for such a discrepancy. In Fig. 10, we have plotted the free surface profile obtained at We = 0.195. A comparison with the experimental measurements carried out by Beavers and Joseph [8] is also shown. Again, a slight discrepancy is observed, similar to that mentioned above. A complex development of flow patterns is observed as a function of We. This is illustrated in Fig. 11, where we have displayed the streamlines in a cross-sectional plane containing the axis of symmetry for various values of We. We believe that such a complex development of secondary motions results from the opposing interaction of normal stress differences and inertia. At low We, we essentially observe a large vortex driven by inertia with two centres, which is more intense than the small vortices located on the bottom of the vessel and under the free surface, which are driven by the normal stresses. The vortex intensities have been normalized with respect to the azimuthal flow rate; one finds that the secondary motions are extremely weak as compared to the primary flow. At low We, the direction of the fluid particles clearly indicates the physical origin of these vortices: the fluid particles move clockwise in a vortex driven by inertia and are identified by negative values of the stream function, while they move
118
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
5.0 IO*
l:o 105
We=.045
We= .109
-i.oIO-5
11010-S
We= .142
! We= i
4.0
10-J
.195
-l.o 1o-5
5.0 lD_5
A.5 10”
We= 240
We = .298
-0.25 10-4
.A.5 104
i
/ I
! We= .322
We = .360
Fig. 11. Streamlines obtained at various values of We. A fluid particle moves counterclockwise along a positive streamline and clockwise along a negative streamline.
counterclockwise in a vortex driven by viscoelasticity, where they are identified by positive values of the stream function. These viscoelastic vortices grow with We; this is especially true for the one located under the free surface. This last vortex eventually becomes
B. Debbaut and B. Hocq/J.
Non-Newtonian Fluid Mech. 43 (1992) 103-126
119
more intense than the original inertial vortex. Meanwhile, a third small viscoelastic vortex is generated under the bulge, the intensity of which is very low. We note here that the two viscoelastic vortices under the free surface show a behaviour similar to that seen in the previous section. At higher We, the original inertial vortex is split into two parts. Both these parts are progressively squeezed by the larger viscoelastic vortex, which results from the assembly of the first two viscoelastic vortices and which still grows. Independent of this phenomenon, the third viscoelastic vortex located under the bulge grows slowly. Calculations have not been pursued beyond We = 0.360, as convergence becomes difficult to achieve. Moreover, as previously shown in Fig. 9, a discrepancy occurs between numerical results and experimental measurements of the climbing at high We. These problems are not yet fully understood. A possible cause can be found in the free surface, which shows at least two stagnation points for high We. 5. The Quel1efik.t
Another axisymmetric swirling flow of interest is the so-called QueZZeffekt in Refs. 4 and 18. Figure 12 illustrates this flow. A cylindrical vessel, the radius of which is R, is filled with a fluid up to filling height H. As a disk, located at the bottom of the container, rotates around its axis, the fluid is set into motion and flows upwards along this axis, like a source. A bulge is then generated at the upper surface of the liquid. This bulge can be characterized by its vertical displacement h along the axis of symmetry. For simulating this flow, we have selected the boundary conditions as
wall
:
Fig. 12. Quelleffekt flow description
and boundary conditions.
120
B. Debbaut and B. Hocq/
1
I
I
I
Fig. 13. Quelleffekt: finite-element
I
J. Non-Newtonian
Fluid Mech. 43 (1992) 103-126
I lttl
mesh.
follows. The fluid sticks to the walls of the container and to the rotating disk; symmetry boundary conditions are imposed along the axis, while eqn. (6) is applied on the upper surface which is stress free. To avoid possible difficulties arising from a discontinuity in the azimuthal velocity boundary condition at the lower right-hand side corner of the domain, vanishing tangential forces are imposed on the bottom within a distance of O.OlR from the wall. As for the rod-climbing simulations of the two previous sections, no boundary conditions are specified for the viscoelastic stress tensor T,, since the flow domain is closed. Since our investigations in Section 3 assess the reliability of the results, we have performed our computations with the finite-element mesh of Fig. 13. In particular, it is characterized by 15006, 1953 and 176 degrees of freedom for T,, v and p respectively. It also contains relatively small elements along the boundary. In this section, we again use the finite-element technique developed by Marchal and Crochet [22]. Biihme et al. [18] reported numerical and experimental results for the Quelleffekt of a 2.5% aqueous polyacrylamide solution. In their experimental device, the radius of the vessel was 4.5 cm. Fluid material data are also found in Ref. 18: p =
1.005 g cme3
~=76P (1 -@A
= 0.776 s,
%$ = 0.22, g = 981 cm se2.
(15)
In view of the non-vanishing second normal stress difference, we have selected a Johnson-Segalman model; for 4, we have taken the value l/9. From (15), it follows that 9 = 5.09. Also, in order to observe the effects of the sole first normal stress difference, we also performed computations with an Oldroyd-B model for which we selected 4 = 0.1; therefore % = 5.03.
B. Debbaut and B. Hocq / J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
121
model, at We = 0.323, and H = 1.
Fig. 14. Streamlines obtained with the Johnson-Segalman
For the experiments discussed in Ref. 18, the height of filling H of the fluid in the cylindrical vessel was an independent parameter. For H ranging from 0.5 to 1.25, numerical predictions and experimental measurements carried out by Biihme et al. 1181quantitatively coincide. Figure 14 shows the contour lines obtained with the Johnson-Segalman model at We = 0.323, for H = 1. Our predictions show a very good quantitative agreement with that obtained in Ref. 18. In Fig. 15, we plot the displacement h of the upper surface along the axis of symmetry as a function of H for both Johnson-Segalman and -3
-3
x 10
x 10
16
16
8
8
4
4
0
0 0.5
a
0.9
1.3
0.5
1.7
b
0.9
1.3
#
1.7
Fig. 1.5. Displacement h of the upper surface along the axis of symmetry as a function of H. predictions of Ref. 18; (+), (a) Oldroyd-B model; (b) Johnson-Segalman model. ( -_), present numerical predictions.
122
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103-126 x
10
-3
2.5
8
2.0
6
1.5 4 1.0 2
0.5 0.0
0 0.0
a
0.2
0.4
0.6
0.0
0.04
0.08
0.12
b
Fig. 16. Graph of h as a function of We *, for H = 0.5. (a) Johnson-Segalman model; (b) Oldroyd-B model. ( + ), Numerical predictions; (), predictions of Ref. 18.
Oldroyd-B models, at We = 0.323 and 0.319, respectively. We also display the predictions of Ref. 18, which showed good agreement with experimental measurements since both first and second normal stress differences were taken into account. Our predictions match well those of Ref. 18. One observes in Fig. 15 that h decreases as H increases, and thus as the effects of the rotating disk move away. Moreover, we find that selecting the Oldroyd-B model leads to a larger displacement of the upper surface along the axis of symmetry, as was predicted [18]. This results from the vanishing second normal stress difference of this model. Figure 16 displays the evolution of h as a function of We2 for both Johnson-Segalman and Oldroyd-B models for H = 0.5. Comparison with results of Ref. 18 is also displayed, and again shows excellent agreement. One observes here that for low values of We, h is a quadratic function of We. Again one finds that selecting the Oldroyd-B model, and thus removing the second normal stress difference, leads to a larger deformation of the upper free surface. At high We, h is no longer quadratic in We, as would be expected, since the predictions of Ref. 18 are only valid for low We. In Fig. 17, we display the free surface profiles for various values of H, obtained with both Johnson-Segalman and Oldroyd-B models, at We = 0.323 and 0.319 respectively. These profiles show good agreement with the predictions of Ref. 18. However, for low H, they are characterized by a slight loss of smoothness near the wall of the vessel. The finite-element discretization is probably not sufficient in this region to capture the abrupt deformation experienced by the free surface. In Fig. 18, we plot the intensity of the secondary motion as a function of We, which has been normalized with respect to the azimuthal flow rate.
B. Debbaut and B. Hocq/ .I. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 0.02
0.02
7
y
y 0.01
0.01
0.00
0.00
-0.01
-0.01
-0.02
a
123
-0.02 0.0
0.2
0.4
0.6
0.6
1.0
b
o.0
0.2
0.4
0.6
0.6
1.0
Fig. 17. Free surface shape for various values of H. (a) Johnson-Segalman We = 0.323; (b) Oldroyd-B model, at We = 0.319.
model, at
Again, we find that the secondary motion is extremely weak as compared to the primary one, so that a fluid particle must do a large amount of rotation around the axis of the vessel before recovering its initial location in the (r, z) plane. 6. Conclusions A finite-element technique has been used for simulating axisymmetric swirling flows of viscoelastic liquids of the differential type. Convergence of the results has been proved on the basis of different meshes and numerical algorithms.
0.12 Vortex
intensity
a
0.08
0.04
We
0.00 0.0
0.2
0.4
0.6
0.8
Fig. 18. Vortex intensity as a function of We. (a) Johnson-Segalman model.
model; (b) Oldroyd-B
124
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126
Our technique has been successfully applied to both rod-climbing flow and the Quelleffekt. Our present results on rod climbing and surface displacement agree well with previous predictions. Moreover, we have been also able to quantify the secondary motion intensity for both these flows and to show that the secondary motions are extremely weak as compared to the main azimuthal flow. For the rod-climbing flow, a full simulation of a realistic experiment has been performed, including both inertia terms and surface tension. The development of complex flow patterns has also been shown, as the rod angular velocity increases. In particular, we have successively found that two viscoelastic vortices grow, split the large original inertial vortex into two parts, and are eventually combined into a larger viscoelastic vortex which progressively squeezes the remaining parts of the split inertial vortex. Numerical simulations performed for the Quelleffekt have shown the opposed effects of the first and second normal stress differences on the surface displacement. We have shown that, as far as bulge displacement is concerned, the second normal stress difference acts against the first one. Further flow problems will be investigated in later work. In particular, it is meaningful to perform transient computations which could possibly predict flow instabilities at high Weissenberg numbers. Extension to differential viscoelastic models with multiple relaxation times will also be attempted. Acknowledgement The authors wish to acknowledge the collaboration of Dr. J. Rosenberg from Polyflow S.A. who improved the final version of the English text. References 1 F.H. Garner and A.H. Nissan, Rheological properties of high viscosity solutions of long molecules, Nature, 158 (1946) 634-635. 2 K. Weissenberg, A continuum theory of rheological phenomena, Nature, 159 (1947) 310-311. 3 H.A. Barnes, J.F. Hutton and K. Walters, An Introduction to Rheology, Elsevier, Amsterdam, 1989, Chapter 5. 4 R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids, Vol. 1, John Wiley, New York, 2nd edn., 1987, pp. 62-65. 5 K. Walters, Rheometry, Chapman and Hall, London, 1975. 6 J. Serrin, Poiseuille and Couette flow of non-Newtonian fluids, Z. Angew. Math. Mech., 39 (1959) 295-299. 7 H. Giesekus, Einige Bemerkungen zum FlieBverhalten elasto-viskoser Fliissigkeiten in stationaren Schichtstromungen, Rheol. Acta, 1 (1961) 404-413.
B. Debbaut and B. Hocq/ J. Non-Newtonian Fluid Mech. 43 (1992) 103- 126 8 G.S. Beavers and D.D. Joseph, The rotating rod viscometer,
9
10 11 12 13 14
15 16 17 18 19
20 21 22 23 24 25
26 27 28 29
125
J. Fluid Mech., 69 (1975) 475-511. D.D. Joseph and G.S. Beavers, Free surfaces induced by the motion of viscoelastic fluids, in R.S. Rivlin (Ed.), The Mechanics of Viscoelastic Liquids, AMD 22, ASME, New York, 1977, pp. 59-99. Y. Yoo, D.D. Joseph and G.S. Beavers, Higher-order theory of the Weissenberg effect, J. Fluid Mech., 92 (1979) 529-590. G.S. Beavers, J.Y. Yoo and D.D. Joseph, The free surface on a liquid between cylinders rotating at different speeds. Part III, Rheol. Acta, 19 (1980) 19-31. D.D. Joseph, Fluid Dynamics of Polymeric Liquids, Springer-Verlag, New York, 1990, Chapter 17. D.D. Joseph, G.S. Beavers, A. Cers, C. Dewald, A. Hoger and P.T. Than, Climbing constants for various liquids. J. Rheology, 28 (1984) 325-345. H.H. Hu, 0. Riccius, K.P. Chen, M. Arney and D.D. Joseph, Climbing constant, second-order correction of Trouton’s viscosity, wave speed and delayed die swell for Ml, J. Non-Newtonian Fluid Mech., 35 (1990) 287-307. D.D. Joseph, Remarks on inertial radii, persistent normal stresses, secondary motions and non-elastic extensional viscosities, J. Non-Newtonian Fluid Mech., 32 (1989) 107-114. B. Debbaut and M.J. Crochet, Extensional effects in complex flows, J. Non-Newtonian Fluid Mech., 30 (1988) 169-184. C.T. Hill, Nearly viscometric flow of viscoelastic fluids in the disk and cylinder system. II. Experimental, Trans. Sot. Rheol., 16 (1972) 213-245. G. BBhme, R. VoR and W. Warnecke, Die freie Oberflache einer Fliissigkeit iiber einer rotierende Scheibe, Rheol. Acta, 24 (1985) 22-33. A. Delgado, S. Berg, R. Kroger and H.J. Rath, Theoretical considerations on the steady rod climbing at vanishing influence of gravity, in D.R. Oliver (Ed.), Proc. of the Golden Jubilee of the British Society of Rheology and 3rd European Rheology Congress, Edinburgh, Elsevier, 1990, pp. 129-131. D.D. Joseph, Stability of Fluid Motion II, Springer tracts in Natural Philosophy, Berlin, 1976, pp. 232-237. M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of non-Newtonian Flow, Elsevier, Amsterdam, 1984. J.M. Marchal and M.J. Crochet, A new mixed finite element for calculating viscoelastic flow, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. B. Debbaut, J.M. Marchal and M.J. Crochet, Numerical simulation of highly viscoelastic flows through an abrupt contraction, J. Non-Newtonian Fluid Mech., 29 (1988) 119-146. M.J. Crochet, Numerical simulation of viscoelastic flow: a review, Rubber Chem. Technol., 62(3) (1989) 425-455. R. Keunings, Simulation of viscoelastic fluid flow, in C.L. Tucker III (Ed.), Fundamentals of Computer Modelling for Polymer Processing, Carl Hanser Verlag, Munchen, 1989, pp. 403-469. S. DuPont and M.J. Crochet, Swirling flows of viscoelastic fluids of the integral type in rheogoniometers, Chem. Eng. Commun., 53 (1987) 199-221. J.G. Oldroyd, On the formulation of rheological equation of state, Proc. R. Sot., London, Ser. A, 200 (1950) 523-541. M.W. Johnson, Jr., and D. Segalman, A model for viscoelastic fluid behavior which allows non-affine deformation, J. Non-Newtonian Fluid Mech., 2 (1977) 255-270. X. Georges et B. Hocq, Simulation numerique de l’effect Weissenberg, Memoire, Unite de Mecanique Appliquee, Louvain-la-Neuve, 1986.
126
B. Debbaut and B. Hocq/J.
30 0. Wambersie,
Non-Newtonian Fluid Mech. 43 (1992) 103-126
Simulation numerique de i’extrusion tridimensionnelle par une methode pseudo-transitoire basCe sur les gradients conjugues, Ph. D. Thesis, Unite de Mecanique Appliquee, Louvain-la-Neuve, 1990. 31 R. Keunings, Progress and challenges in computational rheology, Rheol. Acta, 29 (1990) 556-570. 32 S.F. Kistler and L.E. Striven, Coating flows, in J.R.A. Pearson and SM. Richardson (Eds.), Computational Analysis of Polymer Processing, Applied Science Publishers, London, 1983, pp. 243-299. 33 A.N. Brooks and T.J.R. Hughes, Streamline-upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible NavierStokes equations, Comput. Methods Appl. Mech. Eng., 32 (1982) 199-2.59.