A boundary element investigation of extrudate swell

Journal ojNon - Newtonian Fluid Mechanics, 18 (1985) 143-162 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

A BOUNDARY

ELEMENT

INVESTIGATION

OF EXTRUDATE

143

SWELL

M.B. BUSH *, R.I. TANNER and N. PHAN-THIEN Department (Australia)

of Mechanical

Engineering,

The University of Sydney, New South Wales, 2006

(Received December 3, 1984)

Summary Numerical studies of die swell have until now dealt primarily with the Maxwell or Oldroyd-B viscoelastic models. However, these models exhibit features that often make them unsuitable for numerical work. Furthermore, they are not realistic representations of actual viscoelastic fluids. In this report a comparison is made between the behaviour of a variety of different viscoelastic models when applied to the die swell problem. A wide range of elongational and shear behaviour is exhibited by the models examined. Both types of behaviour are shown to be important in the die swell problem, and the observed swelling is related to these characteristics of the models.

1. Introduction The viscoelastic extrudate swell problem is of practical interest in the polymer processing industry. However, many attempts to solve this problem numerically have been unsuccessful in the sense that practical values of the Weissenberg number and swelling ratio were not reached. Reddy and Tanner [l] considered the plane extrusion of a second order fluid, while most recent investigations have concentrated on the Maxwell fluid [2-81 and Oldroyd-B fluid [9]. A variety of methods have been employed in these investigations to solve the constitutive equations. Finite element schemes using differential constitutive equation forms were the most popular [2-6,8-91, although the integral form of the Maxwell model has also been applied [7]. In the latter case the memory integral is evaluated along the streamlines (particle paths). * Present Address: Department of Mechanical Engineering, The University of Western Australia, Nedlands, Western Australia, 6009, Australia. 0377-0257/85/$03.30

0 1985 Elsevier Science Publishers B.V.

144 More recently the Maxwell constitutive equations have been solved by integration of the differential form of the model along streamlines [8]. This is made possible by the fact that the streamlines represent characteristic curves of the differential operator. A wide range of models can be treated by this method and it is not necessary to first write the model in memory integral form. Both streamline integration approaches have the disadvantage that the streamlines, or particle paths, must be known before the extra stresses can be computed. However, they represent a saving in storage space, program complexity and computation time over conventional finite element formulations. We will consider the streamline integration method further in this report. Unfortunately, the upper convected Maxwell model (or Oldroyd-B model) is not a good representation of real viscoelastic fluid behaviour in the extrudate swell problem. The behaviour is unrealistic in both steady shear flow (the viscosity is independent of shear rate) and in steady elongational flow (the Trouton viscosity is unbounded). Both types of flow are represented in the die swell problem. The relative simplicity of the Maxwell model makes it attractive for use in the development of numerical techniques. However, more realistic models will be required to treat practical extrusion problems. With this goal in mind the present report deals with a comparison between several different models applied to the die swell problem: Oldroyd-B model, Phan-Thien-Tanner model [lo-111, modified Phan-Thien-Tanner model [12], Leonov model [13-141 and modified Leonov model [15]. These models are summarized in Section 2. A wide range of elongational and shear behaviour is represented by this group of models, and we relate this behaviour to the performance of the models when applied to the die swell problem. Our comparisons show that the steady elongational and shear behaviour of the models strongly influence the response of the extruded fluid. A feature of the published results for the Maxwell fluid is the lack of quantitative agreement between the solutions reported by different research groups [8]. This indicates a strong dependence of the numerical solution on the discretization scheme, functional approximations and the entry or exit lengths employed. To avoid similar characteristics interfering with the comparisons reported here, we have employed a single mesh and solution technique throughout. Because of the large number of computations performed in this analysis we elected to employ a relatively coarse grid. Detailed analysis of the most promising models using refined meshes has been reserved for future work. Boundary element methods are used to solve the momentum and continuity equations while the streamline integration technique is used to solve the constitutive equations [8]. Details of the formulations are presented in Section 5.

145 Another phenomenon requiring further study is the question of solution stability. A common feature of all the investigations to date is the presence of a relatively low Weissenberg number limit, above which the iterative solution scheme fails to converge. This value is strongly dependent upon the discretization scheme employed. The reasons for this breakdown are not fully clear. However, it is likely that the failure is caused by a combination of numerical instability, due to characteristics of the discretization and numerical formulation, and inherent instability of the constitutive model itself [16]. Phan-Thien [12,17,18] has illustrated that above a certain Weissenberg number limit the Maxwell model is inherently unstable to small perturbations in a number of simple flows. Such perturbations will be applied during a numerical solution, and it is therefore likely that inherent instability is important. Phan-Thien shows that the Oldroyd-B model is unstable in the same situations, but that the instability sets in at a higher value of the Weissenberg number. Crochet and Keunings [9] have applied the Oldroyd-B model to the extrudate swell problem and were able to achieve solutions at values of the Weissenberg number up to 4.5 (relaxation time x shear rate at the die wall) with the order of 100% swelling. These values are considerably greater than those achieved previously with the Maxwell model, and this behaviour may be attributed, at least in part, to the phenomenon reported by Phan-Thien. We make no attempt to study the stability phenomenon in detail in the present report, except to note the different Weissenberg number limits encountered using different models. 2. Governing equations The flows considered in this report are assumed to be steady, isothermal, creeping and incompressible. Under these conditions the momentum and continuity equations take the form: v*o=Q, v.u=o,

(1)

(2) where (I is the total stress tensor, u is the velocity vector and 0 represents the null vector. The stress tensor in a viscoelastic fluid having a single relaxation mode can be written in the general form: u=

-p6+277,0+7,

(3) where p represents the pressure, S is the unit tensor, II, is the viscosity of a Newtonian solvent, D is the rate-of-strain tensor and T is the extra stress tensor. The constitutive equations governing the extra stress will be of the type:

Xg+R=O,

146 where AT/A? indicates the upper convected derivative of the stress tensor, X is the relaxation time and R depends on the particular model chosen. The upper convected derivative is given by AT &r z=7&+“’

VT-Lr-TLT,

where L is the velocity ered: (i) Maxwell/

gradient

tensor.

The following

models

are consid-

Oldroyd-B model

. In this case we define a constant choose

molecular

contributed

viscosity,

n,,,, and

(6)

R=v2qmD. The Oldroyd-B model arises if a non-zero In this work we have chosen v,,.& = 8.

value of n, is adopted

in eqn. (3).

(ii) Phan- Thien- Tanner model (PTT) The model choosing

proposed

by Phan-Thien

and Tanner

T + Xl{ DT + TO’} - 27~~0.

[lO,ll]

is obtained

by

(7)

In this, q, is again constant (in this case the zero shear rate viscosity of the model) and tr( T) represents the trace of T. The par&meters c and 5 are of the order 0.01 and 0.1 respectively. The elongational behaviour of the model is determined by C, while [ affects the shear behaviour (given that E is small). The values E = 0.01 and I= 0.1 were adopted for use with the PIT model. Unfortunately, the PIT model does not behave realistically in steady shear flow. One finds that the shear component of T approaches zero at high rates of shear. The effect of this can be avoided in an ad hoc manner by ensuring that the term ns in eqn. (3) is non-zero, so that the total shear stress contains a Newtonian component. It is easily shown that if c is of order 0.01 then the total shear stress grows monotonically with the shear rate when n,,,/~ < 8. In the present work we have chosen nm/ns = 8. However, this alteration has a detrimental effect on the behaviour of the recoverable shear, S, (S, = first normal stress difference + shear stress). If the recoverable shear is plotted against shear rate, one finds that a maximum occurs at a finite value of the shear rate. Thus, the relative significance of the first normal stress difference decreases at higher values of the shear rate. A more

147 pleasing solution to this problem was recently suggested The “modified PTT model” is obtained by writing I+e$u(r)

R=

7+X{{&+&}

i

77,(?)

[12].

(8)

-2Tj&)D,

)

where y = /m

by Phan-Thien

and

1+ 3(2 - S)X’$’ = Vo (I + ~2j/2)(w2

*

For given E and 5, the parameter 0 -C n G 1 provides control of the steady shear behaviour. However, the model shows unrealistic behaviour in shear flow (such as shear thickening) unless 5 = 0. The values c = 0.01 and [ = 0 were adopted for use with the modified P’IT model. Throughout the present report the original Phan-Thien-Tanner model (eqn. 7) will be referred to as PTT, while the modified PTT model (eqn. 8) will be represented by MPTT. (iii) Leonov model The Leonov R=1/29

model [ 13,141 is obtained

$rr-21*s+$r

2Ptr(r-‘)-6

t

by choosing tr(7)

, I)

[

where p is a rheological parameter (2Xp = zero shear rate viscosity of the model) and I#I= 1. Although the Leonov model fits steady shear data extremely well, it is incapable of representing elongational flow realistically; the steady Trouton viscosity is virtually independent of strain rate. Larson [15] recently suggested a modification to the model which allows the elongational behaviour to be altered without greatly affecting the shear behaviour. In this the parameter C#Iin (10) is replaced by a function of the form $8= (1 + y

(11)

arctan(O.O5IV/p)]-I,

where W is an “elastic W=1/2~

potential”

&tr(r)+2ptr(r-l)-6). 1

function

given by

(12)

The significance of the function C$is discussed in detail by Larson and, for the purpose of brevity, we will do no more than state the model here. If /I = 0 we recover the original Leonov model. The effect of p on the properties of the model is discussed in the next section.

148 3. Steady shear behaviour The rate-of-strain

tensor in steady simple shear, with shear rate Jo, is given

by (13) The steady shear viscosity is simply qs + ri2/jl, where ri2 is the “extra” shear stress. A measure of the relative importance of the first normal stress difference, N,, in shear flow is the recoverable shear, S,, written as S, = N,/2r,,

.

(14)

In some cases (Maxwell fluid, Oldroyd-B fluid and Leonov fluid with p = 0) it was possible to obtain analytical expressions for the steady shear viscosity and recoverable shear [14,19]. In other cases, it was necessary to solve the governing equations numerically. This was achieved by integrating eqn. (4) through time, after the sudden application of the required steady shearing motion to a previously static body of fluid. A fourth order Runge-Kutta algorithm was employed for this purpose. The integration was continued until steady state conditions were reached. The results are shown in Figs. 1 and 2, where the steady shear viscosity and S, are plotted against the dimensionless shear rate, hq. The viscosity of the Maxwell model and Oldroyd-B model is independent of shear rate. A more realistic shear-thinning behaviour is exhibited by the MPTT, PTT and Leonov models. The parameter /3 has a negligible effect on the viscosity of the Leonov model (eqn. 11). In contrast, a wide range of shear-thinning behaviour can be obtained using the PTT and MPTT models. Note that the steepest curve corresponds to the PTT model with 77,= 0. Indeed, the curve is so steep that the slope of the corresponding shear stress versus shear rate curve becomes zero at a finite shear rate, and becomes negative at higher shear rates. This is corrected by ensuring that qrn/vs < 8; the case of q,/n, = 8 is illustrated in Fig. 1. However, this also imposes a lower limit on the viscosity at high shear rates. The behaviour of the recoverable shear (Fig. 2) is acceptable in all cases except for the PTT model with qs # 0. In the case of the Maxwell model we find that S, = Xq, while S, = q,,,h?/(q, + qs) for the Oldroyd-B model. At low shear rates the recoverable shear corresponding to the PTT model ( qs = 0) is approximateIy S, = A?, and deviates only slightly from this at high shear rates. In the case of the Leonov model, it can be shown that S, is proportional to /m at large values of the shear rate. An analytical expression for all Xi/ can also be found for the case p = 0, by using the exact solution for steady shear flow given in [14]. The effect of /3 is clearly to

149 3

3.4

2

-i _--

1

-

1

Z$O ;g “5 5

-1

% s

-2 -2

-1

0

At,

ht

1

2

(POWERS

3

OF 10)

Fig. 1. Steady shear viscosity and Trouton viscosity as a function of the dimensionless shear rate, A+, or strain rate, Ai: (1) Maxwell (q, = 0, qm = 1, A = 1); (2) Oldroyd-B (11, = 0.112, nm=0.888, X=1); (3) PTT (c=0.01,:~=0.1, qs=o, n,=l, h=l); (4) PTT (c=O.Ol, 5= 0.1, qls = 0.112, nm = 0.888, X =l); (5) MP’IT (n =l, c = 0.01, S = 0, no =l, qS = 0, 1 =l); (6) MPTT (n = 0.5, E = 0.01, S = 0, no = 1, qS = 0, X =l); (7) Leonov (p = 0.5, n, = 0, A= 1, fi = 0); (8) Leonov (c = 0.5, nS = 0, h =l, /!I = 9); - - - - - -, planar elongation of Leonov fluid. The Trouton

viscosity

of the MPTT fluid approaches

2/c

at high strain rates.

3

2 s! IL 0

u-l

1

% 3 8 Y

0

IiY

-1

-2 -2

-1

0

A$

1

(POWERS

Fig. 2. Recoverable shear, S,, see caption to Fig. 1.

2

3

OF 10)

as a function

of the dimensionless

shear rate, A+. For legend,

150 articularly increase S, at a given shear rate. However, the effect is not great; for large X?, we find that S, is increased by a factor of ,!e_ (1 + /3) . We will see later that p has a much greater effect on the elongational flow characteristics of the Leonov model. In the case of the P’IT model with q, # 0, the relative importance of the first normal stress difference diminishes at high shear rates. Thus one might expect that any flow phenomenon attributable to the presence of first normal stress differences will diminish as some measure of the “effective” shear rate (Weissenberg number) is increased. This was found to be so in the die-swell problem. 4. Steady elongational behaviour The rate-of-strain-tensor can be written in the form

00 I*

1 0

D=i [

in steady

0 m

0 0 -(1+m)

elongational

flow, with strain

rate i,

(15

The direction of stretching (with principle rate of strain <) is the xi direction. The parameter “m ” determines the type of elongational flow. The case m = - l/2 produces axisymmetric (uniaxial) elongation, while m = 0 corresponds to planar elongation. The Trouton viscosity, qr, is defined as ?1T= N,/i,

(16)

where Ni = rli - rj3 is the first normal stress difference. An analytical expression for qT in both planar and uniaxial elongation of the Maxwell and Oldroyd-B fluids is well known; see [19]. It appears to be less well known that the Trouton viscosity corresponding to the planar elongation of the Leonov fluid (/3 = 0) can also be found explicitly. The surprising result is that TjT =

8hE.l.

(17)

That is, qT is independent of the strain rate. Thus the Trouton viscosity of the Leonov fluid (/3 = 0) in steady planar elongation behaves exactly like that of a Newtonian fluid. In all other cases it was necessary to evaluate the Trouton viscosity numerically. As before, eqn. (4) was integrated through time, after the sudden application of the required elongational flow. Integration continued until steady state conditions were reached. The results are plotted in Fig. 1, as a function of the dimensionless strain rate Xi. To avoid cluttering the

151 figure, only the case of uniaxial elongation is shown. These results are also characteristic of the planar case. For reference, however, the planar case for the Leonov fluid has been included. The Maxwell and Oldroyd-B models show unsatisfactory behaviour in the sense that vr is unbounded at Xi = l/2. On the other hand, the Leonov model (j3 = 0) is unrealistic in that qr varies only slightly (if at all) from the zero strain rate value. The effect of p is significant. It can be shown that the high strain rate asymptote is increased by a factor of (1 + p). Recall that in steady shear, p has very little effect on the steady shear viscosity and altered S, at high shear rates by a factor of /m. Thus the parameter p provides significant freedom to fit the Leonov model to elongational flow data, without greatly affecting the fit to shear data. This is illustrated by Larson [ 151. On the basis of the foregoing discussion it can be concluded that the Leonov and MPTT models represent the most promising fluids in the sense that they can potentially fit both steady shear and elongational data successfully. We will find that these models also display the most promising behaviour when applied to the die swell problem. 5. Numerical

formulation

Details of the numerical formulation have been presented previously [8,20,21]. Only the important details are repeated here. To begin we must rewrite the total stress tensor (eqn. 3) in the equivalent form: u=

-ps+2?@+e,

(18)

where 271LD represents the total linear part of the stress tensor and e represents the non-linear part. The values adopted for n,_ and e were: Oldroyd-B (Maxwell) and PIT models: q, = ns + q,, e = T - 271~0. MPTT model: ni_ = ns + no, e = r - 27@. Leonov model: vi_ = nI, + 2Xp, e = r - 4ApLD. This ensures that the linear part of the model is non-zero, even if ns = 0. An integral equation representation of eqns. (l), (2) and (18) can be written in Cartesian tensor form:

In this, the domain of the problem is represented by 3, the boundary is represented by P, and P and Q represent arbitrary points in a or on I. The

152 terms ejk( Q) and uj( Q) represent the components of e and u at the point Q. The term tj(Q) represents the components of the boundary traction vector at Q (fj = ujknk, where nk is the outward directed normal vector at the boundary). The weighting field u$( P, Q) represents the velocity field at Q due to a “Stokeslet” directed in the xi direction at P. The Stokeslet is taken to act in an infinite expanse of fluid with viscosity qr,. Likewise, the term t$( P, Q) represents the traction field at Q due to a Stokeslet at P. For details of the Stokeslet, see Ladyzhenskaya [22]. Due to the nature of the singularity contained in the Stokeslet at the point P, the integral involving ti”;.(P, Q) .must be interpreted in the sense of its “Cauchy principle value” when P lies on the boundary. The term Cij( P) accounts for this. If P lies in 52 then Cij( P) = 6,,, the Kronecker delta function. If P lies on l? then Cij( P) depends on the nature of the surface at P. In pacticular, Cij( P) = l/26,, if the surface at P has a defined normal vector and tangent plane (that is, P does not lie on a sharp corner). If the stress distribution is known, eqn. (19) can be discretized and combined with the prescribed boundary conditions to generate a system of simultaneous linear equations which can be solved for the unknown boundary values of velocity and traction. Once these are known, eqn. (19) can be reapplied to compute the velocity components at any point in the domain. In reality the stress distribution is not known a priori and an iterative scheme must be employed in which the velocity field and stress field are alternately estimated. Although the above formulation applies equally to both planar and three dimensional flows, we will restrict discussion to planar flows from this point. Discretization of the integral equation was achieved by replacing the boundary with a series of geometrically linear “boundary elements” and by dividing the domain into an array of triangular “domain cells’. On each boundary element the velocity and traction were taken to be constant. These were defined at a centrally located node. On each domain cell the velocity and extra stresses were defined at the cell verticies (the domain nodes), thereby providing a linear variation across the cell. For a given velocity field in 0, the velocity gradients were evaluated by fitting quadratic polynomials to blocks of nine domain nodes and differentiating the resulting polynomials. The integrals over individual boundary elements and domain cells were evaluated using ordinary Gaussian quadrature rules, except when the singular point P was situated on the element or cell in question. In the case of a boundary element, the singular integral was evaluated exactly. In the case of a domain cell, transformation to polar coordinates centred at the point P cancelled the singularity in the kernel of the integral and allowed ordinary product integration rules to be applied. The upper convected constitutive model (eqn. 4) is a system of first order

153 partial differential equations for the extra stresses. In the case of steady two dimensional and axisymmetric flows it is easily shown that the “characteristic curves” of these equations correspond to streamlines, if the velocity components are known. Thus, for a given steady velocity field, the extra stresses can be evaluated by integrating the constitutive model as a system of ordinary differential equations along the streamlines. Integration will proceed from an upstream point, at which the extra stresses are known, to the downstream point at which the solution is required. If the velocity at the required point is zero, the system of differential equations reduces in general to a system of non-linear simultaneous equations (in the case of the Maxwell and Oldroyd-B models a system of linear equations is produced, which is easily solved). Such a non-linear system usually permits multiple solutions, and there is no guarantee that a numerical solution to these equations will always yield the correct physical solution. An alternative means of solution is to retain the time derivative in eqn. (5) and treat the case u = 0 as a time dependent problem. That is, eqn. (5) was integrated through time after the sudden application of the required velocity gradients to a previously static body of fluid. Integration was continued until steady state conditions were achieved. This procedure can make use of the same integration package used to integrate along streamlines. It was always found to efficiently produce the correct physical solution. To summarize, the total numerical solution was obtained by the following steps: (i) For a given starting stress distribution and free surface geometry (with zero stress boundary conditions prescribed on the free surface), the unknown values of velocity and boundary traction were obtained by solution of eqn. (19). Standard Gaussian elimination was employed for this purpose. The free surface was then adjusted accordingly. (ii) The velocity field in the domain was then obtained by reapplication of eqn. (19). (iii) The stream function was then evaluated by solution of Poisson’s differential equation using a similar boundary integral formulation. The grid pattern was adjusted so that the domain cell vertices lay on given streamlines. (iv) The extra stresses at the vertices were then obtained, using a Runge-Kutta integration package. This represents one cycle of iteration and was repeated until little change was observed. It should be noted that a major advantage of the boundary element formulation is the relatively small system of equations to be solved. This depends solely on the boundary discretization and not on the number of domain cells. Furthermore, it is not necessary to compute the pressure at

154 interior locations, unless this is specifically required. The streamline integration scheme for computing the extra stresses is very efficient, since it requires a simple integration package. This uses very little storage space and only a fraction of the time required for each iteration is dedicated to evaluation of the extra stresses (- 4% in the case of the die swell problem described in the next section). Furthermore, the method is easily adapted to any constitutive equation of the form (4). However, the need to compute the stream function (or particle paths) adds to the overall computational effort required. Despite this, the ability to implement the boundary element method on a small computer (with negligible computer time costs) makes the approach attractive. 6. The die swell problem The geometry of the planar die swell problem is illustrated in Fig. 3. Also shown are the velocity and stress boundary conditions. In addition to these, suitable boundary conditions for the stream function were prescribed. The geometry was chosen such that h, = 1, x~,~,~ = - 2 and ~,,~~x = 3. The velocity and stress profiles adopted at the upstream boundary were the fully developed profiles for pressure driven flow. In the case of the Maxwell and Oldroyd-B models, the velocity profile is the well known Newtonian profile (ie. quadratic). The corresponding stress field is easily evaluated. This is also approximately true of the MFTT fluid, with n = 1, for Wi c 1.5 (in this range the shear viscosity is approximately constant). We also found that it is possible to solve exactly the equations for the Leonov model (/? = 0, $, = 0) to give:

where

P is the pressure

gradient.

Again,

the stress

field

is then

easily

Xl.HAX

Fig. 3. The geometry fluid.

and boundary

conditions

for extrusion

of a plane sheet of viscoelastic

155 evaluated. In other cases the velocity profile was obtained numerically. The governing equations were cast in a form that allowed the velocity gradient (shear rate) to be evaluated at any value of x2. This was then integrated numerically, with the initial condition ul( h,) = 0, to yield the required velocity (and stress) profile. Two measures of the importance of elastic effects were used in this work. The first, the Weissenberg number Wi, was defined as Wi=Xqw,

(21)

where qb is the shear rate at the wall of the die far upstream:This definition can be applied to all elastic fluids, however it is often not a reliable measure of the true importance of elasticity. For example, Wi is not affected by the addition of more Newtonian solvent to an Oldroyd-B fluid, provided that the total viscosity is kept constant, despite the fact that the first normal stress difference has been reduced by this action. A more informative measure to use in these circumstances is the recoverable shear at the die wall, S,, (see eqn. 14). Both of these measures were used by Crochet and Keunings [9] when studying the die swell of an Oldroyd-B fluid. The relationship between Wi and S,, for the various fluids is given in Fig. 2. Note that S,, = Wi for the Maxwell fluid. Clearly, S, is not a suitable measure to use with the PIT fluid ( nls # 0), and we must rely on Wi in this case. The parameters of all models were chosen to give a zero shear rate viscosity of unity. For convenience, in the case of the Maxwell and Oldroyd-B models the upstream velocity profile was chosen to give a value of u1 = 1 on the centreline, and the Weissenberg number was varied by altering X. In other cases A was set to unity and Wi was adjusted by changing the flow rate. The discretization scheme employed is shown in Fig. 4. This consists of 58 boundary elements and 238 domain cells. As previously described, the grid pattern is adjusted during solution so that all the domain cell vertices lie on given streamlines. Thus, the grid lines extending in the main flow direction

Fig. 4. The mesh pattern cells.

used for the numerical

solution:

58 boundary

elements,

238 domain

156 finally correspond to streamlines. This is exemplified in Fig. 4. Note that two of the grid lines are terminated within the solution domain. This is merely to prevent overcrowding of domain cells at downstream locations, where a fine discretization is not required. A grid similar to this was used to treat the extrusion of a Maxwell fluid in [8]. The results were compared with a comprehensive summary of other reported numerical results. The comparison will not be repeated here. All simulations were run on a Perkin-Elmer 3220 computer, in double precision. The Weissenberg number was increased in regular steps; the solution at the previous Wi was used to start the next series of iterations. The step used depended on the model. For example, Wi was increased in steps of 0.25 for the Oldroyd-B model, while steps of 0.5 were used for the PTT model and steps of 1.0 used for the Leonov model. Generally, smaller steps were used at higher values of Wi, if the numerical scheme showed signs of instability. Each case was run for 20 iterations. In general, satisfactory convergence was obtained within lo-15 iterations. If significantly more than 20 iterations were required the Weissenberg number was not increased further. In the case of the Leonov and PTT models the process was stopped at Weissenberg numbers of 5 and 3 respectively, not because of poor convergence but because it was deemed that no further useful information could be gained by continuing the process. Each iteration took approximately 30 minutes on the Perkin-Elmer computer. 7. Results The computed values of the swelling ratio x = h,/h, are plotted in Figs. 5 and 6. There is clearly a wide range of behaviour exhibited by the models tested. The results corresponding to the Maxwell and Oldroyd-B models are identical when plotted against Wi. When increasing Wi in steps of 0.25 we were unable to achieve convergence with the Maxwell model above Wi = 0.75, as illustrated in Fig. 6. The Oldroyd-B model failed to converge above Wi = 1.2 (with a step of 0.2). The reason for the failure in both cases was the emergence of numerical instability (due to the numerical formulation and, perhaps, inherent instability of the model, as discussed in the Introduction) and, in particular, the occurrence of large negative values of Eli. This can be understood by reference to the Trouton viscosity shown in Fig. 1. It was noticed that the flow field near the lip became more like an elongational flow as Wi was increased. That is, the diagonal components of the rate-of-strain tensor increased in size, while the shear components diminished. By analogy with the pure elongational flow case illustrated in Fig. 1, the strain rate exceeded the critical value at which qr is undefined. Beyond this point 7i1 becomes negative. This problem was apparently not encountered by Crochet

157

1.30

X 1.25

1.10 0

1

3

2

4

5

Wi Fig. 5. Computed Wi. For legend, Keunings [9].

1.101 0.0

0.5

values of the swelling ratio, x, as a function of the Weissenberg number, see caption to Fig. 1. Also: (9) Oldroyd-B (q,,/q, = 8), Crochet and

1.0

1.5

2.0

2.5

3.C

SR

Fig. 6. Computed values of the swelling ratio, x, as a function of the recoverable shear, S,,.,. For legend, see caption to Fig. 1. Also: (9), Oldroyd-B (11,/q, = 8), Crochet and Keunings 191.

158 and Keunings [9] since they were able to reach Wi = 4.5 using the Oldroyd-B model. It is likely that the different solution method (finite element) and the discretization employed are responsible for the different experience. For reference, we have also plotted the results of Crochet and Keunings on the figures. There is clearly a disagreement between the present results and those of Crochet and Keunings. This demonstrates the effect of different discretization and solution schemes on the results (even at low JR, see also [S]) and illustrates the need to use a single mesh and solution method in the model comparison work reported here. The MPTT model (with n = 1) produced results very similar to those of the Oldroyd-B model in the range of Wi shown. Again this can be understood by reference to the steady shear and elongational behaviour of the model. Firstly, the steady shear behaviour of the MPTT (n = 1) and Maxwell/Oldroyd-B models is nearly identical in the range A+ < 2. Secondly, the behaviour of the elongational viscosity of the models is similar, in the sense that qr rises sharply to a high value as the strain rate is increased (for this choice of model parameters). The similarity between the performance of the MPTT model and Maxwell/Oldroyd-B model is, therefore, expected. It is important to note that ,numerical instability began to interfere with the solution for Wi > 1.2. However, with some perseverence solutions up to wi = 10 have been obtained. These are shown in a later figure. We encountered no problems associated with the appearance of non-physical solutions, as we did with the Maxwell/Oldroyd-B models. The response of the PTT fluid is disappointing, but not unexpected. Tanner [23] has demonstrated that die swell can be predicted (semi-empirically) on the basis of the recoverable shear at the die wall. The recoverable shear is therefore important in die swell, as well as the elongational behaviour near the lip. By reference to Fig. 2 it can be seen that the rate of increase of S,,( q7s# 0) with shear rate begins to decrease significantly for Wi > 1.5. The slope of the swelling ratio curve is also seen to steadily decrease in this range. Although solutions for Wi > 3 could be obtained, this was not deemed productive. Notice that values of x for the PTT model at a given value of Wi are consistently less than the corresponding x for the Maxwell, Oldroyd-B and MPTT models. In their computer study of die swell of an inelastic power-law fluid, Tanner et al. [24] demonstrated that shearthinning effects tend to reduce the expansion ratio relative to the Newtonian case. In the case of the PTT fluid, the small reduction of x relative to a similar fluid that does not exhibit shear-thinning is expected. The behaviour of the Leonov fluid (j3 = 0) is less easily understood. The recoverable shear varies as m at large Wi, while for the other models it varies as Wi. However, in the range Wi -c 1 the choice of model does not whereas the expansion ratio for the Leonov model significantly affect S,,

159 behaves very differently to that of other models. Note also that shear-thinning behaviour is not unusual in this range of Wi. Therefore, the behaviour of the expansion ratio cannot be explained on the basis of S,, and shear-thinning alone. The explanation can be found by returning to the elongational behaviour. By considering the stiffness of elongated filaments on the outside of the extrudate, Tanner [25] has demonstrated that die swell can also be predicted (semi-empirically) on the basis of elongational flow. Again, this is not the whole explanation of swelling since shear effects are neglected. The important result is that if the viscosity in elongation is independent of the strain rate (as is the case for the Leonov fluid with j3 = 0) then there is no swelling relative to the base flow. If the fluid is Newtonian, there is no swelling relative to the Newtonian base flow. In the case of the Leonov fluid, the theory states that there is no swelling relative to the equivalent shear-thinning base flow. At high shear rates the viscosity of the Leonov fluid follows the power-law behaviour: my-‘. Thus the base flow at high shear rate is a plug flow; there is no expansion. At lower shear rates the shear-thinning behaviour is less severe, and the observed expansion of approximately 12% is not unexpected. In the range Wi > 2 the expansion ratio is observed to increase, indicating that other effects (such as the influence of S,,) are becoming important. It was not instructive to seek solutions for Wi s- 5, although we believe they could be obtained without difficulty. Further evidence of the importance of the elongational characteristics of a fluid in the die swell problem can be found by observing the case p = 9. The parameter /3 has a remarkable effect on the swelling ratio. Recall that j3 does not greatly affect the shear-thinning behaviour of the Leonov model, but it at large Wi by a factor of J(l+p) and increases the increases S, maximum ql- by a factor of (1 + p). The increase of S,, with /I cannot wholly explain the larger die swell. As seen in Fig. 6, at a given value of S,, the die swell with j3 = 9 is still considerably greater than that for the case /?= 0. Thus the improved elongational behaviour has a marked effect on the die swell, and the Leonov model with p # 0 shows promising results. Again, we have not proceeded beyond Wi = 3, although we believe solutions at higher values of Wi could have easily been obtained. We now have a great deal of freedom to fit the model to elongational and shear data, and to use the model to predict realistic swelling ratios. The effect of the different models on the exit pressure loss, AP is illustrated in Fig. 7. The exit pressure loss is defined as the difference between the computed upstream pressure and the pressure drop that would occur without the presence of the free surface. This has been non-dimensionalised with respect to the wall shear stress TV. It is clear that the models exhibiting the greatest swelling also produce the greatest exit pressure loss at a given value of Wi.

OLlI 0

1

2

3

Wi

Fig. 7. Computed values of the dimensionless pressure loss, AP/q, as a function of the Weissenberg number, Wi. For legend, see caption to Fig. 1. Also: (9) Oldroyd-B (q,/qi = 8), Crochet and Keunings [9].

Fig. 8. Computed values of the expansion ratio for the MPTT model (5 = 0, n = 1, TJ~= 0, q, = 1) for c = 0, 0.005 and 0.01. Note that c = 0 corresponds to the Maxwell model. Also shown are the results of Crochet and Keunings [9] for the Oldroyd-B fluid and the theoretical result by Tanner [23]. - - - - - - Finite element results for PTT model (I = 0, nS = 0.112, q, = 0.888, c = 0.015).

161 TABLE

1

Swelling of MPTT fluid Wi( X?) SR Xii/h, Swelling(%)

0 0 0 18

0.25 0.25 0.083 17

0.5 0.5 0.167 21

0.75 0.75 0.25 26

1.0 1.0 0.333 32

2.0 1.9 0.667 46

3.0 4.0 5.0 7.0 2.7 3.3 3.9 4.8 1.00 1.33 1.67 2.33 57 65 71 77

10.0 5.9 3.33 84

Parameters used here (eqns. 8 and 9) are n = 1, { = 0.0, c = 0.01, ns = 0.0, q, = 1.0. The value of the dimensionless number Xii/h, is also given; here ii is the mean speed in the duct.

Finally, in Fig. 8 and Table 1, we show results obtained using the MPTT model at higher Wi ( Wi < 10). These simulations were run on a Vax 11/780 and an Apricot microcomputer. Up to 60 iterations were required to achieve satisfactory convergence at high values of Wi. Note that the present numerical results confirm the trend predicted by Tanner’s theory [23]. Also shown are some finite element results obtained using the program AXFINR [8,24] and the PTT model with 5 = 0 and 6 = 0.015. These results also confirm the trend. The boundary element results for the case E = 0 (Maxwell model) at Wi > 1 were obtained by employing the solution for c = 0.005 as a starting point. We believe these are the first reliable high Weissenberg number results for the Maxwell model. Non-physical solutions did not arise at these values of Wi, indicating that solutions at higher Wz’ using the Maxwell (or Oldroyd-B) model can be obtained by employing very small increments in Wi. It is more practical to use a model that does not suffer from this limitation. 8. Conclusions The following conclusions have been drawn from this work: (i) The performance of a viscoelastic model on the numerical die swell problem can be related to the steady elongational and shear characteristics of the model. Both types of behaviour are important in the die swell problem. (ii) The Maxwell and Oldroyd-B models are capable of predicting significant die swell ( - loo%, see [9]), however the experience gained in the present work indicates that these models are unsuitable for high Weissenberg number die swell calculations. (iii) The MPTT and modified Leonov models show the most promising behaviour. It is expected that significant die swell can be obtained with these models in the future, as well as realistic modelling of actual rheological fluids. The most important of these conclusions is the last one. Having identified

162

the encouraging performance of these models at low wi, a more intensive study of the die swel1 of MF’TT and Leonov fluids is now required. Acknowledgments

This work was supported by the Australian Research Grants Scheme. We are grateful for this support. We also wish to thank Dr. F. Sugeng who assisted with some of the computations. References 1 ‘2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

K.R. Reddy and RI. Tanner, J. Rheology, 22 (1978) 661. P-W. Chang, T.W. Patten and B.A. Finlayson, Computers and Fluids, 7 (1979) 284. M.J. Crochet and R. Keunings, J. Non-Newt. Fluid Mechanics, 7 (1980) 199. M.J. Crochet and R. Keunings, J. Non-Newt. Fluid Mechanics, 10 (1982) 85. C.J. Coleman, J. Non-Newt. Fluid Mechanics, 8 (1981) 261. N.Y. Tuna and B.A. Finlayson, J. Rheology, 28 (1984) 79. B. Caswell and M. Viriyayuthakorn, J. Non-Newt. Fluid Mechanics, 12 (1983) 13. M.B. Bush, J.F. Milthorpe and RI. Tanner, J. Non-Newt. Fluid Mechanics, 16 (1984) 37. M.J. Crochet and R. Keunings, J. Non-Newt. Fluid Mechanics, 10 (1982) 339. N. Phan-Thien and R.I. Tanner, J. Non-Newt. Fluid Mechanics, 2 (1977) 353. N. Phan-Thien, J. Rheology, 22 (1978) 259. N. Phan-Thien, J. Non-Newt. Fluid Mechanics, 16 (1984) 329. AI. Leonov, Rheol. Acta, 15 (1976) 85. A.I. Leonov, E.H. Lipkina, E.D. Paskhin and A.N. Prokunin, Rheol. Acta, 15 (1976) 411. R.G. Larson, Rheol. Acta, 22 (1983) 435. M.J. Crochet, A.R. Davies and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier Science Publishers, Amsterdam, 1984. N. Phan-Thien, J. Non-Newt. Fluid Mechanics, 13 (1983) 325. N. Phan-Thien, J. Non-Newt. Fluid Mechanics, 17 (1985) 37. R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymeric Liquids: Vol. 1, Fluid Mechanics, J. Wiley and Sons, New York, 1977. M.B. Bush and R.I. Tanner, Int. J. Numerical Methods in Fluids, 3 (1983) 71. M.B. Bush and N. Phan-Thien, J. Non-Newt. Fluid Mechanics, 16 (1984) 303. O.A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Flow. Gordon and Breach, U.S.A., 1963. R.I. Tanner, J. Polymer Science, A8 (1970) 2067. RI. Tanner, R.E. Nickel1 and R.W. Bilger, Computer Meth. Appl. Mech. Engin., 6 (1975) 155. R.I. Tanner, J. Non-Newt. Fluid Mechanics, 6 (1980) 289.