- Email: [email protected]
On the die-swell of an axisymmetric Newtonian jet
Computers and b'luids Vol. 8, pp. 275-289 © Pergamon Press Ltd., 1980. Printed in Great Britain
0045-7930180]0901-02751502.00]0
ON THE DIE-SWELL OF AN AXISYMMETRIC NEWTONIAN JET BERNARDJ. OMODEI Basser Department of Computer Science, Universityof Sydney,Sydney,N.S.W. 2006,Australia (Received September 1978;in revisedform 20 February 1979) Abstract--This paper investigates the behaviour of an axisymmetric Newtonian jet of incompressible viscous fluid. Numerical solutions are obtained via the Galerkin formulation of the finite element method for values of the Reynoldsnumber between 0.0001 and 1000,and for various values of the surface tension parameter. Attention is focussed on the shape of the free surface profile and the die-swell ratio. The numerical solutions are subjected to tests for validity based upon the momentumbalance and mass balance equations. The satisfactionof stress boundary conditionsis also tested. The numericalresults are compared with experimental results.
Re S ,f Xp f s n 81 82 z TN tN ts
NOMENCLATURE Reynoldsnumber surface tension parameter numericalsolution for the die-swell ratio predicteddie-swell ratio from formula (5.4) velocityflux which deviates most from unity shear stress on the centreline which deviates most from zero normedstress on CD in Fig. 1 which deviates most from zero maximumvariation in ff(-L0 divided by ff(-L0 maximumvariation in 153(L2)divided by r53(L2) z-componentof the midside and vertex nodes on the free surface prescribednormal stress on free surface numericalsolution for the normal stress on the free surface numericalsolution for the shear stress on the free surface
1. INTRODUCTION In Omodei[9], the steady flow of a plane Newtonian jet of incompressible fluid is investigated, attention being focussed on the shape of the free surface profile and the die-swell ratio. In this paper, we present analogous numerical results for the axisymmetric Newtonian capillary jet. The axisymmetric case obviously has considerable industrial application (for example, in extrusion processes), whereas the plane case is of more theoretical interest. Since the mathematical formulation of the problem and the finite element program used are described in Omodei [9], we shall avoid unnecessary duplication of the material contained therein. Before continuing, let us define the die-swell ratio. Let ho be the radius of the cylindrical tube out of which the liquid is extruded and let ht be the final radius of the jet when the velocity field has become uniform (gravitational forces are ignored); then the die-swell ratio X = hi/ho. A numerical method which accurately predicts the die-swell ratio for a Newtonian jet would provide a firm foundation for extension to non-Newtonian jets where the die-swell ratio is of considerable industrial importance. In Section 2 the mathematical formulation of the problem is given and the Galerkin method is introduced as the basic numerical method for finding an approximate solution. The numerical results of this paper were obtained by using a finite element program called "PROG". The origin and some of the details of " P R O G " are given in Omodei[9], and so we shall redescribe only some of the basic features in Section 3. In Section 4 the numerical results are presented for values of the Reynolds number ranging from 10-4 to 1000, and for various choices of the surface tension parameter. For reasons given in[9], the numerical results obtained should be viewed with a certain degreq of scepticism. In Section 5 the results of Section 4 are subjected to a series of tests for validity, and the outcome of these tests proved favourable. A number of experimenters have obtained experimental results for the die-swell of a Newtonian capillary jet for various values of the Reynolds number and the surface tension parameter. In particular, there are the papers of Batchelor and Horsfall[1], Middleman and Gavis [7], and Goren and Wronski[2]. The results of Goren and Wronski[2] give the shape of the free surface and the die-swell ratio to a CAF Vol. 8, No. 3--A
275
276
B.J. OMODEI
higher level of precision than the results of the other experimenters, and therefore, a comparison is made between the numerical results and the experimental results of[2]. 2. M A T H E M A T I C A L F O R M U L A T I O N
For a steady incompressible axisymmetric Newtonian jet, the governing system of nonlinear elliptic partial differential equations in the interior of the region 1"~defined by Fig. 1 is given by tii:i - pvivi;i = 0,
i = 1, 3
(2.1) (2.2)
vj: i = O,
where i , j = 1,2,3.
t q = - p S q + p.(vi;~+ vj:i),
(2.3)
Note that the usual tensor notion is being used with j = l, 2, 3 corresponding to the cylindrical polar coordinates r, 0, z, respectively, and ;j represents covariant differentiation (see, for example, Huilgol[5]). tii are the components of the stress tensor, p is the fluid density, vi are the components of the velocity vector (vt = u, v3 = w), p is the pressure, jz is the coefficient of viscosity, and 8ii are the components of the unit tensor. Since v2 = 0 and the operator (a/d0) = 0, (2.1)-(2.2) reduce to three differential equations in the three unknowns v~, v3 and p, and there is no dependence on 0. it has been assumed that there is no gravitational force. In imposing boundary conditions, the following assumptions have been made: (i) there is zero pressure in the inviscid fluid surrounding the jet, (ii) fully developed Poiseuille flow exists at a distance L~ upstream of the exit plane, and (iii) there is a uniform velocity field at a distance L2 downstream of the exit plane. The boundary conditions that are prescribed and the notation used are illustrated in Fig. 1, where TN is the imposed normal stress in the direction of the outward normal to the boundary 00 of fl and Ts is the imposed shear stress on c~II, h(z) is the r-coordinate of the free surface, ~r is the coefficient of surface tension, and J is the mean curvature of the free surface and is defined by J(z)
=
h"(z) 1 [1 Jr h'(z)2] 312 h(z)[1 + h'(z)2] 112" -
(2.4)
-
There is no loss of generality in assuming that the tube radius ho = 1, that w = 2(1-/.2) at z= - Lt (this corresponds to a flow rate of unity), and that the viscosity/~ = I. Hence it can be shown that the Reynolds number Re = 2p. A Galerkin approximate solution of (2.1)--(2.3) is sought in the form vi = v'f + ~ , aired~ira,
i = 1, 3,
(2.5)
rPl=l
r
Ts -0, T~-o-J
'-I
D
~,---- u=O,w=O F///////////////////, k Pre~rlbed
'~
velocity
k - 4 u.o
|
~___~ w . 2 . 0 ( i - r ~)
/ LO
u=O
h, ,
h(z)
l
1 h I'
---u-O, r s - 0
L~
~
Z
L2
Fig. 1. Boundary conditions and notation.
H
.
On the die-swellof an axisymmetricNewtonianjet
277
M
/~ = ~
bmrm,
(2.6)
m=l
and this leads to the equations
/ ~j;A~,.r dr dz = O,
m = 1, 2 . . . . . M,
f ( ~j4~,.;j + p~s~,:,4',.)r dr dz = f~a T,qs,,r ds,
(2.7) (2.8)
i = 1, 3 (not summed), m = 1, 2 . . . . . Mj, where T~ are the stress components determined from the prescribed stress boundary conditions. See Omodei[9] for the details. 3. THE FINITE ELEMENT PROGRAM The numerical results for this paper were obtained using a finite element program called "PROG" which has been developed from a collection of subroutines called "AXFINR". A detailed description of "AXFINR" may be found in[8, 10--13]. (a) The nonlinear inertia term pfi~;s in (2.8) in linearized in the (n + 1)th iteration by replacing t~j by v~") which is the velocity obtained as the solution to the nth iteration. (b) The region ~ in Fig. 1 is partitioned into quadrilaterals as shown in Fig. 3, and each quadrilateral is partitioned into four 15 degree-of-freedom triangles. For each triangle, the 15 nodal parameters are defined at vertices or midside nodes as shown in Fig. 2. The trial functions ¢h,,, ¢h,, and ~b,~in (2.5)-(2.6) and the functions v$, v7 in (2.5) are defined in an analogous fashion to those in[9]. Let the unknowns urn(1 -< m -< M~), w,,(1 -< m -< M3) and pro(1-
(3.1)
where the matrix M would be absent if the inertia term were absent from (2.1). K is symmetric and M is unsymmetric (see Nickell et al. [8]). (c) We have not required of the boundary conditions in Fig. 1 that there be zero velocity flux across the proposed free surface. Hence, in order to ensure that there is zero velocity flux across the free surface, we iterate in the following fashion. Let h(")(z) be the approximate free surface for the nth iteration where h(")(z) consists of linear segments as shown in Fig. 3. Let zj, 1 -< j -< N, denote the fixed z-coordinates of the vertex nodes on the free surface with zl = 0 and zN = Lz. The updated free surface h("+~)(z) is defined by linear interpolation through the points (zj, h("+1)(zs)), 1 -< j -< N, where l.(.+l)l.
h("+~)(zj) = ,
~. fzj t~(.).
t,j_ls~- I
j = 2 . . . . . N.
-~-~oz,
*-" Z./_ I
{~,w3,p3.)
(%,
)
{u,,~,#r)
-4 (u4,~)
~, (ua,wz,s~ )
Fig. 2. Triangle nodal parameters.
(3.2)
278
B.J. OMODEI Tube wall
(3 5, I0)
Lip
(0,0, I0)
Exit plane
(-3.5,00)
Centreline Free surface
Lip
(0.0,0.0)
]
\
Exit plane
(00,00) I I 1 ~
(35,00)
~
~
\
~CentrelJne \ surface
~
-~(3.5,0.0)
Iee Centreline
(17.0,0.0
Fig. 3. The mesh used for Re ~ 36.
with ht"+l)(0) = 1, and (t2c"), ~,t,)) is the velocity on ht")(z) and is determined from the solution of eqns (3.1) at the nth iteration. Note that the iteration for the free surface and the iteration for the nonlinear inertia term are performed concurrently. (d) Let us now describe the manner in which surface tension forces are incorporated as surface traction boundary conditions. The motivation for this procedure is in Omodei[9]. Let sC")(z), 0 < z < L2, be the natural cubic spline (see Greville [3]) which satisfies the interpolatory conditions sC")(zj) = ht")(zj),
] = 1, 2 . . . . .
N.
(3.3)
The surface traction vector T(z) on hC")(z) is represented by the continuous piecewise linear function of z which interpolates the vectors crJ(zj)fij, 1 <- j < N, where hi is the outward unit normal to the surface s t") at zj, and st")"(zj) 1 J(zi) - [1 + sC")'(zj)]3/2 s<,)(zi)[l + s<,j,(zi)]l/2.
(3.4)
(e) The manner in which the mesh is chosen is described in[9] so we shall not describe it here. (f) The solution of the system of linear algebraic equations (3.1) is obtained by Gaussian elimination without pivoting. The entries in K and M are approximated by a 7-point Gaussian quadrature scheme over the relevant triangular elements. The details and observations in Omodei[9] apply equally to the solution of eqns (3.1) for the axisymmetric case, and consequently the reader is referred to the important comments made in [9]. 4. NUMERICAL RESULTS In this section we present some of the numerical results that have been obtained for values of the Reynolds number ranging from 10-4 to 103, and for various values of the surface tension parameter which is defined by S = crpho/~ 2.
(4.1)
Let us make a few remarks about the generation of the results. (i) The strategy for choosing a starting approximation (xt°), h °)) in (3.1) is the same as that described in Omodei[9].
On the die-swellof an axisymmetricNewtonianjet
279
(ii) Convergence of the iterative process (3.1) combined with iteration for the free surface was determined in terms of the behaviour of the free surface h ~"). Two tolerances ~ and ¢2 were chosen and the free surface was deemed to have converged when max 2~j<-N
(IhC"+~)(zj) - ht'°(zj)J/zj)
< el
(4.2)
and Ih~"+~)(L2)- h¢"J(L2)l< ~2.
(4.3)
For the results of this paper, ~ and ¢2 were given the values 0.0005 and 0.0002, respectively. For a given mesh configuration, convergence of the free surface h ¢~J according to (4.2) and (4.3) implicitly tests for the convergence of the iterative process (3.1) associated with nonlinear inertia term since the free surface cannot settle down while x¢"J remains significantly different from x~"-~). For the range of values of Re and S considered in this paper, any tendency for the iterations to diverge can be attributed to the free surface configuration rather than the nonlinear inertia term. This conclusion follows from the fact that rapid convergence of the iterative process (3.1) was always observed whenever a proposed free surface was kept fixed for all iterations. (iii) As the ratio S/Re increases, the iterative process (3.1) soon becomes divergent. In order to regain convergence for a large range of values of S/Re, the iteration for the free surface described in Section 3(c) was modified so that the change in the free surface position is not too large; in particular, h t"+~) in (3.2) is replaced by h~ n+~) such that h~, +l) = h2,) + to(h ¢"+l)_ h2-))
(4.4)
where to is called the "damping factor". Relevant factors in the choice of to are described in Omodei[9]. It was found that a good guideline for the choice of to is given by the formula to = min [1.0, 2RdS].
(4.5)
(iv) For values of Reynolds number satisfying Re <-36, it was found that the mesh illustrated in Fig. 3 was completely satisfactory with L~=-3.5 and L2 = 17. However for Re > 36, it was necessary to go farther downstream to ensure that a uniform velocity field was attained. It was found that for 36 < Re-< 1000, the value L2 = 300 was satisfactory, where further values of zj were added to the mesh in Fig. 3 at 28, 46, 75, 120, 190 and 300. The validity of the assumed values of LI and L2 is further discussed in Section 5(v). Results. For some of our numerical results, the free surface profiles have been illustrated graphically (see Figs. 4-12) for a range of values of Re and S. Many of the values of Re chosen are the same as those values used in the experiments of Goren and Wronski[2] so that comparisons can be made. It should be noted that we have illustrated the surface profiles sufficiently far downstream that they are within at least 0.001 units of the final die-swell ratio. All the surfaces are natural cubic splines interpolating the surface vertex nodes which are denoted on the graphs by a vertical stroke crossing the curve. It can be observed from the surface profiles that as the ratio S/Re becomes significantly larger than unity, oscillations in the curvature of free surface become an increasingly significant feature of the surface profiles. Since these surface oscillations tend to be in phase with the position of the vertex nodes on the free surface, the oscillations must be regarded as a property of the numerical solution, and one would not expect this feature in the analytical solution. In this paper, we have only presented numerical results for which the oscillations in the curvature of the free surface do not detract grossly from the quality of the numerical solution. More extreme examples of this type of oscillation have been illustrated in the surface profiles in Omodei[9]. Since we are primarily concerned with an investigation of the die-swell phenomenon, Table 1 contains a complete list of the values of the die-swell ratio (denoted by ,~) obtained from the numerical solutions. Note that values of Re and S have been included for
2110
B. J. OMODEI
./~/~
]2~
$=0 $=00002
/ L
/
/
,
1/
0 8 :L
•
'
....
-
/ ,,
$- =0 O00Q4 ,
~
S=000095
.,/~
S=000008
//," ~
",," /
S:O0001
........
--
~
S:O00014
/
//// /i../ ,,///J~/
O6
,04
/'/{, 7
,p,//
i02 f/ Lip
i02
4
I0000
04
,-%8
i2
~
20
24
12
29
36
Z
Fig. 4. Free surface profiles for Re = 0.0001.
which the surface profiles are not illustrated graphically. However, all the values of ~ tabulated are meaningful in the sense that the value of ,~ determined at the second-last vertex node on the free surface is the same to within three decimal places of accuracy as the value of ,~ at the final vertex node at z = L2. In other words, we have rejected those numerical results for which the oscillatory behaviour has prevented the free surface from settling down to a constant radius. Let us make some general remarks about the numerical results that we have obtained. The results obtained for R e = 10.4, neglecting surface tension, are essentially identical to the results obtained for creeping flow where the inertia forces are neglected. Without surface tension, the axisymmetric Newtonian jet exhibits a 13.0% expansion for creeping flow, and there is neither expansion nor contraction when R e = 14.4. For low Reynolds number flow (Re < 10), the effect of surface tension is to reduce the expansion of the jet, and for high Reynolds number flow, surface tension reduces the contraction of the jet. In general, the ratio SIRe gave an indication of the effect that surface tension would have on altering the position of the free surface. Consequently, it was found that for large values of Re, one required equally large values of S to achieve a significant change. Hence for values of R e >- 65, the effect of surface tension can be regarded as negligible in the sense that the value of S required to move the free surface .''''r
....
I ....
I ....
I
. . . .
[
I ....
. . . .
r ....
J
'''f
....
t
'''f
....
r ....
i ....
I ....
i
....
~ . . . .
I
[
. . . .
] ]
1.12
LI2
.9=0 I.K)
HO
.-/'''I~
.9=1.2
~.oe
?:4 /
J04
/
/
,oe
s-3.6
I
I.~
~
1,02
i
1.04
! Llp ~ l,O01r,,,,i
102
! ....
I .... 0.4
L ....
i,,..J 0.8
....
J ,,,~ 1.2
....
E ,..,
....
16
I. . . . . . 2,0
t .......
24
Z
Fig..5. Free surface profilesfor Re = 4.2.
28
i ....
I .... 32
L ....
I,,,` 36
i00
281
On the die-swell of an axisymmetric Newtonian jet
Table l(a). Numerical data for 10-4 -< Re <- 10.1 (see nomenclature for the definition of the symbols)
Re IE-4 IE-4 1E4 IE-4 IE-4 IE-4 IE-4 IE-4 2.0 3.0 4.2 4.2 4.2 4.2 4.2 5.0 6.0 7.2 7.2 7.2 7.2 8.5 10.1 10.1 10.1 10.1
S 0 2E-5 4E-5 6E-5 8E-5 IE-4 1.4E-4 1.8E-4 0 0 0 1.2 2.4 3.6 4.8 0 0 0 3.6 7.2 10.8 0 0 5.05 10.1 15.15
X 1.130 1.107 1.091 1.080 1.070 1.063 1.052 1.044 1.121 1.113 I.102 1.088 1.077 1.068 1.060 1.093 1.082 1.068 1.060 1.053 1.045 1.053 1.037 1.036 1.034 1.032
XP
1.032 1.032 1.031 1.030
f 1.0007 1.0007 0.9995 1.0008 1.0010 1.0011 0.9991 1.0023 1.0005 1.0006 0.9995 1.0003 1.0009 1.0012 1.0019 1.0005 0.9996 0.9996 1.0007 1.0022 1.0014 0.9996 0.9996 0.9994 0.9989 0.9987
s
n
-0.0081 -4E-7 -0.0080 3E-6 -0.0080 9E-6 -0.0079 2E-5 -0.0078 6E-5 -0,0078 IE-4 -0.0077 4E-4 -0.0076 7E-4 -0.0070 -6E-7 -0.0063 -7E-7 0.0060 -9E-7 0.0061 IE-7 0.0061 3E-5 0.0062 8E-5 0.0062 2E-4 0.0057 -IE-6 0.0057 -IE-6 0.0057 -2E-6 0.0058 2E-5 0.0059 9E-5 0.0059 4E-4 0.0057 -2E-6 0.0054 -2E-6 0.0055 IE-5 0.0056 8E-5 0.0057 3E-4
81
82
9E-6 9E-6 9E-6 9E-6 9E-6 9E-6 9E-6 8E-6 5E-6 4E-6 2E-6 2E-6 2E-6 2E-6 2E-6 3E-6 4E-6 6E-6 6E-6 6E-6 5E-6 8E-6 9E-6 9E-6 9E-6 8E-6
2E-7 3E-6 7E-6 2E-5 4E-5 8E-5 2E-4 4E-4 3E-7 4E-7 4E-7 IE-5 3E-5 9E-5 2E-4 4E-7 5E-7 6E-7 3E-5 IE-4 6E-4 7E-7 8E-7 2E-5 I E-4 4E-4
significantly would be unrealizable in a laboratory situation. It is interesting to note that as Re is increased starting from creeping flow, one reaches a critical value of Re slightly greater than 10 for which the initial effect of introducing surface tension is to increase the die-swell ratio even though the jet exhibits an expansion greater than 3% without surface tension. For a range of values of Reynolds number greater than this critical value, the initial effect of surface tension is to increase the die-swell ratio but as S is increased further, the die-swell ratio reaches a Table l(b). Numerical data for 11.8 -< Re < 1000 (see nomenclature for the definition of symbols)
Re 11.8 ! 1.8 !1.8 !1.8 14.0 14.0 14.0 14.0 14.0 14.8 16.6 16.6 18.0 21.0 24.0 24.0 29.0 36.0 47.4 47.4 65.0 95.0 150.0 250.0 500.0 1000.0
S
X
Xp
0 5.9 11.8 17.7 0 7.0 14.0 21.0 28.0 0 0 16.6 0 0 0 24.0 0 0 0 47.4 0 0 0 0 0 0
1.021 !.022 !.023 1.023 1.003 1.007 i.009 1.010 1.011 0.997 0.987 0.994 0.978 0.964 0.954 0.964 0.940 0.926 0.913 0.922 0.901 0.890 0.881 0.876 0.871 0.869
1.018 1.020 1.021 1.022 1.002 1.006 1.008 1.010 1.011 0.997 0.987 0.994 0.979 0.965 0.955 0.965 0:941 0.928 0.914 0.922 0.902 0.891 0.882 0.876 0.871 0.869
/ 0.9994 0.9992 0.9987 0.9986 0.9993 0.9990 0.9986 0.9984 0.9983 0.9993 0.9996 0.9991 0.9995 0.9991 0.9991 0.9990 0.9991 0.9990 0.9990 1.0009 0.9991 0.9995 0.9993 !.0004 0.9996 !.0005
s
n
81
82
0.0049 0.0049 0.0051 0.0052 -0.0051 -0.0051 -0.0051 -0.0052 -0.0052 -0.0051 -0.0051 -0.0051 -0.0052 -0.0055 -0.0054 -0.0054 -0.0056 -0.0055 -0.0058 -0.0058 -0.0070 -0.0064 -0.0067 -0.0080 -0.0(08 -0.0075
-3E-6 4E-6 5E-5 2E-4 -3E-6 -3E-6 3E-5 2E-4 3E-4 -4E-6 -4E-6 7E-6 -4E-6 -5E-6 -6E-6 -1E-5 -1E-5 4E-5 1E-9 6E-9 3E-10 IE-9 3E-10 IE-9 IE-7 5E-6
IE-5 IE-5 IE-5 9E-6 IE-5 IE-5 IE-5 9E-6 9E-6 IE-5 9E-6 8E-6 7E-6 7E-6 IE-6 9E-6 2E-5 3E-5 4E-5 4E-5 5E-5 5E-5 4E-5 2E-5 2E-5 1E-5
1E-6 IE-5 8E-5 3E-4 1E-6 2E-6 4E-5 2E-4 4E-4 1E-6 IE-6 2E-5 1E-6 IE-6 lE-6 8E-6 4E-6 IE-5 IE-9 3E-8 IE-9 IE-9 IE-9 IE-8 5E-7 2E-4
B. J. OMODE[
282
i
112
'''~
....
I ....
I ....
I ....
t ....
I ....
~
....
I ....
r ....
[ ....
q ....
I ....
' ....
I ....
1 ....
I ....
'
.... I'''~ I 112
F IIO~-
~ I.IO
r ,08~-
7108
:-
S:O
/f
,o,~-
-~
,_s°36
i,o~
,
104
~104
-~102
O0
08
[ .... 3.2
2.4
1.6
i ....
I .... 4.0
, ....
t .... 4.8
i ....
I .... 56
L ....
I .... 6 4
, ....
I,.,,ilO 72
0
Z
Fig. 6. Free surfaceprofilesfor Re = 7.2. maximum value and then starts to decrease. This behaviour is related to the fact that for Reynolds number in the approximate range 7 -< Re <- 24 and without surface tension, the free jet profile immediately after extrusion is concave up but becomes concave down farther downstream; the point of inflection in the surface profile moves from the exit lip to far downstream as Re increases within the range. For Reynolds number in the approximate range 10 -< Re <- 24 and S = 0, the radius of the jet decreases immediately after extrusion, a minimum value of the radius is attained, and then the radius increases as one looks farther downstream. Mesh refinement: By taking a finer mesh, the spatial oscillations in the curvature of the free surface increase in amplitude. Since the surface tension boundary conditions are derived from the curvature, it is not surprising that refining the mesh makes it more difficult to achieve convergence of the free surface. 5. VALIDITY OF NUMERICAL RESULTS In this section we present evidence in support of the statement that the numerical results of Section 4 are indeed good approximations to the solution of the system of partial differential equations (2.1)-(2.3) with the given boundary conditions, for the given choices of Reynolds number Re and surface tension parameter S. We shall also compare the numerical results generated by "PROG" with the experimental results of Goren and Wronski[2]. Let us firstly .'"'1
....
I ....
I ....
1 ....
I
.
.
.
.
I ....
~ ....
I ....
i ....
I ....
' ....
S=O S :5.05
1.04
I
-
'1
'
I ....
r ....
I ....
r ....
I"'~
4
'~
1.04
~s=Jo.i ~ ~ ~ ,
1.03
1,o 3
S=
5.05
1.02
S
0
1.01
1,02
r
1.01
lip
1.00
1.00
0.99
0.99
0.96
..,1 .... 0.0
I .... 0.8
I
....
I ....
1.6
J
....
I .... 2.4
J ....
I .... 3.2
I ....
I .... 4.0
I ....
i .... 4.8
i ....
J .... 5.6
Z
Fig. 7, Free surfaceprofilesfor Re = 10.1.
I ....
] ....
6.4
I...
J ....
7.2
0.98
283
On the die-swell of an axisymmetric Newtonian jet
1.04
1.03 ~-
S=0
1.03
-
S=O 1.02 .
~
.02
J
1.01
1.0[ lip 1.00
1.00
0.99
0.99
0.98 0.0
D.98
0.8
1.6
2.4
3.2
40
48
5.6
64
72
Z
Fig. 8. Free surface profiles for Re • , .... , ~ .... l.... ~ ....l .... ~ .... l .... ~ l .
.
.
]
.
.
.
.
.
l
.
.
~
.
.
.
.
l .
.
~ .
11.8.
=
.
.
l
.
.
.
.
.
~
.
.
.
.
l
.
.
.
J .
.
.
.
.
.
.
.
I
.
.
.
.
.
1 . 0 4
1.04
1,0:5
1.03
1.02
1.02
1.01
-
-
1,01
S lip 1.00
1.00
0.99
0.99
0.98
,,p
....
i ....
i ....
0.8
0.0
i ....
i ....
1.6
l ....
t ....
J ....
2.4
J ....
i ....
3.2
J ....
4.0
i ....
i ....
4.8
i ....
l ....
56
i ....
i ....
6.4
1,
.
0.98
7.2
Z
Fig. 9. Free surface profiles for Re
, , . . i
....
i ....
i ....
i
" ,
•
I ....
i..
l
....
, ....
i ....
lip 1.00
,
" 1
....
14.0.
=
, ....
r ....
i ....
i " '
, ....
i " " .
1.00
Re =16.6, S=16.6 p
~
*
,
Re= 16.6p S=O
, 0.98
0.98 Re:24~
S=24
,
O. 96
0.96
O. 94
0.94
O. 92
~
Re=47.4, S = 0
0.92
,
0.90
0.90
.=
0.88 0.0
....
J ....
0.8
[ ....
I ....
1,6
] ....
J ....
2.4
I ....
I ....
L ....
3,2
I ....
4.0
, ....
~ ....
4.8
* ....
I ....
5.6
= ....
I ....
64
Z
Fig. 10. Free surface profiles for Re = 16.6, 24.0 and 4?.4.
I ....
I...
7.2
0.88
B.J. OMODEI
284
•
1.00 i
I
1.00
098~k
-~0.90
096L
096
o,9z ~
~ 0.92
1o.o 0.88
I,...L.~.~ 0.0 20
.~...=~..~..~
4.0
60
......... l~,J.... : .... , .... J , . , . . . . 1.... , .... ~.... , .... ~.,, ,! 0 . 8 8 8.0 10.0 120 14.0 16.0 18.0 Z
Fig. 1 I. Free surface profiles without surface tension for Re = 65, 95, 150 and 250.
subject the numerical results of Section 4 to a series of tests; however, it should be noted at the outset that the success of these tests does not formally establish the validity of the numerical results. (i) Analytical results. There are no known analytical solutions to the free jet problem. However, there is one analytical result which can be used to test the validity of the numerical results for high Reynolds number and no surface tension. This is the result attributed to Harmon [4] which says that the limit of the die-swell ratio as the Reynolds number approaches infinity is given by x® =
V(3)/2
=
0.8660.
(5.1)
Observe from Table 1 that the numerical die-swell ratio ~ approaches the theoretical limit (5.1) as the Reynolds number becomes very large. (ii) Momentum balance. The test which we shall now discuss is based upon the momentum balance equation (2.1) in the axial direction, and it is in the author's opinion the most significant test of the validity of the numerical results of Section 4. Using the momentum balance analysis of Joseph [6], one can establish the following formula
lip I. 0 0
I00
0.96
0.96 Re=1000 0.92
0.92
0.88
,Re=500
-~
2~
T-
~
Re=500
0.88
Re= 1 0 0 0
0866
i
0.84
0.84
0.80
0.80
0.76
0.76 I0
0
Fig.
12.
20
30
40
50
60
70
80
Free surface profiles without surface tension for Re = 500 and 1000.
90
On the die-swell of an axisymmetric Newtonian jet
for the die-swell ratio in terms of the pressure p ( - L t ) at z = - L~, the shear stress -L~ -< z -< 0, along the capillary wall, and the slope hl of the free surface as z--* 0+:
-~ 2 x3 + [-~ +~ -
+f~ Ts(z) dz LI
~ ] :X - ~P = o• (1 + h~2)"mJ
285
Ts(z),
(5.2)
When cr = 0, (5.2) becomes
x=[~+P(-LO +2 f°L Ts(z)dz]-C'm.
(5.3)
Details of the derivation of (5.2) can be found in Section 50 of Huilgol[5]. Recall that we have assumed, without loss of generality, that the tube radius, the coefficient of viscosity, and the input flow rate are all unity; also recall assumptions (ii) and (iii) in Section 2 in relation to the boundary conditions. The analytical formula (5.2) can be used to test the "consistency" of the numerical results of Section 4 in the following sense: for given values of Re and S, the numerical solution for the flow is used to calculate ~(-LO, f°i.~ Ts(z)dz and h~ and these values are substituted into (5.2) to yield the following formula for the predicted die-swell ratio
XP: ~x,, +
+
+~a-L,"~s(z)dz (1 + ~2)tvz) ,r ,1(~-p--2= O.
(5.4)
When ~r# 0, the value of XP is then determined as the appropriate root of the cubic polynomial (5.4). Xe is now compared with ,f, the numerical die-swell ratio, and explicit values of ,~ and )O, are given in Table 1. See Section 5(v) for further discussion o f / ~ ( - L 0 . The integration of the shear stress term in (5.4) is done by the composite Simpson's rule using the shear stress at the midside as well as the vertex nodes on the tube wall. For low Reynolds number flow, the pressure term and the shear stress term in (5.4) are equally large in magnitude and of opposite sign; consequently, the calculation of XP is dominated by cancellation error, and one cannot expect XP to compare favourably with ,~. For example, when R = 4.2 and S = 1.2, f f ( - L 0 = 14.51 and 2
F Ts(z) LI
dz = 14.71.
Because of the cancellation error, Table 1 does not give values of )¢e for low Reynolds number, namely Re < 10. However, as Re increases, the cancellation effect diminishes and, as can be seen from Table 1, Xp and ;f become very close. If XP and ;~ are in agreement at a certain level of accuracy for a particular Reynolds number with S = 0, then the agreement at that level persists when nonzero values of S are introduced. (iii) Mass balance. The test which we are about to describe is based upon the equation of continuity (2.2). As can be seen from Fig. 3, the mesh defines many lines joining the centreline to either the tube wall or the free surface. Since the input flow rate at z= - L l is assumed to be unity, the velocity flux across any of the conical surfaces generated by the lines mentioned above must be equal to unity if the continuity equation (2.2) is satisfied exactly. For each line, the numerical solution is used to compute the velocity flux associated with it. Let the velocity flux which deviates most from unity be denoted by f, and the values of f for all of the numerical results of Section 4 are given in Table I. This test is always satisfied at a high level of accuracy, and indeed the worst case in Table 1 is when f = 1.0023. As one iterates for the free surface, the value of f approaches unity as the interations converge. Hence, in conjunction with the tests for convergence described in Section 4(ii), one can also use the value of f to test whether the iterations have converged. This approach is particularly useful when the value of the ratio S/Re is becoming large and it is becoming difficult to test for convergence with a high level of confidence (see Section 4(iii) in Omodei[9] for further discussion). (iv) Boundary conditions. Let us investigate the extent to which the boundary conditions illustrated in Fig. 1 are satisfied. Obviously the velocity boundary conditions are satisfied
286
B.J. OMODEI
exactly. Because of the mathematical formulation of the Galerkin method, the imposed stress boundary conditions are only approximately satisfied by the numerical solution. The first stress boundary that we check is the shear stress along the centreline. The shear stress computed from the numerical solution is compared with the imposed value of zero. The shear stress is computed at the vertex and midside nodes along the centreline and the value of the shear stress which deviates most from zero is denoted by s. The values of s for the numerical results of Section 4 are given in Table 1. Similarly, the downstream condition on the normal stress, TN = - tr/ht on CD (see Fig. 1), is checked. Using midside and vertex nodes, let n denote the computed value of the normal stress on CD which deviates most from - o / h r. Values of n are given in Table 1. Finally, and most importantly, let us test the stress boundary conditions associated with the free surface. The imposed surface traction vector T(z) on the piecewise linear surface htn~(z) is defined in Section 3(d). Note that the direction of the normal to the surface at the vertex nodes is given in terms of the interpolatory cubic spline as in Section 3(d), and the direction of the normal to the surface at the midside nodes is obtained by averaging the two unit normal vectors at the adjacent vertex nodes. Let TN = IT[ denote the imposed normal stress, and the associated imposed shear stress equals zero. Let tN and ts denote the normal and shear stress, respectively, on the free surface, where tN and ts are calculated from the numerical solution in the following way. At the midside nodes, tN and ts are determined from the computed velocity and pressure fields in the appropriate triangular element. At each vertex node, tN and ts are obtained by a weighted average of the appropriate stresses in each of the four triangular elements having the vertex node in common. The weights are chosen so that their sum is unity and each is proportional to the angle subtended at the vertex node by the corresponding triangle. The method used for the calculation of tN and ts at midside and vertex nodes is also used in the determination of s and n defined above. Some of the results of this test are given in Table 2 where we have selected a representative sample of values of Re and S. Note that when S = 0, it is not necessary to tabulate TN since TN = 0. Every odd-numbered node in Table 2 is a midside node, and every even-numbered node is a vertex node. The values of (Re, S) represented in Table 2 are (104, 0), (10-4, 0.00006), (7.2, 3.6) and (24, 24). The results of this test are less favourable than those of the other tests in this section; this can be attributed to the tendency of the free surface to exhibit oscillatory curvature, especially when surface tension is significant. It is clear from Table 2 that the test fails in the immediate vicinity of the exit lip at which a stress singularity is known to exist (see Joseph[6]). It is also observed that the stress boundary conditions on the free surface become more closely satisfied as the Reynolds number increases; this can be attributed to the curvature oscillations becoming weaker. (v) Upstream and downstream assumptions. Recall from Section 2 that fully developed Poiseuille flow was assumed to exist at z = - L~. In order to test the validity of this assumption, we check from the numerical solution that the pressure at z= - L~,/~(-L~), is indeed a constant value across the flow from the centreline to the tube wall. For a given choice of Re and S, let 8~ denote the maximum variation in/~(-Li) divided by/~(-L~). The values of/;~ are given in Table 1. Let us now test the validity of the assumption that there is a uniform velocity field at z = L2. Let z33(L2) denote the axial velocity at z -- L2, and let 82 denote the maximum variation in t~3(L2) across the flow divided by z33(L2). The values of/;2 are given in Table 1 for all of the numerical results of Section 4. It is clear that both upstream and downstream assumptions are completely validated by the tabulated values of 8t and 82. Experimental results. Using a glycerol-water mixture, Goren and Wronski[2] performed a series of experiments to record the shape of the free surface of a Newtonian jet for a range of values of the Reynolds number. The surface tension parameter was roughly constant at 0.3. Our aim is to generate numerical solutions for the free jet profile for all the values of Reynolds number recorded graphically in Goren and Wronski[2], the surface tension parameter being kept constant at 0.3. The numerical solutions are then compared with the experimental results. In Fig. 13, we illustrate the numerical solutions for the surface profiles by the continuous lines with Reynolds number ranging from 4.2 to 47.4. The corresponding experimental results for recorded points on the free surfaces are indicated by marked points in Fig. 13 (different symbols are used to distinguish different jet profiles). In reproducing Fig. 2 of Goren and Wronski[2] we believe that the scale that they used for their axial distances in Fig. 2 to be incorrect by a factor
On the die-swell of an axisymmetric Newtonian jet
I
I
287
I
i
I
It . . . .
. . . : ~ . , . ; , , . : ~ . , . : . . . : . . : ~ ~ ~ ~
I
i
I
II
i
i
I
I
I
I I I
i
I
I
I
I
i
I
I
I
I
0
I
II
I I I l l
I
i
II ~
.
.
.
~
.
,
~
.
.
I
11
~
. . . . . . . . . . . .
I
f
I I I
i
I
i
I
f
I 0e-.
.
.
.
.
.
.
.
z
M
I
i
I
I
~
~
~
~
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . .
II
~
I I I I I
I
M M ~ d ~
~
I
I
~
II
.
.
.
I r l
.
I
II
II
I
I
It
i
II
.
.
.
I
.
I l l
.
.
.
.
I
I
I~1
.
I I I I
.
.
.
.
.
I I I I I
.
288
B.J. OMODEI Re=4.2
1.10
Re=7.2 1.0S • • "
x x
x
x
x
-Re=lO.
1
-Re=ll.8 x
1.00
×
× x x
x
*
-Re=14 - - Re= 1 4 . 8 Re=16.6 Re=18
r
0.95
Re=24
Re=47.4
0.90 0.0
1.0
2.0
3.0
4.0
.0
.z
Fig. 13. Free surface profiles for S = 0.3 and a range of values of the Reynolds number; experimental results denoted by either "x" or ".", and numericalresults denoted by continuous curves. of precisely one half (this may be a consequence of using the diameter of the capillary instead of the radius).t After making this change of scale, we found that our numerical solutions were in very close agreement with the experimental results with respect to the shape of the free surfaces and the values of the die-swell ratio. It is interesting to note that Goren and Wronski[2] found that the Reynolds number for which the final radius is the same as the initial jet radius is approximately 14.4, and this is the same value as our numerical solutions predict. Using a coarse mesh, Tanner and Reddy[13] have used " A X F I N R " to generate a numerical solution for the free Newtonian jet with R e =4.2 and S = 0.3. This solution is in close agreement with the corresponding solution that we have generated using " P R O G " . Again using a coarse mesh, Reddy and Tanner[10] have illustrated graphically the shape of the free surface for the following values of the parameters (Re, S): (0,0), (14.8, 0.3), (20, 0), (20, 1), (47.3, 0.3), (50, 0) and (50, 1). Given the limitations on accuracy associated with the coarse mesh, the results of[10] are in general agreement with the corresponding results presented in this paper.
Acknowledgements--The research for this paper began while the author was an A.R.G.C.Research Fellow in the Schoolof
Mathematical Sciences at the Fiinders University of South Australia, and he would like to acknowledge the financial assistance of the AustralianResearch Grants Committee.
REFERENCES 1. J. Batchelorand F. Horsfall, Die-swellin elastic and viscous fluids. Rubber and Plastics Res. Assoc. of Great Britain, Res. Rep. 189 (1971). 2. S. L. Goren and S. Wronski,The shape of low-speedcapillary jets of Newtonian liquids. J. Fluid Mech. 25, 185-198 (1~). 3. T. N. E. Greville (Editor), Theory and Application of Spline Functions. AcademicPress, New York (1969). 4. D. B. Harmon, Drop sizes from low speed jets. J. Franklin Inst. 259, 519 (1955). 5. R. R. Huiigol,Continuum Mechanics of Viscoelastic Fluids. John Wiley, New York (1975). 6. D. D. Joseph, Slow motion and viscometricmotion; stability and bifurcationof the rest state of a simple fluid. Arch. Rat. Mech. Anal. 56, 99--157(1974). 7. S. Middlemanand J. Gavis, Expansionand contraction of capillaryjets of viscoelasticliquids.Phys. Fluids 4, 963-969 (1961). 8. R. E. Nickell, R. I. Tanner and B. Caswell, The solution of viscous incompressiblejet and free-surface flows using finite element methods. 1. Fluid Mech. 65, 189-206(1974). 9. B. J. Omodei,Computer solutionsof a plane Newtonianjet with surface tension. Computers and Fluids 7, 79--96(1979). 10. K. R. Reddy and R. I. Tanner, Finite element solutionof viscous jet flows with surface tension. Computers and Fluids 6, 83-91 (1978). tPrivate communicationwith Goren has not yet establisheddefinitivelythat such an error was made.
On the die-swell of an axisymmetric Newtonian jet
289
11. R. I. Tanner, Finite element computing methods in fluid mechanics and rheology, operating manual for AXFINR program. Brown University (1975). 12. R. I. Tanner, R. E. Nickell and R. W. Bilger, Finite element methods for the solution of some incompressible non-Newtonian fluid mechanics problems with free surfaces. Comp. Meth. Appl. Mech. Eng. 6, 155-174 (1975). 13. R. I. Tanner and K. R. Reddy, Finite element method for the solution of some non-Newtonian fluid mechanics problems; in Finite Element Methods in Engineering (Edited by Y. K. Cheung and S. G. Hutton). University of Adelaide (1976).
0045-7930180]0901-02751502.00]0
ON THE DIE-SWELL OF AN AXISYMMETRIC NEWTONIAN JET BERNARDJ. OMODEI Basser Department of Computer Science, Universityof Sydney,Sydney,N.S.W. 2006,Australia (Received September 1978;in revisedform 20 February 1979) Abstract--This paper investigates the behaviour of an axisymmetric Newtonian jet of incompressible viscous fluid. Numerical solutions are obtained via the Galerkin formulation of the finite element method for values of the Reynoldsnumber between 0.0001 and 1000,and for various values of the surface tension parameter. Attention is focussed on the shape of the free surface profile and the die-swell ratio. The numerical solutions are subjected to tests for validity based upon the momentumbalance and mass balance equations. The satisfactionof stress boundary conditionsis also tested. The numericalresults are compared with experimental results.
Re S ,f Xp f s n 81 82 z TN tN ts
NOMENCLATURE Reynoldsnumber surface tension parameter numericalsolution for the die-swell ratio predicteddie-swell ratio from formula (5.4) velocityflux which deviates most from unity shear stress on the centreline which deviates most from zero normedstress on CD in Fig. 1 which deviates most from zero maximumvariation in ff(-L0 divided by ff(-L0 maximumvariation in 153(L2)divided by r53(L2) z-componentof the midside and vertex nodes on the free surface prescribednormal stress on free surface numericalsolution for the normal stress on the free surface numericalsolution for the shear stress on the free surface
1. INTRODUCTION In Omodei[9], the steady flow of a plane Newtonian jet of incompressible fluid is investigated, attention being focussed on the shape of the free surface profile and the die-swell ratio. In this paper, we present analogous numerical results for the axisymmetric Newtonian capillary jet. The axisymmetric case obviously has considerable industrial application (for example, in extrusion processes), whereas the plane case is of more theoretical interest. Since the mathematical formulation of the problem and the finite element program used are described in Omodei [9], we shall avoid unnecessary duplication of the material contained therein. Before continuing, let us define the die-swell ratio. Let ho be the radius of the cylindrical tube out of which the liquid is extruded and let ht be the final radius of the jet when the velocity field has become uniform (gravitational forces are ignored); then the die-swell ratio X = hi/ho. A numerical method which accurately predicts the die-swell ratio for a Newtonian jet would provide a firm foundation for extension to non-Newtonian jets where the die-swell ratio is of considerable industrial importance. In Section 2 the mathematical formulation of the problem is given and the Galerkin method is introduced as the basic numerical method for finding an approximate solution. The numerical results of this paper were obtained by using a finite element program called "PROG". The origin and some of the details of " P R O G " are given in Omodei[9], and so we shall redescribe only some of the basic features in Section 3. In Section 4 the numerical results are presented for values of the Reynolds number ranging from 10-4 to 1000, and for various choices of the surface tension parameter. For reasons given in[9], the numerical results obtained should be viewed with a certain degreq of scepticism. In Section 5 the results of Section 4 are subjected to a series of tests for validity, and the outcome of these tests proved favourable. A number of experimenters have obtained experimental results for the die-swell of a Newtonian capillary jet for various values of the Reynolds number and the surface tension parameter. In particular, there are the papers of Batchelor and Horsfall[1], Middleman and Gavis [7], and Goren and Wronski[2]. The results of Goren and Wronski[2] give the shape of the free surface and the die-swell ratio to a CAF Vol. 8, No. 3--A
275
276
B.J. OMODEI
higher level of precision than the results of the other experimenters, and therefore, a comparison is made between the numerical results and the experimental results of[2]. 2. M A T H E M A T I C A L F O R M U L A T I O N
For a steady incompressible axisymmetric Newtonian jet, the governing system of nonlinear elliptic partial differential equations in the interior of the region 1"~defined by Fig. 1 is given by tii:i - pvivi;i = 0,
i = 1, 3
(2.1) (2.2)
vj: i = O,
where i , j = 1,2,3.
t q = - p S q + p.(vi;~+ vj:i),
(2.3)
Note that the usual tensor notion is being used with j = l, 2, 3 corresponding to the cylindrical polar coordinates r, 0, z, respectively, and ;j represents covariant differentiation (see, for example, Huilgol[5]). tii are the components of the stress tensor, p is the fluid density, vi are the components of the velocity vector (vt = u, v3 = w), p is the pressure, jz is the coefficient of viscosity, and 8ii are the components of the unit tensor. Since v2 = 0 and the operator (a/d0) = 0, (2.1)-(2.2) reduce to three differential equations in the three unknowns v~, v3 and p, and there is no dependence on 0. it has been assumed that there is no gravitational force. In imposing boundary conditions, the following assumptions have been made: (i) there is zero pressure in the inviscid fluid surrounding the jet, (ii) fully developed Poiseuille flow exists at a distance L~ upstream of the exit plane, and (iii) there is a uniform velocity field at a distance L2 downstream of the exit plane. The boundary conditions that are prescribed and the notation used are illustrated in Fig. 1, where TN is the imposed normal stress in the direction of the outward normal to the boundary 00 of fl and Ts is the imposed shear stress on c~II, h(z) is the r-coordinate of the free surface, ~r is the coefficient of surface tension, and J is the mean curvature of the free surface and is defined by J(z)
=
h"(z) 1 [1 Jr h'(z)2] 312 h(z)[1 + h'(z)2] 112" -
(2.4)
-
There is no loss of generality in assuming that the tube radius ho = 1, that w = 2(1-/.2) at z= - Lt (this corresponds to a flow rate of unity), and that the viscosity/~ = I. Hence it can be shown that the Reynolds number Re = 2p. A Galerkin approximate solution of (2.1)--(2.3) is sought in the form vi = v'f + ~ , aired~ira,
i = 1, 3,
(2.5)
rPl=l
r
Ts -0, T~-o-J
'-I
D
~,---- u=O,w=O F///////////////////, k Pre~rlbed
'~
velocity
k - 4 u.o
|
~___~ w . 2 . 0 ( i - r ~)
/ LO
u=O
h, ,
h(z)
l
1 h I'
---u-O, r s - 0
L~
~
Z
L2
Fig. 1. Boundary conditions and notation.
H
.
On the die-swellof an axisymmetricNewtonianjet
277
M
/~ = ~
bmrm,
(2.6)
m=l
and this leads to the equations
/ ~j;A~,.r dr dz = O,
m = 1, 2 . . . . . M,
f ( ~j4~,.;j + p~s~,:,4',.)r dr dz = f~a T,qs,,r ds,
(2.7) (2.8)
i = 1, 3 (not summed), m = 1, 2 . . . . . Mj, where T~ are the stress components determined from the prescribed stress boundary conditions. See Omodei[9] for the details. 3. THE FINITE ELEMENT PROGRAM The numerical results for this paper were obtained using a finite element program called "PROG" which has been developed from a collection of subroutines called "AXFINR". A detailed description of "AXFINR" may be found in[8, 10--13]. (a) The nonlinear inertia term pfi~;s in (2.8) in linearized in the (n + 1)th iteration by replacing t~j by v~") which is the velocity obtained as the solution to the nth iteration. (b) The region ~ in Fig. 1 is partitioned into quadrilaterals as shown in Fig. 3, and each quadrilateral is partitioned into four 15 degree-of-freedom triangles. For each triangle, the 15 nodal parameters are defined at vertices or midside nodes as shown in Fig. 2. The trial functions ¢h,,, ¢h,, and ~b,~in (2.5)-(2.6) and the functions v$, v7 in (2.5) are defined in an analogous fashion to those in[9]. Let the unknowns urn(1 -< m -< M~), w,,(1 -< m -< M3) and pro(1-
(3.1)
where the matrix M would be absent if the inertia term were absent from (2.1). K is symmetric and M is unsymmetric (see Nickell et al. [8]). (c) We have not required of the boundary conditions in Fig. 1 that there be zero velocity flux across the proposed free surface. Hence, in order to ensure that there is zero velocity flux across the free surface, we iterate in the following fashion. Let h(")(z) be the approximate free surface for the nth iteration where h(")(z) consists of linear segments as shown in Fig. 3. Let zj, 1 -< j -< N, denote the fixed z-coordinates of the vertex nodes on the free surface with zl = 0 and zN = Lz. The updated free surface h("+~)(z) is defined by linear interpolation through the points (zj, h("+1)(zs)), 1 -< j -< N, where l.(.+l)l.
h("+~)(zj) = ,
~. fzj t~(.).
t,j_ls~- I
j = 2 . . . . . N.
-~-~oz,
*-" Z./_ I
{~,w3,p3.)
(%,
)
{u,,~,#r)
-4 (u4,~)
~, (ua,wz,s~ )
Fig. 2. Triangle nodal parameters.
(3.2)
278
B.J. OMODEI Tube wall
(3 5, I0)
Lip
(0,0, I0)
Exit plane
(-3.5,00)
Centreline Free surface
Lip
(0.0,0.0)
]
\
Exit plane
(00,00) I I 1 ~
(35,00)
~
~
\
~CentrelJne \ surface
~
-~(3.5,0.0)
Iee Centreline
(17.0,0.0
Fig. 3. The mesh used for Re ~ 36.
with ht"+l)(0) = 1, and (t2c"), ~,t,)) is the velocity on ht")(z) and is determined from the solution of eqns (3.1) at the nth iteration. Note that the iteration for the free surface and the iteration for the nonlinear inertia term are performed concurrently. (d) Let us now describe the manner in which surface tension forces are incorporated as surface traction boundary conditions. The motivation for this procedure is in Omodei[9]. Let sC")(z), 0 < z < L2, be the natural cubic spline (see Greville [3]) which satisfies the interpolatory conditions sC")(zj) = ht")(zj),
] = 1, 2 . . . . .
N.
(3.3)
The surface traction vector T(z) on hC")(z) is represented by the continuous piecewise linear function of z which interpolates the vectors crJ(zj)fij, 1 <- j < N, where hi is the outward unit normal to the surface s t") at zj, and st")"(zj) 1 J(zi) - [1 + sC")'(zj)]3/2 s<,)(zi)[l + s<,j,(zi)]l/2.
(3.4)
(e) The manner in which the mesh is chosen is described in[9] so we shall not describe it here. (f) The solution of the system of linear algebraic equations (3.1) is obtained by Gaussian elimination without pivoting. The entries in K and M are approximated by a 7-point Gaussian quadrature scheme over the relevant triangular elements. The details and observations in Omodei[9] apply equally to the solution of eqns (3.1) for the axisymmetric case, and consequently the reader is referred to the important comments made in [9]. 4. NUMERICAL RESULTS In this section we present some of the numerical results that have been obtained for values of the Reynolds number ranging from 10-4 to 103, and for various values of the surface tension parameter which is defined by S = crpho/~ 2.
(4.1)
Let us make a few remarks about the generation of the results. (i) The strategy for choosing a starting approximation (xt°), h °)) in (3.1) is the same as that described in Omodei[9].
On the die-swellof an axisymmetricNewtonianjet
279
(ii) Convergence of the iterative process (3.1) combined with iteration for the free surface was determined in terms of the behaviour of the free surface h ~"). Two tolerances ~ and ¢2 were chosen and the free surface was deemed to have converged when max 2~j<-N
(IhC"+~)(zj) - ht'°(zj)J/zj)
< el
(4.2)
and Ih~"+~)(L2)- h¢"J(L2)l< ~2.
(4.3)
For the results of this paper, ~ and ¢2 were given the values 0.0005 and 0.0002, respectively. For a given mesh configuration, convergence of the free surface h ¢~J according to (4.2) and (4.3) implicitly tests for the convergence of the iterative process (3.1) associated with nonlinear inertia term since the free surface cannot settle down while x¢"J remains significantly different from x~"-~). For the range of values of Re and S considered in this paper, any tendency for the iterations to diverge can be attributed to the free surface configuration rather than the nonlinear inertia term. This conclusion follows from the fact that rapid convergence of the iterative process (3.1) was always observed whenever a proposed free surface was kept fixed for all iterations. (iii) As the ratio S/Re increases, the iterative process (3.1) soon becomes divergent. In order to regain convergence for a large range of values of S/Re, the iteration for the free surface described in Section 3(c) was modified so that the change in the free surface position is not too large; in particular, h t"+~) in (3.2) is replaced by h~ n+~) such that h~, +l) = h2,) + to(h ¢"+l)_ h2-))
(4.4)
where to is called the "damping factor". Relevant factors in the choice of to are described in Omodei[9]. It was found that a good guideline for the choice of to is given by the formula to = min [1.0, 2RdS].
(4.5)
(iv) For values of Reynolds number satisfying Re <-36, it was found that the mesh illustrated in Fig. 3 was completely satisfactory with L~=-3.5 and L2 = 17. However for Re > 36, it was necessary to go farther downstream to ensure that a uniform velocity field was attained. It was found that for 36 < Re-< 1000, the value L2 = 300 was satisfactory, where further values of zj were added to the mesh in Fig. 3 at 28, 46, 75, 120, 190 and 300. The validity of the assumed values of LI and L2 is further discussed in Section 5(v). Results. For some of our numerical results, the free surface profiles have been illustrated graphically (see Figs. 4-12) for a range of values of Re and S. Many of the values of Re chosen are the same as those values used in the experiments of Goren and Wronski[2] so that comparisons can be made. It should be noted that we have illustrated the surface profiles sufficiently far downstream that they are within at least 0.001 units of the final die-swell ratio. All the surfaces are natural cubic splines interpolating the surface vertex nodes which are denoted on the graphs by a vertical stroke crossing the curve. It can be observed from the surface profiles that as the ratio S/Re becomes significantly larger than unity, oscillations in the curvature of free surface become an increasingly significant feature of the surface profiles. Since these surface oscillations tend to be in phase with the position of the vertex nodes on the free surface, the oscillations must be regarded as a property of the numerical solution, and one would not expect this feature in the analytical solution. In this paper, we have only presented numerical results for which the oscillations in the curvature of the free surface do not detract grossly from the quality of the numerical solution. More extreme examples of this type of oscillation have been illustrated in the surface profiles in Omodei[9]. Since we are primarily concerned with an investigation of the die-swell phenomenon, Table 1 contains a complete list of the values of the die-swell ratio (denoted by ,~) obtained from the numerical solutions. Note that values of Re and S have been included for
2110
B. J. OMODEI
./~/~
]2~
$=0 $=00002
/ L
/
/
,
1/
0 8 :L
•
'
....
-
/ ,,
$- =0 O00Q4 ,
~
S=000095
.,/~
S=000008
//," ~
",," /
S:O0001
........
--
~
S:O00014
/
//// /i../ ,,///J~/
O6
,04
/'/{, 7
,p,//
i02 f/ Lip
i02
4
I0000
04
,-%8
i2
~
20
24
12
29
36
Z
Fig. 4. Free surface profiles for Re = 0.0001.
which the surface profiles are not illustrated graphically. However, all the values of ~ tabulated are meaningful in the sense that the value of ,~ determined at the second-last vertex node on the free surface is the same to within three decimal places of accuracy as the value of ,~ at the final vertex node at z = L2. In other words, we have rejected those numerical results for which the oscillatory behaviour has prevented the free surface from settling down to a constant radius. Let us make some general remarks about the numerical results that we have obtained. The results obtained for R e = 10.4, neglecting surface tension, are essentially identical to the results obtained for creeping flow where the inertia forces are neglected. Without surface tension, the axisymmetric Newtonian jet exhibits a 13.0% expansion for creeping flow, and there is neither expansion nor contraction when R e = 14.4. For low Reynolds number flow (Re < 10), the effect of surface tension is to reduce the expansion of the jet, and for high Reynolds number flow, surface tension reduces the contraction of the jet. In general, the ratio SIRe gave an indication of the effect that surface tension would have on altering the position of the free surface. Consequently, it was found that for large values of Re, one required equally large values of S to achieve a significant change. Hence for values of R e >- 65, the effect of surface tension can be regarded as negligible in the sense that the value of S required to move the free surface .''''r
....
I ....
I ....
I
. . . .
[
I ....
. . . .
r ....
J
'''f
....
t
'''f
....
r ....
i ....
I ....
i
....
~ . . . .
I
[
. . . .
] ]
1.12
LI2
.9=0 I.K)
HO
.-/'''I~
.9=1.2
~.oe
?:4 /
J04
/
/
,oe
s-3.6
I
I.~
~
1,02
i
1.04
! Llp ~ l,O01r,,,,i
102
! ....
I .... 0.4
L ....
i,,..J 0.8
....
J ,,,~ 1.2
....
E ,..,
....
16
I. . . . . . 2,0
t .......
24
Z
Fig..5. Free surface profilesfor Re = 4.2.
28
i ....
I .... 32
L ....
I,,,` 36
i00
281
On the die-swell of an axisymmetric Newtonian jet
Table l(a). Numerical data for 10-4 -< Re <- 10.1 (see nomenclature for the definition of the symbols)
Re IE-4 IE-4 1E4 IE-4 IE-4 IE-4 IE-4 IE-4 2.0 3.0 4.2 4.2 4.2 4.2 4.2 5.0 6.0 7.2 7.2 7.2 7.2 8.5 10.1 10.1 10.1 10.1
S 0 2E-5 4E-5 6E-5 8E-5 IE-4 1.4E-4 1.8E-4 0 0 0 1.2 2.4 3.6 4.8 0 0 0 3.6 7.2 10.8 0 0 5.05 10.1 15.15
X 1.130 1.107 1.091 1.080 1.070 1.063 1.052 1.044 1.121 1.113 I.102 1.088 1.077 1.068 1.060 1.093 1.082 1.068 1.060 1.053 1.045 1.053 1.037 1.036 1.034 1.032
XP
1.032 1.032 1.031 1.030
f 1.0007 1.0007 0.9995 1.0008 1.0010 1.0011 0.9991 1.0023 1.0005 1.0006 0.9995 1.0003 1.0009 1.0012 1.0019 1.0005 0.9996 0.9996 1.0007 1.0022 1.0014 0.9996 0.9996 0.9994 0.9989 0.9987
s
n
-0.0081 -4E-7 -0.0080 3E-6 -0.0080 9E-6 -0.0079 2E-5 -0.0078 6E-5 -0,0078 IE-4 -0.0077 4E-4 -0.0076 7E-4 -0.0070 -6E-7 -0.0063 -7E-7 0.0060 -9E-7 0.0061 IE-7 0.0061 3E-5 0.0062 8E-5 0.0062 2E-4 0.0057 -IE-6 0.0057 -IE-6 0.0057 -2E-6 0.0058 2E-5 0.0059 9E-5 0.0059 4E-4 0.0057 -2E-6 0.0054 -2E-6 0.0055 IE-5 0.0056 8E-5 0.0057 3E-4
81
82
9E-6 9E-6 9E-6 9E-6 9E-6 9E-6 9E-6 8E-6 5E-6 4E-6 2E-6 2E-6 2E-6 2E-6 2E-6 3E-6 4E-6 6E-6 6E-6 6E-6 5E-6 8E-6 9E-6 9E-6 9E-6 8E-6
2E-7 3E-6 7E-6 2E-5 4E-5 8E-5 2E-4 4E-4 3E-7 4E-7 4E-7 IE-5 3E-5 9E-5 2E-4 4E-7 5E-7 6E-7 3E-5 IE-4 6E-4 7E-7 8E-7 2E-5 I E-4 4E-4
significantly would be unrealizable in a laboratory situation. It is interesting to note that as Re is increased starting from creeping flow, one reaches a critical value of Re slightly greater than 10 for which the initial effect of introducing surface tension is to increase the die-swell ratio even though the jet exhibits an expansion greater than 3% without surface tension. For a range of values of Reynolds number greater than this critical value, the initial effect of surface tension is to increase the die-swell ratio but as S is increased further, the die-swell ratio reaches a Table l(b). Numerical data for 11.8 -< Re < 1000 (see nomenclature for the definition of symbols)
Re 11.8 ! 1.8 !1.8 !1.8 14.0 14.0 14.0 14.0 14.0 14.8 16.6 16.6 18.0 21.0 24.0 24.0 29.0 36.0 47.4 47.4 65.0 95.0 150.0 250.0 500.0 1000.0
S
X
Xp
0 5.9 11.8 17.7 0 7.0 14.0 21.0 28.0 0 0 16.6 0 0 0 24.0 0 0 0 47.4 0 0 0 0 0 0
1.021 !.022 !.023 1.023 1.003 1.007 i.009 1.010 1.011 0.997 0.987 0.994 0.978 0.964 0.954 0.964 0.940 0.926 0.913 0.922 0.901 0.890 0.881 0.876 0.871 0.869
1.018 1.020 1.021 1.022 1.002 1.006 1.008 1.010 1.011 0.997 0.987 0.994 0.979 0.965 0.955 0.965 0:941 0.928 0.914 0.922 0.902 0.891 0.882 0.876 0.871 0.869
/ 0.9994 0.9992 0.9987 0.9986 0.9993 0.9990 0.9986 0.9984 0.9983 0.9993 0.9996 0.9991 0.9995 0.9991 0.9991 0.9990 0.9991 0.9990 0.9990 1.0009 0.9991 0.9995 0.9993 !.0004 0.9996 !.0005
s
n
81
82
0.0049 0.0049 0.0051 0.0052 -0.0051 -0.0051 -0.0051 -0.0052 -0.0052 -0.0051 -0.0051 -0.0051 -0.0052 -0.0055 -0.0054 -0.0054 -0.0056 -0.0055 -0.0058 -0.0058 -0.0070 -0.0064 -0.0067 -0.0080 -0.0(08 -0.0075
-3E-6 4E-6 5E-5 2E-4 -3E-6 -3E-6 3E-5 2E-4 3E-4 -4E-6 -4E-6 7E-6 -4E-6 -5E-6 -6E-6 -1E-5 -1E-5 4E-5 1E-9 6E-9 3E-10 IE-9 3E-10 IE-9 IE-7 5E-6
IE-5 IE-5 IE-5 9E-6 IE-5 IE-5 IE-5 9E-6 9E-6 IE-5 9E-6 8E-6 7E-6 7E-6 IE-6 9E-6 2E-5 3E-5 4E-5 4E-5 5E-5 5E-5 4E-5 2E-5 2E-5 1E-5
1E-6 IE-5 8E-5 3E-4 1E-6 2E-6 4E-5 2E-4 4E-4 1E-6 IE-6 2E-5 1E-6 IE-6 lE-6 8E-6 4E-6 IE-5 IE-9 3E-8 IE-9 IE-9 IE-9 IE-8 5E-7 2E-4
B. J. OMODE[
282
i
112
'''~
....
I ....
I ....
I ....
t ....
I ....
~
....
I ....
r ....
[ ....
q ....
I ....
' ....
I ....
1 ....
I ....
'
.... I'''~ I 112
F IIO~-
~ I.IO
r ,08~-
7108
:-
S:O
/f
,o,~-
-~
,_s°36
i,o~
,
104
~104
-~102
O0
08
[ .... 3.2
2.4
1.6
i ....
I .... 4.0
, ....
t .... 4.8
i ....
I .... 56
L ....
I .... 6 4
, ....
I,.,,ilO 72
0
Z
Fig. 6. Free surfaceprofilesfor Re = 7.2. maximum value and then starts to decrease. This behaviour is related to the fact that for Reynolds number in the approximate range 7 -< Re <- 24 and without surface tension, the free jet profile immediately after extrusion is concave up but becomes concave down farther downstream; the point of inflection in the surface profile moves from the exit lip to far downstream as Re increases within the range. For Reynolds number in the approximate range 10 -< Re <- 24 and S = 0, the radius of the jet decreases immediately after extrusion, a minimum value of the radius is attained, and then the radius increases as one looks farther downstream. Mesh refinement: By taking a finer mesh, the spatial oscillations in the curvature of the free surface increase in amplitude. Since the surface tension boundary conditions are derived from the curvature, it is not surprising that refining the mesh makes it more difficult to achieve convergence of the free surface. 5. VALIDITY OF NUMERICAL RESULTS In this section we present evidence in support of the statement that the numerical results of Section 4 are indeed good approximations to the solution of the system of partial differential equations (2.1)-(2.3) with the given boundary conditions, for the given choices of Reynolds number Re and surface tension parameter S. We shall also compare the numerical results generated by "PROG" with the experimental results of Goren and Wronski[2]. Let us firstly .'"'1
....
I ....
I ....
1 ....
I
.
.
.
.
I ....
~ ....
I ....
i ....
I ....
' ....
S=O S :5.05
1.04
I
-
'1
'
I ....
r ....
I ....
r ....
I"'~
4
'~
1.04
~s=Jo.i ~ ~ ~ ,
1.03
1,o 3
S=
5.05
1.02
S
0
1.01
1,02
r
1.01
lip
1.00
1.00
0.99
0.99
0.96
..,1 .... 0.0
I .... 0.8
I
....
I ....
1.6
J
....
I .... 2.4
J ....
I .... 3.2
I ....
I .... 4.0
I ....
i .... 4.8
i ....
J .... 5.6
Z
Fig. 7, Free surfaceprofilesfor Re = 10.1.
I ....
] ....
6.4
I...
J ....
7.2
0.98
283
On the die-swell of an axisymmetric Newtonian jet
1.04
1.03 ~-
S=0
1.03
-
S=O 1.02 .
~
.02
J
1.01
1.0[ lip 1.00
1.00
0.99
0.99
0.98 0.0
D.98
0.8
1.6
2.4
3.2
40
48
5.6
64
72
Z
Fig. 8. Free surface profiles for Re • , .... , ~ .... l.... ~ ....l .... ~ .... l .... ~ l .
.
.
]
.
.
.
.
.
l
.
.
~
.
.
.
.
l .
.
~ .
11.8.
=
.
.
l
.
.
.
.
.
~
.
.
.
.
l
.
.
.
J .
.
.
.
.
.
.
.
I
.
.
.
.
.
1 . 0 4
1.04
1,0:5
1.03
1.02
1.02
1.01
-
-
1,01
S lip 1.00
1.00
0.99
0.99
0.98
,,p
....
i ....
i ....
0.8
0.0
i ....
i ....
1.6
l ....
t ....
J ....
2.4
J ....
i ....
3.2
J ....
4.0
i ....
i ....
4.8
i ....
l ....
56
i ....
i ....
6.4
1,
.
0.98
7.2
Z
Fig. 9. Free surface profiles for Re
, , . . i
....
i ....
i ....
i
" ,
•
I ....
i..
l
....
, ....
i ....
lip 1.00
,
" 1
....
14.0.
=
, ....
r ....
i ....
i " '
, ....
i " " .
1.00
Re =16.6, S=16.6 p
~
*
,
Re= 16.6p S=O
, 0.98
0.98 Re:24~
S=24
,
O. 96
0.96
O. 94
0.94
O. 92
~
Re=47.4, S = 0
0.92
,
0.90
0.90
.=
0.88 0.0
....
J ....
0.8
[ ....
I ....
1,6
] ....
J ....
2.4
I ....
I ....
L ....
3,2
I ....
4.0
, ....
~ ....
4.8
* ....
I ....
5.6
= ....
I ....
64
Z
Fig. 10. Free surface profiles for Re = 16.6, 24.0 and 4?.4.
I ....
I...
7.2
0.88
B.J. OMODEI
284
•
1.00 i
I
1.00
098~k
-~0.90
096L
096
o,9z ~
~ 0.92
1o.o 0.88
I,...L.~.~ 0.0 20
.~...=~..~..~
4.0
60
......... l~,J.... : .... , .... J , . , . . . . 1.... , .... ~.... , .... ~.,, ,! 0 . 8 8 8.0 10.0 120 14.0 16.0 18.0 Z
Fig. 1 I. Free surface profiles without surface tension for Re = 65, 95, 150 and 250.
subject the numerical results of Section 4 to a series of tests; however, it should be noted at the outset that the success of these tests does not formally establish the validity of the numerical results. (i) Analytical results. There are no known analytical solutions to the free jet problem. However, there is one analytical result which can be used to test the validity of the numerical results for high Reynolds number and no surface tension. This is the result attributed to Harmon [4] which says that the limit of the die-swell ratio as the Reynolds number approaches infinity is given by x® =
V(3)/2
=
0.8660.
(5.1)
Observe from Table 1 that the numerical die-swell ratio ~ approaches the theoretical limit (5.1) as the Reynolds number becomes very large. (ii) Momentum balance. The test which we shall now discuss is based upon the momentum balance equation (2.1) in the axial direction, and it is in the author's opinion the most significant test of the validity of the numerical results of Section 4. Using the momentum balance analysis of Joseph [6], one can establish the following formula
lip I. 0 0
I00
0.96
0.96 Re=1000 0.92
0.92
0.88
,Re=500
-~
2~
T-
~
Re=500
0.88
Re= 1 0 0 0
0866
i
0.84
0.84
0.80
0.80
0.76
0.76 I0
0
Fig.
12.
20
30
40
50
60
70
80
Free surface profiles without surface tension for Re = 500 and 1000.
90
On the die-swell of an axisymmetric Newtonian jet
for the die-swell ratio in terms of the pressure p ( - L t ) at z = - L~, the shear stress -L~ -< z -< 0, along the capillary wall, and the slope hl of the free surface as z--* 0+:
-~ 2 x3 + [-~ +~ -
+f~ Ts(z) dz LI
~ ] :X - ~P = o• (1 + h~2)"mJ
285
Ts(z),
(5.2)
When cr = 0, (5.2) becomes
x=[~+P(-LO +2 f°L Ts(z)dz]-C'm.
(5.3)
Details of the derivation of (5.2) can be found in Section 50 of Huilgol[5]. Recall that we have assumed, without loss of generality, that the tube radius, the coefficient of viscosity, and the input flow rate are all unity; also recall assumptions (ii) and (iii) in Section 2 in relation to the boundary conditions. The analytical formula (5.2) can be used to test the "consistency" of the numerical results of Section 4 in the following sense: for given values of Re and S, the numerical solution for the flow is used to calculate ~(-LO, f°i.~ Ts(z)dz and h~ and these values are substituted into (5.2) to yield the following formula for the predicted die-swell ratio
XP: ~x,, +
+
+~a-L,"~s(z)dz (1 + ~2)tvz) ,r ,1(~-p--2= O.
(5.4)
When ~r# 0, the value of XP is then determined as the appropriate root of the cubic polynomial (5.4). Xe is now compared with ,f, the numerical die-swell ratio, and explicit values of ,~ and )O, are given in Table 1. See Section 5(v) for further discussion o f / ~ ( - L 0 . The integration of the shear stress term in (5.4) is done by the composite Simpson's rule using the shear stress at the midside as well as the vertex nodes on the tube wall. For low Reynolds number flow, the pressure term and the shear stress term in (5.4) are equally large in magnitude and of opposite sign; consequently, the calculation of XP is dominated by cancellation error, and one cannot expect XP to compare favourably with ,~. For example, when R = 4.2 and S = 1.2, f f ( - L 0 = 14.51 and 2
F Ts(z) LI
dz = 14.71.
Because of the cancellation error, Table 1 does not give values of )¢e for low Reynolds number, namely Re < 10. However, as Re increases, the cancellation effect diminishes and, as can be seen from Table 1, Xp and ;f become very close. If XP and ;~ are in agreement at a certain level of accuracy for a particular Reynolds number with S = 0, then the agreement at that level persists when nonzero values of S are introduced. (iii) Mass balance. The test which we are about to describe is based upon the equation of continuity (2.2). As can be seen from Fig. 3, the mesh defines many lines joining the centreline to either the tube wall or the free surface. Since the input flow rate at z= - L l is assumed to be unity, the velocity flux across any of the conical surfaces generated by the lines mentioned above must be equal to unity if the continuity equation (2.2) is satisfied exactly. For each line, the numerical solution is used to compute the velocity flux associated with it. Let the velocity flux which deviates most from unity be denoted by f, and the values of f for all of the numerical results of Section 4 are given in Table I. This test is always satisfied at a high level of accuracy, and indeed the worst case in Table 1 is when f = 1.0023. As one iterates for the free surface, the value of f approaches unity as the interations converge. Hence, in conjunction with the tests for convergence described in Section 4(ii), one can also use the value of f to test whether the iterations have converged. This approach is particularly useful when the value of the ratio S/Re is becoming large and it is becoming difficult to test for convergence with a high level of confidence (see Section 4(iii) in Omodei[9] for further discussion). (iv) Boundary conditions. Let us investigate the extent to which the boundary conditions illustrated in Fig. 1 are satisfied. Obviously the velocity boundary conditions are satisfied
286
B.J. OMODEI
exactly. Because of the mathematical formulation of the Galerkin method, the imposed stress boundary conditions are only approximately satisfied by the numerical solution. The first stress boundary that we check is the shear stress along the centreline. The shear stress computed from the numerical solution is compared with the imposed value of zero. The shear stress is computed at the vertex and midside nodes along the centreline and the value of the shear stress which deviates most from zero is denoted by s. The values of s for the numerical results of Section 4 are given in Table 1. Similarly, the downstream condition on the normal stress, TN = - tr/ht on CD (see Fig. 1), is checked. Using midside and vertex nodes, let n denote the computed value of the normal stress on CD which deviates most from - o / h r. Values of n are given in Table 1. Finally, and most importantly, let us test the stress boundary conditions associated with the free surface. The imposed surface traction vector T(z) on the piecewise linear surface htn~(z) is defined in Section 3(d). Note that the direction of the normal to the surface at the vertex nodes is given in terms of the interpolatory cubic spline as in Section 3(d), and the direction of the normal to the surface at the midside nodes is obtained by averaging the two unit normal vectors at the adjacent vertex nodes. Let TN = IT[ denote the imposed normal stress, and the associated imposed shear stress equals zero. Let tN and ts denote the normal and shear stress, respectively, on the free surface, where tN and ts are calculated from the numerical solution in the following way. At the midside nodes, tN and ts are determined from the computed velocity and pressure fields in the appropriate triangular element. At each vertex node, tN and ts are obtained by a weighted average of the appropriate stresses in each of the four triangular elements having the vertex node in common. The weights are chosen so that their sum is unity and each is proportional to the angle subtended at the vertex node by the corresponding triangle. The method used for the calculation of tN and ts at midside and vertex nodes is also used in the determination of s and n defined above. Some of the results of this test are given in Table 2 where we have selected a representative sample of values of Re and S. Note that when S = 0, it is not necessary to tabulate TN since TN = 0. Every odd-numbered node in Table 2 is a midside node, and every even-numbered node is a vertex node. The values of (Re, S) represented in Table 2 are (104, 0), (10-4, 0.00006), (7.2, 3.6) and (24, 24). The results of this test are less favourable than those of the other tests in this section; this can be attributed to the tendency of the free surface to exhibit oscillatory curvature, especially when surface tension is significant. It is clear from Table 2 that the test fails in the immediate vicinity of the exit lip at which a stress singularity is known to exist (see Joseph[6]). It is also observed that the stress boundary conditions on the free surface become more closely satisfied as the Reynolds number increases; this can be attributed to the curvature oscillations becoming weaker. (v) Upstream and downstream assumptions. Recall from Section 2 that fully developed Poiseuille flow was assumed to exist at z = - L~. In order to test the validity of this assumption, we check from the numerical solution that the pressure at z= - L~,/~(-L~), is indeed a constant value across the flow from the centreline to the tube wall. For a given choice of Re and S, let 8~ denote the maximum variation in/~(-Li) divided by/~(-L~). The values of/;~ are given in Table 1. Let us now test the validity of the assumption that there is a uniform velocity field at z = L2. Let z33(L2) denote the axial velocity at z -- L2, and let 82 denote the maximum variation in t~3(L2) across the flow divided by z33(L2). The values of/;2 are given in Table 1 for all of the numerical results of Section 4. It is clear that both upstream and downstream assumptions are completely validated by the tabulated values of 8t and 82. Experimental results. Using a glycerol-water mixture, Goren and Wronski[2] performed a series of experiments to record the shape of the free surface of a Newtonian jet for a range of values of the Reynolds number. The surface tension parameter was roughly constant at 0.3. Our aim is to generate numerical solutions for the free jet profile for all the values of Reynolds number recorded graphically in Goren and Wronski[2], the surface tension parameter being kept constant at 0.3. The numerical solutions are then compared with the experimental results. In Fig. 13, we illustrate the numerical solutions for the surface profiles by the continuous lines with Reynolds number ranging from 4.2 to 47.4. The corresponding experimental results for recorded points on the free surfaces are indicated by marked points in Fig. 13 (different symbols are used to distinguish different jet profiles). In reproducing Fig. 2 of Goren and Wronski[2] we believe that the scale that they used for their axial distances in Fig. 2 to be incorrect by a factor
On the die-swell of an axisymmetric Newtonian jet
I
I
287
I
i
I
It . . . .
. . . : ~ . , . ; , , . : ~ . , . : . . . : . . : ~ ~ ~ ~
I
i
I
II
i
i
I
I
I
I I I
i
I
I
I
I
i
I
I
I
I
0
I
II
I I I l l
I
i
II ~
.
.
.
~
.
,
~
.
.
I
11
~
. . . . . . . . . . . .
I
f
I I I
i
I
i
I
f
I 0e-.
.
.
.
.
.
.
.
z
M
I
i
I
I
~
~
~
~
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . .
II
~
I I I I I
I
M M ~ d ~
~
I
I
~
II
.
.
.
I r l
.
I
II
II
I
I
It
i
II
.
.
.
I
.
I l l
.
.
.
.
I
I
I~1
.
I I I I
.
.
.
.
.
I I I I I
.
288
B.J. OMODEI Re=4.2
1.10
Re=7.2 1.0S • • "
x x
x
x
x
-Re=lO.
1
-Re=ll.8 x
1.00
×
× x x
x
*
-Re=14 - - Re= 1 4 . 8 Re=16.6 Re=18
r
0.95
Re=24
Re=47.4
0.90 0.0
1.0
2.0
3.0
4.0
.0
.z
Fig. 13. Free surface profiles for S = 0.3 and a range of values of the Reynolds number; experimental results denoted by either "x" or ".", and numericalresults denoted by continuous curves. of precisely one half (this may be a consequence of using the diameter of the capillary instead of the radius).t After making this change of scale, we found that our numerical solutions were in very close agreement with the experimental results with respect to the shape of the free surfaces and the values of the die-swell ratio. It is interesting to note that Goren and Wronski[2] found that the Reynolds number for which the final radius is the same as the initial jet radius is approximately 14.4, and this is the same value as our numerical solutions predict. Using a coarse mesh, Tanner and Reddy[13] have used " A X F I N R " to generate a numerical solution for the free Newtonian jet with R e =4.2 and S = 0.3. This solution is in close agreement with the corresponding solution that we have generated using " P R O G " . Again using a coarse mesh, Reddy and Tanner[10] have illustrated graphically the shape of the free surface for the following values of the parameters (Re, S): (0,0), (14.8, 0.3), (20, 0), (20, 1), (47.3, 0.3), (50, 0) and (50, 1). Given the limitations on accuracy associated with the coarse mesh, the results of[10] are in general agreement with the corresponding results presented in this paper.
Acknowledgements--The research for this paper began while the author was an A.R.G.C.Research Fellow in the Schoolof
Mathematical Sciences at the Fiinders University of South Australia, and he would like to acknowledge the financial assistance of the AustralianResearch Grants Committee.
REFERENCES 1. J. Batchelorand F. Horsfall, Die-swellin elastic and viscous fluids. Rubber and Plastics Res. Assoc. of Great Britain, Res. Rep. 189 (1971). 2. S. L. Goren and S. Wronski,The shape of low-speedcapillary jets of Newtonian liquids. J. Fluid Mech. 25, 185-198 (1~). 3. T. N. E. Greville (Editor), Theory and Application of Spline Functions. AcademicPress, New York (1969). 4. D. B. Harmon, Drop sizes from low speed jets. J. Franklin Inst. 259, 519 (1955). 5. R. R. Huiigol,Continuum Mechanics of Viscoelastic Fluids. John Wiley, New York (1975). 6. D. D. Joseph, Slow motion and viscometricmotion; stability and bifurcationof the rest state of a simple fluid. Arch. Rat. Mech. Anal. 56, 99--157(1974). 7. S. Middlemanand J. Gavis, Expansionand contraction of capillaryjets of viscoelasticliquids.Phys. Fluids 4, 963-969 (1961). 8. R. E. Nickell, R. I. Tanner and B. Caswell, The solution of viscous incompressiblejet and free-surface flows using finite element methods. 1. Fluid Mech. 65, 189-206(1974). 9. B. J. Omodei,Computer solutionsof a plane Newtonianjet with surface tension. Computers and Fluids 7, 79--96(1979). 10. K. R. Reddy and R. I. Tanner, Finite element solutionof viscous jet flows with surface tension. Computers and Fluids 6, 83-91 (1978). tPrivate communicationwith Goren has not yet establisheddefinitivelythat such an error was made.
On the die-swell of an axisymmetric Newtonian jet
289
11. R. I. Tanner, Finite element computing methods in fluid mechanics and rheology, operating manual for AXFINR program. Brown University (1975). 12. R. I. Tanner, R. E. Nickell and R. W. Bilger, Finite element methods for the solution of some incompressible non-Newtonian fluid mechanics problems with free surfaces. Comp. Meth. Appl. Mech. Eng. 6, 155-174 (1975). 13. R. I. Tanner and K. R. Reddy, Finite element method for the solution of some non-Newtonian fluid mechanics problems; in Finite Element Methods in Engineering (Edited by Y. K. Cheung and S. G. Hutton). University of Adelaide (1976).