Extrudate swell and isothermal melt spinning analysis of linear low density polyethylene using the Wagner constitutive equation

lkdmks ELSEVIER

J. Non-Newtonian Fluid Mech., 69 (1997) 113 136

Extrudate swell and isothermal melt spinning analysis of linear low density polyethylene using the Wagner constitutive equation R. Fulchiron a,,, p. Revenu b, B.S. K i m b, C. Carrot b, j. Guillet b ~'Laboratoire des Mat~riaux Plastiques et des Biomat~riaux, URA CNRS 507, UniversitO Claude Bernard Lyon 1, 43 Bd du 11 Novembre 1918, 69622, Villeurbanne COdex, France b Laboratoire de Rh~ologie des Mati~res Plastiques, Facult~ des Sciences et Techniques, UniversitO Jean Monnet, 23 rue du Dr. Paul Michelon, 42023, Saint-Etienne COdex 2, France Received 13 December 1995; revised 20 September 1996

Abstract The rheological behaviour of a linear low density polyethylene melt is studied using a wide range of techniques: shear and elongational stress growth at constant strain rate, steady shear flow, extrudate swell and isothermal melt spinning. The analysis of the behaviour is based on the Wagner constitutive equation which appears to be suitable. Moreover, the extrudate swell and the isothermal spinning are simulated by taking into account the past shear and elongation deformation in the reservoir, the contraction and the die of a capillary rheometer. For this purpose, the Finger deformation tensor is evaluated using a Protean coordinate system. This approach is shown to be very useful to test the validity of the constitutive equation and to determine its characteristic parameters. © 1997 Elsevier Science B.V. Keywords: Elongation; Extrudate swell; Isothermal melt spinning; Linear LDPE; Protean coordinates; Shear; Viscoelasticity; Wagner model

1. Introduction

The analysis of the viscoelastic behaviour of polymer melts is of importance in polymer processing operations. Among the constitutive equations used to describe flow phenomena, the major difficulty is the choice of an equation relevant to a wide range of deformations and strain rates, both in shear and in elongation. Furthermore, especially in elongation, it remains difficult to obtain experimental data with homogeneous deformation at large constant strain rates. That is why many researchers have used isothermal melt spinning experiments to overcome this * Corresponding author. 0377-0257/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 7 - 0 2 5 7 ( 9 6 ) 0 1 527-3

114

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

problem. However, in this case, other difficulties occur when comparing experimental data with the model predictions. These difficulties arise from the inhomogeneity of the deformation, from the non-constant strain rate along the spin-line and, moreover, from the extrudate swell phenomenon. Indeed, in a spinning operation, before the extensional deformation within the spin line, the polymer melt undergoes both elongation and shear strain in the converging die entrance and in the die itself. Consequently, because of the elasticity of the material, these previous complex deformations influence the extensional behaviour of the melt after the die exit. Our work is based on the Wagner integral constitutive equation [1] which was widely tested in simple elongation and shear [2-9]. Recently, different research groups have published calculation methods with integral models of the K-BKZ type to simulate exit flows from a contraction [10-17]. Some papers deal with the simulation of the melt spinning process [18-20], where the researchers took into account the extensional strain along the spin line, but they generally considered quite unrealistic deformations before the die exit so that they cannot take into account the extrudate swell. Considering the extrudate swell problem, a few papers deal with the flow simulation of memory integral fluids at the exit of long and short dies [10-13]. Some of the cited papers present comparisons between experimental and numerical results. Goublomme et al. [11] gave numerical results for the extrudate swell problem for a long capillary die. In their calculation, they used a K-BKZ constitutive equation involving a spectrum of relaxation times and a damping function either of the Wagner [4] or of the Papanastasiou [6] type. If the calculation provided consistent results for a long capillary die, the assumption of a short die, which necessitated including the upstream converging section in the calculations, resulted in largely overestimated values of the swell ratio. The numerical data presented in a second paper [12] proved that, if extrudate swell was practically unaffected by non-isothermal flow conditions, it was considerably modified by change of parameters in the K-BKZ equation. A non-zero constant ratio of the second and first normal stress differences led to a significantly lower computed extrudate swell ratio than the experimental data, as a consequence of modification of the extensional behaviour of the fluid. The experimental and numerical results reported by Luo and Tanner [10] are also of particular interest since they are related to a thoroughly characterized IUPAC low density polyethylene (LDPE) melt. In order to ensure good predictions in mixed flow situations, they relaxed the separability assumption in the K-BKZ form adopted in their study, at least in elongation, allowing a parameter of the damping function to vary according to each relaxation time of the memory function. Predictions close to experimental data were obtained at low shear rates (~',~ 1 s 1), even with very short dies. However, when the shear rate was increased, the numerical results were found to be larger than the experimental results. As pointed out later by Goublomme and Crochet [12], neither the non-isothermal flow assumption nor the arbitrary decrease in normal stress ratio can affect the numerical swell predictions sufficiently to eliminate completely the discrepancy between numerical and corresponding experimental swell ratios. A similar study is reported in the recent paper of Barakos and Mitsoulis [13]. It is worth mentioning the simulations for the same IUPAC LDPE melt have been undertaken using the same constitutive equation with parameters determined by Luo and Tanner [10]. Comparisons are made with extrusion experiments given by Meissner [21]. Concerning calculations for short and long dies, the trends are similar to those observed by Luo and Tanner [10]. If the agreement

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113 136

115

is generally good for short dies and low shear rates, meaningful overestimation becomes larger and larger as the die length-to-diameter ratio and the apparent shear rate increase. For the purpose of comparisons between theory and experiments, in these papers the need is pointed out for significant experimental studies related to extrudate swell ratio measurements for polymer melts characterized in shear and in elongation. Progress in flow computation in domains involving the upstream and downstream channels together with the exit region is also necessary. In order to obtain a better insight into these major issues, the purpose of the present work is to report some experimental and numerical results on flows at the exit of an axisymmetric contraction for a linear LDPE (LLDPE) melt. The analysis is carried out either with a melt spinning process (tensile force different from zero) or with an extrudate swell experiment (without tensile force). The whole past deformation before the die exit is considered in order to assess both the velocity profile and the force in the spin line. A memory-integral Wagner-type model is set up to describe the rheological properties in simple flow. From a numerical viewpoint, comparisons with experiments are made with results provided by a numerical method involving the Protean coordinate system developed by Duda and Vrentas [22].

2. Calculation method

2.1. Constitutive equation The Wagner constitutive equation [1] is widely used to describe non-linear viscoelastic experimental data. This model is an improvement of the Lodge rubber-like liquid theory [23] which is based on the assumption that the number of network strands is constant whatever the strain magnitude is. Therefore, this model is expressed only with a time-dependent memory function which takes into account the probability of a network strand surviving between the creation and the observation instants (respectively t' and t). In the case of the Wagner model, a second function is included called the damping function, which is the probability of a network strand surviving the relative deformation between t and t'. Following a separability assumption, these two decay mechanisms for network strands are considered to be independent, the first being only time dependent and the second only strain dependent. According to the Wagner constitutive equation, the stress tensor is expressed by

~r(t) = - p l + I'

m(t -- t')h(I~, I2)Cc '(t') dt'

(1)

d

where p is the isotropic pressure, I the unit tensor, C t l(t') the relative Finger strain tensor between t and t', r e ( t - t ' ) the memory function and h(It, 12) the damping function which depends on the first and the second invariants of the Finger tensor (respectively I1 and 6). The memory function is obtained from linear viscoelastic experiments. Considering a discrete relaxation spectrum with n components, this function is defined by

re(t) =

g~ exp i=1

Ti

-

dt

,

(2)

R. Fulehiron et al. / J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

116

where vi are the relaxation times, gi are the modulus contributions of each relaxation time and G(t) is the relaxation modulus. Many mathematical expressions have been used to describe the damping function for shear or elongational flows [1-3]. Papanastasiou et al. [6] proposed a sigmoidal form of the damping function which was later modified by Soskey and Winter [7]. In this function, the deformation is expressed by the general invariant I which is a linear combination of the first (I0 and the second (12) invariants of the Finger tensor: 1

h(I) = 1 + a ( I - - 3)b/2

(3)

I = flit + (1 - fl)I2,

(4)

with

where I, = 12 = 3 + 72 for simple shear flows and I~ = e2*:+ 2e-': and 12 = e-2,:+ 2e,; in simple elongational flows. This damping function is interesting since it may be used both in shear and in elongation. In the case of shear, where the first and the second invariant are the same, the fl parameter in Eq. (4) becomes irrelevant. So, this parameter is obtained only from elongation experiments. This function was applied to PDMS, PE and PS by Papanastasiou et al. [6] and to PP by Fulchiron et al. [9]. Moreover, as already mentioned, the damping function characterizes the deformation-dependent survival probability of a strand. Thus, this function decreases when the deformation increases. However, the rupture of a network strand is clearly an irreversible phenomenon. Obviously, reducing the deformation after a strand disappearance will not lead to a re-entanglement. This irreversibility assumption of network disentanglement, first introduced by Wagner and Stephenson [24], is expressed by a functional ~ ( I ) instead of the damping function h(I) in Eq. (1): t" _

t

)(f[I(t", t')] = ,,m=m {h[I(t", t')]}.

(5)

This functional takes the smallest value of h (I) between the times t' and t and thus the value of h(I) which corresponds to the highest value of the general invariant/. 2.2. The spinning problem The isothermal spinning device is schematically shown in Fig. 1. The molten polymer is extruded from the cylinder (diameter De) through the die (diameter Dd) using a constant speed capillary rheometer. In the converging part, an elongational deformation appears in addition to shear. Then, in the die, the shearing deformation is developed. At the die exit, extrudate swell occurs because of the elastic strain recovery of the material and the rearrangement of the velocity field. Using a take-up device, the filament is drawn down through an isothermal chamber with a take-up speed greater than the extrudate velocity. A constant tensile force and a stationary diameter profile are set up in the thread line. Afterwards, the filament is quenched in a water bath. In order to simulate a spinning operation, we focus on the diameter of the filament inside the isothermal chamber and on the tensile force F.

R. Fulchiron et a l . / J . Non-Newtonian Fluid Mech. 69 (1997) 113-136

117

2.3. Kinematic assumptions In the present work, the computation concerns the velocity profile in the thread line. Nevertheless, since an integral constitutive equation is used, an unbounded integration occurs in the expression of the stress. In other words, the calculated velocity profile is not only a consequence of the drawing ratio but also a consequence of the past deformation above the die exit. In the previous works of Papanastasiou et al. [18,19] and Fulchiron et al. [20], the extrudate swell was not taken into account. Indeed, the 'historical' effects of the flow in the die, necessary to calculate the stress at the inlet, were drastically approximated. In order to take into account a more realistic past deformation, the Finger strain tensor must be evaluated including the deformation in the converging entry of the die and in the die itself. In the present work, the past velocity field in this zone is prescribed, i.e. it does not result from a calculation. However, these kinematic assumptions seem to be realistic. The prescribed velocity field above the die exit is based on the Poiseuille tube flow of a pseudoplastic power law fluid (Fig. 2). The origin of the z coordinate is located at the exit of the capillary. Here, another material parameter is introduced: the power law index n. In the cylinder the converging part and the capillary, the axial velocity is expressed by

u(r, z) - 3n +______~1 Q n + 1 ~R(z) ~

{

1

(6)

-

LR(z)I

Y O~

Lo

Ld (Diameter Dd)

,,okm~I I / / ~'so'"errna'

j ,came ke-upDevice

~[~

Scales

Fig. 1. Sketch of the isothemal spinning process.

118

R. Fulchiron et al. /J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

Fig. 2. Velocity field description. w h e r e Q is the v o l u m e t r i c flow rate, R ( z ) = Rc is the cylinder radius f o r z < - ( L c + Ld), R ( z ) = -- (z + Ld) t a n ( 7 / 2 ) + Rd is the e v o l v i n g radius in the c o n v e r g e n c e f o r - (Lc + Ld) < z < --Lo (~ is the c o n e angle, Lc a n d Ld are respectively the c o n v e r g e n c e a n d the die lengths) a n d R ( z ) = Ro is the die radius for - Lo < z < 0. In the c o n v e r g e n t p a r t , such a n a s s u m p t i o n f o r axial velocity m a y be a subject o f c o n t r o v e r s y b e c a u s e the h y d r o d y n a m i c l u b r i c a t i o n a p p r o x i m a t i o n [25] is n o t valid. A c t u a l l y , the c o n e a n g l e is t o o large (in o u r case 90°). H o w e v e r , the velocity field in this d o m a i n is o n l y used to Table 1 Relaxation spectrum of linear low density polyethylene at 160° T~(s)

g, (Pa)

1.28 × 10 4 6.12× 10 3

1.849 x 10 + 6 2.199 x 10+5

4.10x

8.202×

10 - 2

2.77 × 10 I 2.01 × 10° 1.57 × 10+j 1.35 × 10 +2

10 + 4

1.694 × 10 + 4 1.854 x 10 +3 1.278 x 10 +2 7.084 × 10 +°

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113 136

119

evaluate the Finger tensor at the exit of the die and not to calculate the pressure drop in convergence. F r o m the expression for the axial velocity and using the continuity equation for incompressible fluids, the radial velocity is obtained: Ou 1 ~,(rv) + ~.z r ~r

0,

v(r,-)- n-~l ~R~) 3 ~zz

(7)

I-L

j

"

(8)

Obviously, the radial velocity equals zero in the cylinder and in the capillary, where d R ( z ) / d z = O. In the convergence, this velocity can be rewritten n + 1 zcR(z) 3tan

1

LR(z)J

(9)

or

v(r, z) = - u(r, z) R--~ tan -~ .

(10)

Along the thread line (z > 0), the axial velocity u(r, z) is considered to be uniform over the cross-section of the fibre (independent of r). Therefore, for z > 0, the continuity equation leads to r du 2 dz"

v(r, z ) -

(ll)

The volumetric flow rate Q being constant along the thread line, the axial velocity is related to the diameter D ( z ) of the filament by

40

u(z) = ~D(z)2 .

(12)

With this velocity field, one can obtain the stream lines which are used to follow the fluid particles. Thus, the location of a given particle only depends on z. These stream lines are defined by r(z) ~P - R (z~) - constant,

(13)

with0<~o
2~

f~

cr(r, z)r dr - Fg(z) + Fi(z) - F = 0,

where F is the tensile force.

(14)

R. Fulehiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

120

10 6 10 5

10 4 10 3 C i0 2 0

i0 1

1© o 1 0 - i10 _ ~

,,1-0.2,,' ,i ........ i 0 -II........ I 0 ° ' ........ i 0 II ........ 1 0 2 ' ....... i0 /

W

aT

(Pad/s)

Fig. 3. Storage and loss moduli of the L L D P E at 160°C: A, experimental data; - -

-, calculated curves.

Fg(z) is the gravitational force expressed by Fg(Z) = pgQ

Ldz'

(15)

: u(z)

and F~(z) is the inertial force: Fi(z) = pQ[u(L) - u(z)],

(16)

where p is the fluid density, g the gravitational acceleration and L the length of the isothermal chamber. o-(r, z) is the first normal stress difference. According to the velocity profile above the die exit, each particle underwent a different historic deformation before the current location z. Consequently, the first normal stress difference depends on the radial position. So, these stresses are s u m m e d over the cross-sectional area to determine the rheological force. Eq. (14) can be rewritten using Eq. (13):

2re ~o a( q~, z )R (z)2~o d~p - Fg(z ) + F~(z) - F = O.

(17)

In Eq. (17), the first normal stress difference ~r(~o, z) is related to the location instead of the time as in Eq. (1). Therefore, all the c o m p o n e n t s of Eq. (1) have to be rewritten as a function of z. Firstly, for each stream line defined by q~, the difference t - t' which appears in the m e m o r y function can be expressed by t-- t ' =

dz"

(18)

, u(q~, z')'

where z' is the location at time t' and z is the location at time t. Thus, the first normal stress difference is written as a(q~, z) = f -

m(cp, z, z')Jt~[I(cp, o~

z , -Tt)][C~__- 1 ((~9, z , _7') -

C ~ r 1((/9, z , _7')] .uT( t(-.dz' ~o ,

"

(19)

R. Fulchiron et a l . / J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

121

In this equation, the main difficulty is to evaluate the Finger tensor in the cylindrical coordinates. This calculation is developed in the next section. 2.4. Calculation of the Finger relative strain tensor

The calculation of the Finger tensor is achieved using the procedure described by Adachi [26,27] which is based on the Protean coordinate system due to Duda and Vrentas [22]. This method is particularly relevant for steady flows in two-dimensional or axisymmetric problems. A material point which occupied the location p' at time t' is assumed to reach the position p at time t. In its general form, the Finger tensor is defined as C - I = g ,k! t~X,,, ~X,,

0x~ c~x'~grog,,,

(20)

where the summations on k, l, m, n (ranging from 1 to 3) are implicit, O 'kl are the metric tensor components at p', g is a base vector at p, and grog, is a dyadic product [25]. The Protean coordinate system is defined as follows: 41

=

"~,

"~2 = ~(r, z),

X3 = 0.

(21)

The corresponding coordinate lines coincide respectively with the stream line, the r coordinate and the 0 coordinate lines. In the case of an incompressible fluid, for axisymmetric flow, the stream function ~//is defined by

ag, Oz

rv,

#r

ru.

(22)

Therefore, the continuity equation is necessarily satisfied. In the following equations, u and v are the axial and radial velocity components at time t (location z, r) and u' and v' are the axial and radial velocity components at time t' (location z', r'). Adachi [26] showed that in this case (axisymmetric flow and Protean coordinate system), only two terms of 0£,,/c~£~, in Eq. (20) are unknown. These terms are expressed by 04~

u

r)4]

u"

041 _ ~ c~u" dz" 0 ~ - u j_, a~" u 'e 022

043

Ohm all other ~-~ = O,

(23)

f

(24)

with =

fl •u" dz" , c~'" u "2 "

122

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

The calculation of the metric tensor g,k~ of Eq. (20) in the particular case of Protean coordinate system (g,k~ becomes g'*~) is achieved by differentiating Eq. (21) and using Eq. (22) which leads to the a2'k/ax'i values: =

It

=

(~X' I

(~X i

all other ~

= O.

r vu t ,

-

-

--

1,

aX;

=

p t

- - r u ~

(25)

ax' 1

In a similar manner, it can be shown that

ax', = 1, as',

ax~ 1 a2'2 r'u" ax2

all other ~ -

ax; = 1, as;

ax~ v' an', u"

--

(26)

0.

This last equation will be used later in this paper. In the case of axisymmetric flows the metric tensor components are [25] t

t

t

g'_-_-= 1, grr = 1, goo = r '2, g'---- = 1, g'r~ = 1, g,OO = 1/r,2,

all other gu = 0, all other g,U = 0.

(27)

Using the tensor transformation rule, Eq. (25) and Eq. (27) lead to the expression of the metric tensor components in the Protean coordinate system: (28)

gtkl = g,~/ OX'k a2}

ax', ax;.' 1

g,kl =

-

r'v'

-- r'v'

r'2(u'2 + v '2)

0

0

07 ~J

(29)

From Eq. (29) and Eq. (23), using Eq. (20), one obtains the Finger tensor components in the Protean coordinate system:

( ~ - I ) 2 2 = r,2(U,2 _{_ V,2), (~-,)33

=

1

__

r~2

,

all o t h e r ( C - l)kt = O.

(30)

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

123

In order to apply Eq. (29) to the spinning problem, the Finger tensor c o m p o n e n t s must be written with respect to the velocity field expressed by Eqs. (6)-(11). Firstly, the integral defined in Eq. (24) is expressed for z' located in each zone of the flow: zone 1, the cylinder z' < - (L~ + Lo); zone 2, the convergence - (Lc + Lo) < z' < - La; zone 3, the capillary - Ld < z' < 0; zone 4, the thread line 0 < z'. To calculate this integral, the term Ou/O~ must be rewritten. Using Eq. (26), one obtains -

-

-

-

0u

0u 1

8~

8r ru"

(31)

Therefore, Eq. (24) becomes =

(32)

, Or" r"u ''3"

Using Eq. (32), the integral term of Eq. (30) can be calculated for each location of z', z being always greater than 0. For example, when z' is located in zone 1, the integral is the sum of four terms: =

:,

Or" r"u "~ +

- - ( L c + L d)

Or" r"u -'--7~+

Ld

8r" r"u ''~3+

8r" r"u ''3

or

f = S1 + $2 + $3 + $4

(33)

Using Eqs. (6)-(11), these terms can be evaluated for a stream line defined by Eq. (1): In zone 1, u is given by Eq. (6) with R ( z ) = Rc so that it is independent of z. It can be easily shown that S1 = K R 2 ( L c + Ld + z'),

(34)

(n + 1)37r2(pl/n- 1 K -- n(3n + 1)2Q2(1 - q) l + ,,,)3-

(35)

with

In zone 2, u is given by Eq. (6) with R ( z ) = expression for $2: $2 = - K

~Cc+ Ld)

R(z) 2 dz = K t a n -

\2/\

In zone 3, u is given by Eq. (6) with R ( z ) = $3 =

-

R d

KR2Ld .

In zone 4, u is independent of r so Ou/& = 0:

-(z+

3

Ld)tan(:¢/2)+ Rd. This leads to the

"

(36)

independent of z, so (37)

124

R. Fulchiron et a l . / J . Non-Newtonian Fluid Mech. 69 (1997) 1 1 3 - 1 3 6

$4 = 0.

(38)

• Finally, if z' is located in zone 1, the integral term is written

f =K [R2(Lc+Ld+z')+tan

\2J\3

'

39,

with K defined by Eq. (35) •

If z' is located in zone 2, one obtains

,40,

f = K{ tan-l(~']FR~\2j[_3 •

If z' is located in zone 3,

f = KR2z'.

(41)

• If z' is located in zone 4,

f

= 0.

(42)

The Finger tensor components in Protean coordinates can now be calculated with respect to Eqs. (39)-(42) for each fluid particle (along a stream line). In order to use the Finger tensor in Eq. (19), one has to transform its components into the cylindrical coordinate system by applying the tensor transformation rule (C-l)~yl = (C- I)kl c~xiOXi c~k 0~t'

(43)

where the terms OXi/O~kare given by Eq. (26) (applied to the location z). Finally, those components are transformed into the physical components using [25] ( C - 1 ) 6 ___(C-l)~yl(giigj/)l/2

(44)

,

where the g~j values are given by Eq. (27). Using Eq. (43) and Eq. (44) one obtains the physical components of the Finger tensor in cylindrical coordinates:

(C-- 1)ZZ= (C-- 1)1I

rv(C-1),1 + (gj- 1),2 (C-1),_- = (C-')-'"

=

(C-1),, = r2v2(~ - ,)lJ +

ru

2rv(C,-1),2 + (0-,)22 r2u 2

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136 (C-1)oo = r=(~-1)33.

125 (45)

Then, the Finger tensor invariants used in Eq. (3) and Eq. (4) are given by It=tr(C-I)

and

12=5, { t r 2 ( C - i )

tr[(C-t)2]}.

(46)

These invariants are used to obtain the generalized invariant I of Eq. (4). It should be noted that, in the computation, the maximum I value between z' and z is retained to define the minimum value of the damping function in this domain (Eq. (5)). 2.5. Velocity profile and tensile force calculations in the spin line 2.5.1. Sp&ning process The known parameters of the problem are as follows. • The spinning conditions are the volumetric flow rate Q, the geometrical values (Re, Ra, Ld, ~, L), the drawing ration between the velocity u(L) at the chamber exit and the velocity Uo at the die exit. The extrusion velocity Uo is given by

Q uo = rcR2.

(47)

• The material parameters are the relaxation spectrum obtained from linear viscoelastic experiments and the damping function obtained from non-linear viscoelastic experiments. The unknowns are the velocity profile from the die exit to the end of the isothermal chamber and the tensile force. The velocity profile in the thread line is computed using a finite difference analysis. The spin line (z > 0) is divided into N elements of length z i - zi_ ~ in which the variation in the velocity u(z) is assumed to be linear. So, the velocity between z~_ 1 and zi is expressed by u(z)-(1

-

~)U~_l+ (1 + ~)u~ , 2 2

with

~ = 2z - (z~ + z~_ l) z~- zs_ l

(48)

The problem is now the calculation of the N - 1 unknown values of u~ at the intersection of the elements since Uo for z = 0 and UN = u(L) for z = L are known. So, for a spinning simulation, there are N unknowns which are the N - 1 unknown velocities and the tensile force F. These values are determined following the procedure previously described by Fulchiron et al. [20]. By applying Eq. (17) to all the locations z~ for i = 0 to i = N, one obtains a system of N non-linear equations. The solution is obtained numerically using the Newton method which converges towards the solution starting from initial values of velocities and force. These initial values can be obtained, for example, from the analytical solution of a newtonian fluid or from the experimental data. The integral terms which cannot be obtained analytically (Eq. (19) with respect to z and Eq. (17) with respect to ~0) are computed using the Gauss-Legendre method. Furthermore, the lower boundary in Eq. (19) is theoretically - o o . Practically, the lower boundary is determined as follows: a primary lower value z] of z' is chosen at the beginning of the convergence to compute the stress at z = 0. Then, a second lower value z~ < z] of z' is

R. Fulchiron et al. ,/J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

126

chosen. The stress at z = 0 is computed again and so on until the lower value of z' does not influence the stress any further.

2.5.2. Swelling experiment The calculation scheme is the same but the last unknown is the velocity /'/N for z = L instead of the force which is known (F = 0). In this case, it should be noted that the drawing ratio is less than I. In addition, the length L is shorter than in a spinning operation (8 mm).

3. Experimental results 3.1. Material The polymer studied is a commercial L L D P E (FN1010) supplied by A T O C H E M , the molecular parameters of which are as follows: weight-average molecular weight Mw = 120 000 g m o l - J ; polydispersity index Ip = 6.3; polymer density at 160°C, 0.76 g cm -3

3.2. Linear viscoelastic behaviour To obtain the relaxation spectrum of the polymer, dynamic shear experiments between parallel plates have been carried out using a Rheometrics mechanical spectrometer (RDA 700). The measurements have been performed from 140 to 200°C in a frequency range from 0.01 to 500 rad s-1 [28]. Time-temperature superposition has been applied to draw master curves of storage (G'(co)) and loss (G"(~o)) moduli at the reference temperature 160°C to enlarge the frequency range. From these experimental data, the discrete relaxation spectrum of the sample at 160°C has been calculated using a procedure described by Carrot et al. [29]. The relaxation times and their modulus contributions are displayed in Table 1. Fig. 3 shows the good agreement between the experimental master curves of G'(o9) and G"(~o) and the curves calculated from the relaxation spectrum using Eq. (49) and Eq. (50). G' and G" values are recovered within an error of less than 3%. 20000 A @ G

Shear

rate

1

s -i

15000 @ @ e L

r-

1 0 0 0 0

--a=OO34b=23 5000

....

0 Time

Fig. 4. Shear stress growth at 1 s - ~ : •

•

a=OO18b=326

I

i

10

20 (S)

• , experimental data; - - ,

, calculated curves.

R. Fulehiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

127

10 4

o ullOab >

---a=O

034

a=0 018

[

b=2 3 b=3 26

. . . . . . . l. O . . -. ., . . . . . . . . . . . . . . . . . . . . . . . . . . . . ShearRate (S

1)

Fig. 5. Steady state shear viscosity vs. shear rate at 160°C: ~ , capillary rheometer: E3, cone and plate rheometer stress growth experiments; *, dynamic viscosity (Cox-Merz rule).

G'(~) = £ gi(~°ri)2 i= 1 1 + ((Dz'i) 2 G"(~o) =

(49)

gi(D.Ci

(50)

i= 1 1 + (~ori)2

3.3. Non-linear viscoelastic behaviour

The parameters of the damping function (Eq. (3)) are obtained from non-linear viscoelastic experiments. Stress growth measurements in shear and in elongation at constant strain rate ) or have been carried out at 160°C [28]. In shear, these experiments were performed on a Rheometrics mechanical spectrometer (RMS 800) using a cone-plate geometry (diameter 25 mm; cone angle, 0.1 rad). In elongation, the experiments have been carried out using an elongational rheometer developed by Muller and Froelich [30] which allows a constant elongation strain rate by applying an exponentially growing velocity to the ends of the sample. From such experiments, the damping function is obtained using the following [2,9]: in shear,

200000 Elongatlon

&

ware

1

65

s

i

150000

• @

i00000

5000C

i

J

1

2

Tlme

/

3

i

4

5

(s)

Fig. 6. Elongation stress growth at 1.65 s-~: • • • , experimental data;

, calculated curves.

128

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

a~2(t)/G(t) - | trl2(t')[m(t')/G2(t')] dr' dO

=

,

with

7=~t;

(51)

in elongation, [~r, l(t) - tr22(t)]/G(t) - | [al,(t') - ~r22(t')][m(t')/G2(t')] dr' h(e)=

,10e2~:,_e-~'

,

with

e=~t.

(52)

However, from a practical point of view, these damping functions may be difficult to obtain. Firstly, as was mentioned before, the fl parameter can be estimated from elongational measurements. Unfortunately, in elongation, the strain domain is generally relatively short. Moreover, the cross-sectional area of the sample decreases when the strain increases, so that the force measurement becomes less and less accurate. This can lead to uncertainties in the stress value. Secondly, the differences in the numerators of Eq. (51) and Eq. (52) can induce spurious results especially when the accuracy of the experimental stress is poor. This lack of reliability of the damping functions led us to test two sets of parameters a, b and ft. The first set was obtained from elongation stress growth experimental data. Indeed, the fits on damping functions obtained from different shear experiments (step strain and stress growth) did not lead to the same parameters. Therefore, the experimental data obtained in elongation were used to determine the a and b parameters. The fl parameter was fixed to 0.02, following the results of Papanastasiou et al. [6], who obtained this value for three different polymers (PDMS, LDPE and PS). Then Eq. (3) was linearized to obtain the two other parameters [28]: a = 0.018,

b = 3.26,

fl = 0.02.

(53)

The second combination of parameters was obtained following a trial and error method. First, the value of a was fixed to 0.034 according to Ref. [28]. Then the values of b and fl were adjusted gradually in order to obtain satisfactory results in the spinning simulations. These results will be shown later. These parameters are a = 0.034,

b = 2.3,

fl = 0.24.

(54)

Furthermore, steady shear flow experiments have been carried out using a capillary rheometer (Instron 3211) with three die diameters (0.75, 1.25, 2 mm) and different length-to-diameter ratios (from 1 to 12.5). The Bagley and Rabinowitsch corrections were applied to the results. These experiments have been performed at three temperatures (140, 160 and 190°C). The flow curve at the reference temperature of 160°C has been built using time-temperature superposition. The validity of the two sets of parameters is tested by calculating the prediction in the different deformation modes and comparing the results with the experimental data. Fig. 4 shows the results for a shear stress growth experiment (i = 1 s - J). The calculated curves are obtained using Eq. (1) applied to a shear stress growth experiment, i.e. o'~2(t) =

fo

m(s)~s

G(t)~t

1 + a(~s) h ds + 1 + a(~/t)h"

(55)

For the two sets of parameters, the agreement between the experimental results and the calculated curves is relatively correct despite the 10% overestimated values for the model. The two damping functions lead to similar results.

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

129

Table 2 Extrudate swell experimental results for linear low density polyethylene ( T = 160°C)

Ld/Dd

Q (m3 s-l)

~"app (S--')

Dmax/'Dd

1 1 1 1 3 3 3 3 5 5 5 5 10 10 10 10 20 20 20 20

7.125 x 10 lo 2.375 x 10 9 7.125 x 10 9 2.375 x 10 -8 7.125 x 10 - " ) 2.375 x 10 9 7.125 × 10 - 9 2.375 x 10 8 7 . 1 2 5 x 1 0 Jo 2.375 × 10 - 9 7.125x10 9 2.375 x 10 s 7 . 1 2 5 x 1 0 Io 2.375 X 10 - 9 7.125 x 10 -9 2.375 x 10 - s 7.125 x 10 - l ° 2.375 x 10 - 9 7.125 × 10 - 9 2.375 x 10 - s

3.716 12.39 37.16 123.9 3.716 12.39 37.16 123.9 3.716 12.39 37.16 123.9 3.716 12.39 37.16 123.9 3.716 12.39 37.16 123.9

1.28 1.368 1.472 1.592 1.224 1.288 1.368 1.472 1.208 1.272 1.336 1.44 1.224 1.28 1.352 1.432 1.224 1.256 1.328 1.388

)'~app is the apparent shear

rate:

);app = 32Q/nD3.

Fig. 5 shows the steady state viscosity at 160°C. The experimental curve exhibits data from three rheometrical measurements, namely steady state shear flow in a capillary rheometer, steady state shear flow in a cone and plate geometry and, using the Cox-Merz rule [31], oscillatory shear flow: r/(~) = It/* (~o)[,,, = #

(56)

The calculated curve is obtained from Eq. (55) applied with t ~ ~ :

rl(;') = fl, ~ 1 +m(s)s a(~s) h ds.

(57)

The same conclusion as previously arises from Fig. 5: both of the damping functions lead to calculated viscosity higher than experimental data in high shear rate range, especially with the parameters of Eq. (54). However, both of the models adequately account for the experimental data. The damping functions are also tested in elongation mode by comparison between experimental data and stress growth calculated curves as shown in Fig. 6. In this case, the prediction of the model is close to the experimental curve after 0.3 s, whatever the damping function is. As a partial result, it can be concluded that it is difficult to distinguish between the two damping functions from the usual shear or elongation experiments since they lead to similar calculated results.

130

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136 1

+ x • o "I

m ~' '~ 1 ] o r~ o 1

L/D= 1 L/D=3 L/D=5 L/D=IO L/D=20

X O

4"

Q × 1

I •

b=3

B=O

26 02

D 1

0

i

I

I

12

14

16

18

experimental

Fig. 7. Experimental and calculated swelling ratios: - parameters are given by Eq. (53).

)

, straight line of equation x = y. The damping function

3.4. Extrudate swell and isothermal melt spinning experiments Capillary experiments were performed at 160°C with a constant speed capillary rheometer (Instron ICR 3211). The axisymmetrically converging test section of half-angle 45 ° is defined as follows. Tile upstream and downstream radii are respectively 4.76 and 0.625 mm so that the entrance contraction ratio is 7.62. In addition, the capillary length-to-diameter ratio La/Da ranges from 1 to 50 and the apparent shear rate can be varied from 3 to 120 s-~, so that flow measurements were generally performed in the shear-thinning region. It should be pointed out that the measurements were limited at a shear rate of 120 s - ~ by the onset of melt fracture. For extrudate swell measurements, the polymer melt was extruded at 160°C downwards directly into air. The diameter of the passing extrudate was then measured at 8 mm from the die exit, using an electro-optical device (ZIMMER OGH). All measurements were made approximately 5 mm back from the lower end so that sagging under gravity could be neglected [32]. Unfortunately, the isothermal chamber had to be removed in order to follow the evolution of the shape of the extrudate at the die exit. Equilibrium extrudate swell data are displayed in Table 2. ~i D

8

@

x+ • 6 _o •

4~ ~D1

0

L/D=1 L/D=3 L/D=5 L/D= 10 L/D=20

x

4

2 xl

a=O b=2 B:O

2

E C~ 1 0 10

+

-2,

r~ °l

~+x ~ +

034 3 24

i

1 2i

Dm~×/D d

Fig. 8. Experimental and calculated swelling ratios. - parameters are given by Eq. (54).

1 4'

16

18

(expePzmental)

• straight line of equation x = y. The damping function

R. Fulchiron et al./'J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

131

In order to obtain complementary information on the evolution of the shape of the extrudate near the die exit, photographs were simultaneously taken for each flow condition. Melt spinning experiments were carried out with the longest die (Ld/Da = 50) in order to minimize the swelling ratio. Then the molten filament passed through a 10 cm long heated chamber kept at 160°C. Two glass windows allow visual observation and measurement of the diameter along the thread line. A water bath was used to quench the material at the exit of the chamber. The filament was finally wound with various take-up speeds in order to ensure different drawing ratios. The tensile force was simultaneously measured with scales supporting the cooling bath (see Fig. 1). Measurement of the diameter along the spin line permits calculation of the axial velocity using Eq. (12).

4. Results

For the computations, the number and the length of the N elements in the thread line were chosen in order to obtain a sufficient accuracy (without too long a calculation time). Practically, near the spinneret, where sharp variation in the velocity takes place, the elements are shorter than in the ending part. As a whole, the computations were carried out with a number N between 14 and 20, using a personal computer (486-33 MHz). The solutions were obtained after three or four iterations.

4. I. Swelling ratios As explained above, the computation of the swelling ratios has been performed by calculating the velocity profile of the extrudate without any tensile force (F = 0), the velocity for z = L being an unknown. Using Eq. (12), the diameter is obtained: { 4Q ~,.'2

D(z) : \~u(z)]

(58) "

20 --t~ "-. E E

,~-

i%

a:O

034

b:2

3

13=0

24

q:l

a:O

034

b:2

3

B:O

24

r,,:O

-a=O

018

b=3

26

[3=C

02

5

n:i

i

o

10 o o M > 0

I

b

I

0

1

2

I

I

3 4 Dlstance

i

5 (mm)

±

6

i

•

8

Fig. 9. E x t r u d a t e swell e x p e r i m e n t : velocity vs. distance f r o m lhe die exit. Q = 2.375 x 10 - ~ m 3 s La/D a = 0.96.

~, D a = 1.3 m m ,

132

R. Fulchiron et a l . / J . Non-Newtonian Fluid Mech. 69 (1997) 113-136 0

6

~. 0 E E ~ 0

5

o'---

. . . .

4

>~ 4J m u o ~ o >

0 2

--a:O 0 2

----a:O ....

©

1

0

0

0

a=O

034 034 018

b:2

3

B:O

24

n=l

b=2

3

B=O

24

n:O

b=3

26

B=O

02

n=l

I

I

I

t

I

I

I

1

2

3

4

5

6

7

Distance

8

(mm)

Fig. 10. Extrudate swell experiment: velocity vs. distance from the die exit. Q = 1.725 × 10 ~o m 3 s Ld/D d = 19.23.

1 Dd =

1.3 ram,

To compare the experimental and the calculated results, the length L of the extrudate is fixed as 8 mm in the computation. The value of the maximum diameter Dm,x is taken at this location. The calculation has been carried out for the two damping functions with a value of the power law index set to 1. The influence of this parameter will be discussed later. In Fig. 7 and Fig. 8, the experimental swelling ratios are compared with the calculated values. In Fig. 7, the values calculated with Eq. (53) appear close to the experimental values but they are relatively scattered. On the contrary, the values calculated with Eq. (54) (Fig. 8) are less scattered but they appear around 10% greater than the measured values. This may be due to an experimental problem. Indeed, the isothermal chamber being removed, the temperature of the extrudate may decrease before the maximum diameter is reached. This should lead to an undermeasurement of the experimental extrudate swell. Nevertheless, in all cases, the predictions of the model are fairly close to experimental data in a very wide range of extrusion conditions (1 < Ld/Dd < 20 and 1 s - l < ~app< 150 s - l ) . In addition, we checked the ability of the model to predict the evolution of the shape of the extrudate at the die exit for strong flow conditions (low length-to-diameter ratio and high flow rate) and for smooth flow conditions (high length-to-diameter ratio and low flow rate). In the case of a short die (Fig. 9), it seems that the correct set of parameters is given by Eq. (53) but the two sets provide close results. On the contrary, in the case of long die and low shear rate (long residence time), for which the difference between the two sets of parameters clearly appears (Fig. 10), the best fit is provided by the values of Eq. (54). For these two last experiments, the diameter of the dies was slightly modified. The new values are indicated in the caption of the figures. Table 3 Isothermal melt spinning: measured and calculated forces with the two damping functions (Eq. (53) and Eq. (54)) Drawing u(L)/u(O) Experimental force (mN)

Calculated force with Eq. (53) Calculated force with Eq. (54) (raN) (mN)

2.3 3.5 4.5

3.8 5.0 5.7

5.3 6.5 7.0

5.2 6.5 7.4

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

• •

Fe×p-7 F~p=6

0 5

mN mN

Fc~ z 5 F ~ = 5

7 0

133

mN mt

,d

•

>"

•

•

4J

H 0 0

• *

~) >

*

~

o

o

• *

**~

o

o

o

b=3 13=0

Distance

26 02

(cm)

Fig. 11. Experimental (O, *, O) and calculated ( ) velocity profiles along the spin line. Q = 7.125 x 10-9 m3 s - - I The damping function parameters are given by Eq. (53). The influence of the power law index value is also shown in Figs. 9 and 10 where the results are displayed for n = 0.5 and n = 1 (for one of the two damping functions). These values of n bound the domain of values which describe the steady shear flow (see Fig. 5). As can be seen in Figs. 9 and 10, the influence o f the power law index is practically not significant (in Fig. 9, the two results are even indistinguishable). This result could appear surprising but it must be recalled that the power law index n only occurs in the velocity field expression above the die exit and its effect on the parabolic shape o f this velocity field is not strong enough to affect drastically the Finger tensor calculation under the die exit. In other words, the deformation in the convergence and in the die (elongation and shear) must be taken into account to describe the extrudate swell correctly but an unrefined expression o f the velocity field is sufficient provided that it is reasonably correct. 4.2. Isothermal melt spinning

For isothermal melt spinning experiments, the velocity profile along the thread line was calculated. As mentioned before, the velocity at the end of the isothermal chamber u(L) is an • ~

Fe×~-7 Fe×~=6

O 5

mN mN

Foal-7 Fc~s=6

4 5

mN m~

u >, M u o @ >

a=O

; Distance

034

b=2

3

•=0

24

lo (cm)

Fig. 12. Experimental (©, *, e) and calculated ( ) velocity profiles along the spin line. Q = 7.125 x 10-9 m3 s - - I The damping function parameters are given by Eq. (54).

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

134

Table 4 Isothermal melt spinning for draw ratios of 2.3 and 4.5: calculated results for two values of the power law index (the damping function is given by Eq. (54)) Drawing ratio Distance from the spinneret (cm) Velocity (cm

2.3

4.5

s

1),

Velocity (cm

s

1) F/ =

Difference (%)

n= 1

0.5

0.3 0.5 0.7 1.1 2 3 4 6 8 9

0.372 0.369 0.371 0.386 0.438 0.506 0.584 0.784 1.067 1.250

0.383 0.380 0.382 0.397 0.448 0.515 0.590 0.788 1.070 1.251

3.0 3.0 3.0 2.8 2.3 1.8 1.0 0.5 0.3 0.1

Force

5.21 mN

5.19 mN

0.4

0.2 0.5 0.7 1.3 2.1 3 4 6 8 9

0.351 0.356 0.365 0.412 0.489 0.592 0.736 1.167 1.878 2.384

0.364 0.369 0.378 0.425 0.501 0.601 0.746 1.175 1.883 2.385

3.7 3.7 3.6 3.2 2.5 1.5 1.4 0.7 0.3 0.04

Force

7.40 mN

7.35 mN

0.7

input in the computation whereas the force is unknown. In Table 3, force measurement data and calculations are compared. For these computations, the value of n was set to 1. The results with the damping function parameters of Eq. (54) are remarkably close to the experimental results whereas the calculated forces with the parameters of Eq. (53) are about 20% lower than the measured forces. Moreover, in Fig. 11 a discrepancy appears between experimental and calculated velocity profiles with the damping function given by Eq. (53). However, the damping function parameters of Eq. (54) lead to calculated velocity profiles which are very close to the experimental profiles (Fig. 12). Actually, these parameters have been adjusted to fit the spinning experimental results but it is remarkable that they provide good results for the three drawing ratios for the velocity profiles as well as for the force values. Moreover, it has been previously shown that these parameters lead to correct results in all other controlled flow situations presented in this work, namely steady shear flow and transient elongational and shear flow. Thus, the parameters which appear the most suitable are given by Eq. (54) which leads to the following damping function:

R. Fulchiron et al. / J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

h (I,, I2) =

1

135 (59)

1 + 0.034(0.241, + 0.7612 -- 3) L ' 5

In order to confirm the non-influence of the power law index, some computations were carried out using a different value of n (n = 0.5). The comparison of the results with n -- 1 and n = 0.5 is displayed in Table 4. Here again, the results are quite similar, showing that the power law index is not very influential (the forces only differ by 0.7%). Furthermore, a complementary computation was carried out to reveal the effect of the past deformation on the spinning results. For this computation, the die length to diameter ratio was set to 5 (instead of 50), all other conditions being the same (particularly, the draw ratio). Fig. 13 shows the difference between the calculated velocity profiles. The calculated forces, indicated in Fig. 13, differ by 10%. F r o m this result, it is obvious that the behaviour in the thread line is influenced by the past deformation and this must be taken into account to compare the experimental results and the prediction of a model (for example, when using isothermal melt spinning as an elongational rheometer).

5. Conclusion

The rheological behaviour of L L D P E has been studied experimentally using diverse techniques. Wagner's constitutive equation appears suitable to describe the behaviour in these different cases of flow (shear and elongation). Nevertheless, the difficulty when using the constitutive equation is to define the right d a m p i n g function. Indeed, it has been shown that different parameter combinations could lead to correct results when calculating the response in the case of the usual experiments such as constant strain rate stress growth or steady shear viscosity measurements. However, for extrudate swell and isothermal spinning experiments, a calculation m e t h o d which takes into account the whole past deformation before the die exit has been developed. The main problem of this calculation, the Finger tensor evaluation, has been solved using a Protean coordinate approach. It should be pointed out that a power law axial velocity profile across the 3

~ 011 0

L/D=50

Fcai:7

L/D=5

F~ai=8

4

mN / 1 m N /

-~

, 5 Dlstance

b = 2

3

S:O

24

iC (cm)

Fig. 13. Calculated velocity profiles for two die lengths (LID = 50 and LID = 5). Q = 7.125 × 10 9 m 3 s - 1; drawing ratio, 4.5.

136

R. Fulchiron et al./J. Non-Newtonian Fluid Mech. 69 (1997) 113-136

section is assumed before the die exit, and a constant axial velocity across the section after the exit. Nevertheless, it is no particular problem to deal with this discontinuity for calculating the Finger tensor. Moreover, the results show that such an assumption is valid. From the simulations of extrudate swell and melt spinning, it has been possible to distinguish clearly the effects of the two damping functions. Thus, it can be concluded that this kind of study could be a helpful tool to determine the characteristic parameters of the non-linear viscoelastic behaviour of polymer melts.

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] [26] [27] [28] [29] [30] [31] [32]

M.H. Wagner, Rheol. Acta, 15 (1976) 136-142. M.H. Wagner, J. Non-Newtonian Fluid Mech., 4 (1978) 39-55. H.M. Laun, Rheol. Acta, 17 (1978) 1-15. M.H. Wagner, Rheol. Acta, 18 (1979) 681-692. T. Raible, S.E. Stephenson, J. Meissner and M.H. Wagner, J. Non-Newtonian Fluid Mech., 11 (1982) 239-256. A.C. Papanastasiou, L.E. Scriven and C.W. Macosko, J. Rheol., 27 (1983) 387-410. P.R. Soskey and H.H. Winter, J. Rheol., 28 (1984) 625-645. M.M.K. Khan and R.I. Tanner, Rheol. Acta, 29 (1990) 281-297. R. Fulchiron, V. Verney and G. Marin, J. Non-Newtonian Fluid Mech., 48 (1993) 49-61. X.L. Luo and R.I. Tanner, Int. J. Numer. Methods Eng., 25 (1988) 9. A. Goublomme, B. Draily and M.J. Crochet, J. Non-Newtonian Fluid Mech., 44 (1992) 171. A. Goublomme and M.J. Crochet, J. Non-Newtonian Fluid Mech., 47 (1993) 281-287. G. Barakos and E. Mitsoulis, J. Rheol., 39 (1995) 193. D.G. Kiriaskidis, H.J. Park, E. Mitsoulis, B. Vergnes and J.F. Agassant, J. Non-Newtonian Fluid Mech., 45 (1992) 63-80. E. Mitsoulis and H.J. Park, Theoretical and Applied Rheology, Proc. 1lth Int. Congr. on Rheology, Brussels, August 17-21, 1992, pp. 385-387. X.-L. Luo and R.I. Tanner, J. Non-Newtonian Fluid Mech., 22 (1986) 61-89. D.G. Kiriakidis and E. Mitsoulis, Adv. Polym. Technol., 12 (1993) 107-117. A.C. Papanastasiou, C.W. Macosko, L.E. Scriven and Z. Chen, AIChE J., 33 (1987) 834-842. Z. Chen and A.C. Papanastasiou, Int. Polym. Process., 2 (1987) 33-40. R. Fulchiron, V. Verney, A. Michel and J.C. Roustant, Polym. Eng. Sci., 35 (1995) 518-526. J. Meissner, Pure Appl. Chem., 42 (1975) 551. J.L. Duda and J.S. Vrentas, Chem. Eng. Sci., 22 (1967) 855. A.S. Lodge, Rheol. Acta, 7 (1968) 379-392. M.H. Wagner and S.E. Stephenson, J. Rheol., 23 (1979) 489-504. R.B. Bird, R.C. Armstrong and O. Hassager, in Dynamics of Polymeric Liquids, Vol. 1, Fluid Mechanics, Wiley, New York, 2nd edn., 1987. K. Adachi, Rheol. Acta, 22 (1983) 326-335. K. Adachi, Rheol. Acta, 25 (1986) 555-563. P. Revenu, Ph.D. Thesis, University of Saint-Etienne, 1992. C. Carrot, J. Guillet, J.F. May and J.P. Puaux, Makromol. Chem., Theory Simul., 1 (1992) 215-231. R. Muller and D. Froelich, Polymer, 26 (1985) 1477-1482. W.P. Cox and E.H. Merz, J. Polym. Sci., 28 (1958) 619-622. J. Guillet and M. Seriai, Rheol. Acta, 30 (1991) 540-548.