Numerical Simulation of Viscoelastic Extrudate Swell Through Elliptical Ring Die

FLUID FLOW AND TRANSPORT PHENOMENA Chinese Journal of Chemical Engineering, 19(1) 10ü17 (2011)

Numerical Simulation of Viscoelastic Extrudate Swell Through Elliptical Ring Die* XU Xingming (䇮᱕᱄)1,2, ZHAO Guoqun (䎫ള㗚)2,**, QIN Shengxue (〜ॽᆜ)3 and WANG Wei (⧁့)1 1 2 3

College of Mechanical and Electrical Engineering, Shandong University of Science and Technology, Tai’an 271019, China Engineering Research Center for Mold & Die Technologies, Shandong University, Jinan 250061, China College of Mechanical and Electrical Engineering, Shandong University of Science and Technology, Qingdao 266510, China

Abstract The numerical simulation of extrudate swell is significant in extrusion processing. Precise prediction of extrudate swell is propitious to the control of melt flow and the quality of final products. A mathematical model of three-dimensional (3D) viscoelastic flow through elliptical ring die for polymer extrusion was investigated. The penalty function formulation of viscoelastic incompressible fluid was introduced to the finite element model to analyze 3D extrusion problem. The discrete elastic viscous split stress (DEVSS) and streamline-upwind PetrovGalerkin (SUPG) technology were used to obtain stable simulation results. Free surface was updated by updating the streamlines which needs less memory space. According to numerical simulation results, the effect of zero-shear viscosity and elongation parameter on extrudate swell was slight, but with the increase of volumetric flow rate and relax time the extrudate swell ratio increased markedly. Finally, the numerical simulation of extrudate swell flow for low-density polyethylene (LDPE) melts was investigated and the results agreed well with others’ work. These conclusions provided quantitative basis for the forecasting extrudate swell ratio and the controlling of extrusion productivity shape. Keywords viscoelastic fluid, extrudate swell, finite element method, polymer extrusion

1

INTRODUCTION

Extrudate swell is a common elastic manifestation existed in polymer processing, which is also known as die swell. The prediction of extrudate swell ratio is of great significance in both theoretical and industrial context [1, 2]. The selection of appropriate numerical scheme has great influence on accurate simulation of flow properties and flow patterns [3]. Much effort has been made in attempting to find robust and stable numerical methods including finite difference method [4, 5], finite volume method [6], finite element method [7], spectral method [8], and boundary integral method [9]. Among them, the finite element method incarnates the advantage of adaptability to the complex geometry boundaries [10]. Based on the numerical simulation methods mentioned above, many researchers have investigated the extrudate swell behavior with different rheological models. Mitsoulis and Heng simulated die swell ratio of Newtonian liquids for diverging, straight and converging dies respectively [11]. They found out that for all resins examined, the diameter swell ratio of converging dies is always higher than that of the diverging dies. Garcia-Rejon et al. did further research on extrudate swell ratio for diverging and converging dies [12]. They found that thickness extrudate swell ratio increases with the increase of die contraction ratio but decreases with the increase of inclination angle and

the length of tapered section. They also compared the extrudate swell ratio of Newtonian liquids with that of non-Newtonian liquids. Mitsoulis worked on extrudate swell ratio of pseudoplastic and viscoplastic fluids [13], and came to the conclusion that the thickness swell ratio increases monotonically as shear-thickening increases, extreme shear-thinning produces no swelling and no exit correction for power law liquids, while for the viscoplastic Bingham model, the thickness swell ratio decreases to zero as viscoplasticity increases. Ganvir et al. proposed a method for the numerical simulation of extrudate swell ratio using an Arbitrary Lagrangian Eulerian (ALE) technique based finite element formulation [14]. They predicted die swell ratio of planar and axisymmetric extrusion with abrupt contraction ahead of the die exit. In actual manufacture, precise geometric dimension is essential to produce qualified final products and die swell is considered one of the primary obstacles in product shape control [15]. As is known to all, when the melt flows out of the die orifice, the melt extends to all directions when die swell happens. However, most of the relevant numerical simulation work done before mainly aimed at two-dimensional or axisymmetric flow problems and most of the cross sectional shape is simply circular, annular or rectangular. The full 3D simulation is necessary to obtain more realistic flow properties and patterns. In the present work, we established the mathematical model of

Received 2009-12-13, accepted 2010-12-12. * Supported by the National Science Foundation for Distinguished Young Scholars of China (50425517) and the Shandong Province Natural Science Foundation (Y2007F59). ** To whom correspondence should be addressed. E-mail: [email protected]

11

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

3D polymer die swell of elliptical ring die with differential constitutive equation. Free surface was updated with decoupled method by updating the streamlines which needs less memory space than Gandharv Bhatara’s method adopted by Guo et al [16]. The numerical simulation was executed with different material and technological parameters and the distribution of velocity, shear stress and first normal stress difference were analyzed.

a given volumetric flow rate Qv . (2) There is no slip along the wall surface ( u 0 ). (3) On the free surface the following conditions are satisfied: (V n)t 0 , (V n)n 0 , and there is also no flow through the surface ( un 0 ). (4) On the downstream boundary surface far from the die, tangential and normal stresses is supposed to vanish as well as tangential velocity [ (V n)n 0 , (V n)t 0 , ut 0 ].

2

4

MATHEMATIC MODELING

The flow of viscoelastic melt extrusion is typically governed by principles of conservation which are expressed in terms of partial differential equations [17]. Considering incompressible, isothermal and creeping steady flow, the governing equations of conservation of mass and momentum are simplified as follows:

’u 0

(1)

’V 0 (2) where ı is Cauchy stress tensor which is defined as

V

 pI  2Ks D  W

(3)

where D is strain rate tensor which is expressed as 1 D ’u  ’T u (4) 2 The Phan-Thien and Tanner (PTT) constitutive equation is adopted to describe the rheological characteristics of the polymer: ’

O W  g (W )W

2K p D

HO tr(W ) Kp

Free surface update

As the position of the free surface is unknown a priori, the surface update scheme turns out to be a challenge as it introduces non-linearity and the number of iteration increases to obtain a solution [21]. The free surface is found in a decoupled way. The contour of the free surface is first assumed, and then it is updated with the streamline equation [22]: dx dy dz (9) ux u y uz Integration in the z-direction only, the new location of the free face is calculated from the velocity field based on the latest velocity field as follows [23]: u dx ³ x dz (10a) uz dy

(5)

where K p is the polymer contributions to the total viscosity of the liquid defined as Kp K0  Ks , and g(IJ) is a stress coefficient function [18]: g (W ) 1 

4.1

NUMERICAL PROCEDURE

(6)

The notation ’ over the stress tensor IJ is the upperconvected time derivative defined as wW  u’W  W L  LTW (7) wt where L ’u  [ p D is called the effective velocity gradient [19].

4.2

3

’  Op ’u  2Ks ’D  ’W

un

f (Qv )

0

(8)

where f (Qv ) is the velocity profile corresponding to

(11)

The discrete elastic viscous split stress (DEVSS) method is adopted by introducing an additional elliptic term and then the momentum equation becomes [24]:

Kr ’  ’u  'uT  Op ’  ’u 2 Kr  Ks ’ ˜ D  ’W

To solve the above governing equations, it is necessary to impose appropriate boundary conditions. The boundary conditions are set as follows [20]: (1) On the inlet boundary surface, the velocity is expressed as follows:

0,

(10b)

Considering both computational efficiency and memory requirement, the penalty finite element algorithm is introduced to solve the nonlinear governing equations [24]. By using the penalty model, the momentum Eq. (2) becomes:

BOUNDARY CONDITIONS

ut

³ u z dz

Finite element formulation

’

W

uy

(12)

A Galerkin finite element discrete form is proposed by taking the same form of weighted function as that of interpretation function. T ³: Kr ’  ’u  'u  Op’  ’u Ni d: ³:  2 Kr  Ks ’D  ’W Ni d: e

e

(13)

where Ni is the isoparametric weighted function which

12

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

is expressed as follows: 1 Ni 1  [i[ 1  KiK 1  ] i] (i 1, 2,3, ˜ ˜ ˜,8) 8 When the convective term becomes dominant, the classical Galerkin weighted residual method will lose its optimal approximation. In order to reduce the adverse effect of convective term u ˜ ’W in the constitutive equation on computational stability, the upwindstreamline Petrov-Galerkin streamline-upwind PetrovGalerkin (SUPG) method is adopted and an asymmetric weighted function Wi is introduced [25, 26]: Wi

Ni 

ku ’N i (uu)

(14)

where the coefficient k is defined by the velocity components at the element center  u[ , uK , u] and the element characteristic length  h[ , hK , h] :

1ª 2 2 2 º1/ 2   (15) u h u h u h    [ [ K K ] ] ¼ 2¬ With the asymmetric weighted function adopted, the finite element formula of simplified Phan-Thien and Tanner (SPTT) model is obtained: 1 § T ³: e O ¨© u’W  W’u  ’uW  2 [ p W  ’u  ’ u  k

 ’u  ’T u W ·¸Wi d:  ¹

³:

e

§ HO · ¨1  K W ¸ tr(W )Wi d: p ¹ ©

T ³: Kp  ’u  ’ u Wi d: e

(16) Table 1

5

RESULTS AND DISCUSSION

5.1

Geometric model and simulation parameters

The geometry of the elliptic ring die is defined by the inner and outer ellipses (Fig. 1). The geometric parameters of computational domain is as follows: the length of longer axis of the outer ellipse a1, the length of shorter axis of outer ellipse b1, the length of longer axis of the inner ellipse a2, the length of shorter axis of inner ellipse b2, the entry length L1 and the extrudate length L2. The origin of coordinates is set at the ellipse center of die outlet. The z axis is along the flow direction, x axis is along shorter axis of elliptical ring cross section and y axis is along longer axis of elliptical ring cross section. The values of geometric parameters and the values of material parameters for low-density polyethylene (LDPE) with the average molecular weight of 25000 are shown in Table 1 [16]. All the numerical simulation results below are calculated with the parameters in Table 1 unless otherwise specified. The finite element meshes adopted in present paper are shown in Fig. 1. As the deformation of the melt around the die exit is relatively severe, the meshes near the die exit are refined. The coordinate direction is shown in the figure, and the origin of coordinate is set on the centre of die orifice. By taking both computational efficiency and calculation accuracy into consideration, only one quarter of the computational domain is used in the simulation due to symmetry. The swell ratio calculated by the finite element simulation is defined as S  S0 B (17) S0

Geometric and material parameters

a1/m

b1/m

a2/m

b2/m

L1/m

L2/m

0.03

0.015

0.02

0.008

0.1

0.16

189 °C Șp/Pa·s

Șs/Pa·s

Șr/Pa·s

Ȝ/s

İ

ȟ

1.0×104

0

1.0×105

0.1

0.1

0

Figure 1 Finite element meshes used in the computation

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

5.2

13

Mesh dependency

In order to study the effects of mesh refinements on the numerical results, four mesh division schemes are carried out as shown in Table 2. Fig. 2 shows the comparison of extrudate swell ratio respectively predicted with three meshes for the flow of PTT viscoelastic fluid. The material parameters are taken from Table 1. A similar distribution tendency of extrudate swell ratio is found from Fig. 2, whereas mesh 1 leads to a lower extrudate swell ratio and mesh 2 leads to a higher extrudate swell ratio. As the meshes reach to a certain number (mesh 3 and mesh 4), the extrudate swell ratio tends to be very much the same. Considering the sufficient accuracy and computational cost-effectiveness, mesh 3 is adopted as our base mesh and all the results to be given below are computed with the mesh 3 unless otherwise stated. Table 2

Mesh characteristics Elements

Nodes

Mesh 1

768

1125

Mesh 2

1080

1500

Mesh 3

1260

1740

Mesh 4

1350

1860

Figure 2 Fitting curves of swell ratio versus z-coordinate with different mesh numbers (Volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s) mesh number:ƶ768;ƻ1080; Ƹ 1260; ͪ 1350

5.3 Distribution of physical quantities 5.3.1 Velocity field Figure 3 shows the deformed finite element meshes after the extrusion. As the geometric model is not axisymmetric, the distribution of the velocity on axial cross section is not equispaced any more as shown in Fig. 4. Fig. 4 shows the distribution of resultant velocity vector of ux and uy on the cross section at z 0.003 m. It can be seen that the velocity near to the wall layer is along the normal direction, and the velocity of a certain layer inside the die is almost zero. It is obvious that the velocity along the external boundary is higher than that along the inner boundary and the velocity at both ends of the longer axis is

Figure 3 Deformed finite element meshes

Figure 4 Resultant velocity vector of ux and uy on z 0.003 m cross section (Volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s)

higher than that at the ends of shorter axis. When the polymer flows out of the die, the velocity of polymer redistributes immediately to show the swelling and restore the elastic deformation generated inside the die. As the curvature radius at the ends of longer axis is smaller, the extensional flow is relatively intensive. So more elastic deformation energy is produced which leads to more unstable flow. 5.3.2 Distribution of shear stress Figure 5 shows the contour of shear stress W xy on z 0.003 m cross section which is similar with the distribution of ux on the same cross section. As is discussed in Section 5.3.1, when the polymer flows out of the die, the polymer tends to develop fully without the constraint of die and the elastic deformation generated inside the die begins to restore. Thus, the velocity varies greatly at the ends of longer axis and velocity gradient is relatively large which leads to serious shearing deformation inevitably. It can be seen from Fig. 5 that the shear stress changes greatly around both the ends of longer axis and its gradient is large as a result, which means that more elastic deformation develops at the ends of longer axis. 5.3.3 Distribution of first normal stress difference Normal stress difference is a characteristic parameter for polymer due to its elastic effect. Usually the value of first normal stress difference is much larger than that of second normal stress difference. Thus, the first normal stress difference is the main normal stress difference making the die swell occur. Fig. 6 shows the contours of first normal stress difference on different cross sections. Due to symmetry, only the right-hand side is drawn. Fig. 6 (a) shows the contours of first normal stress difference near to the die orifice on the yz plane, which suggest that the flow has already reached steady state inside the die, but the first normal stress difference changes greatly around the

14

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

Figure 5 Contour of IJxy (unit: Pa) on z 0.003 m cross section (Volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s)

(a) x 0 m axial section (b) y 0 m axial section Figure 6 Contour of first normal stress difference (unit: Pa) on different cross sections (Volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s)

die orifice. When the deformation of the melt is restored outside the die, the first normal stress difference vanishes quickly. The contour of first normal stress difference near the die orifice on the xz plane is similar as shown in Fig. 6 (b), but the stress gradient around the die orifice is a little smaller. 5.4

Effect of zero-shear viscosity

Viscosity is an important rheological parameter for polymers. As the solvent contributions to the total viscosity of the liquid are negligible, we set it to zero and the value of Șp is equal to the zero-shear viscosity. Different zero-shear viscosities are selected in numerical simulation and we compare the extrudate swell ratio at the outlet surface, as shown in Fig. 7. The swell ratio decreases from 0.433 to 0.425 as the zero-shear viscosity increases from 10000 Pa·s to 50000 Pa·s, but the change of extrudate swell ratio is slight. 5.5

Effect of volumetric flow rate

The fluid is elongated by extensional flow along flow direction and elastic deformation occurs. Once

Figure 7 Effect of zero-shear viscosity on swell ratio (Volumetric flow rate: 1.0×105 m3·s1)

the melt goes out of the die, elastic deformation begins to restore and die swell happens. When volumetric flow rate increases, severer elastic deformation occurs. Fig. 8 shows the contours of shear stress IJxy on z 0.003 m cross section for three different volumetric flow rates respectively. It is obvious that the contour of shear stress IJxy is similar but the maximum value of IJxy increases as the volumetric flow rate increases. Similar behavior is observed in our simulation work for die swell ratio in Fig. 9: when volumetric flow rate increases, die swell ratio increases as well.

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

15

(a) Volumetric flow rate of 1.0×105 m3·s1

(b) Volumetric flow rate of 3.0×105 m3·s1

(c) Volumetric flow rate of 5.0×105 m3·s1 Figure 8 Contour of IJxy (Pa) on z 0.003 m cross section (Zero-shear viscosity: 1.0×104 Pa·s)

the extrudate swell ratio at the outlet surface changes from 0.43 to 0.95, as shown in Fig. 10. It is known that relax time has much to do with polymer elasticity and the increase of relax time results in the increases of polymer elasticity.

Figure 9 Effect of volumetric flow rate on swell ratio (Zero-shear viscosity: 1.0×104 Pa·s) volumetric flow rate/m 3 · s 1 : ƶ 1.0×10  5 ; ƻ 3.0×10  5 ; 5 Ƹ5.0×10

5.6

Effect of relax time

From the stand point of rheology, relax time is a characteristic parameter that reflects elasticity of polymers. When relax time increases from 0.1 s to 1.0 s,

Figure 10 Fitting curves of swell ratio versus z-coordinate with different relax times (volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s) relax time/s:ƶ0.1;ƻ0.5;Ƹ1.0

16

5.7

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

Effect of elongation parameter İ

The material parameter İ in PTT constitutive equation adopted in present paper controls elongation viscosity. Usually, elongation viscosity is 102103 times of shear viscosity. Shear viscosity is of shear thinning while the variation of elongation with shear rate is complex. As the parameter İ increases, the extrudate swell ratio decreases, but the change is slight, as shown in Fig. 11. The variation of extrudate swell ratio with İ is similar with the variation of extrudate swell ratio with zero-shear viscosity.

5.9

Comparison with other works

It is necessary to compare the results of our work with that of others. With the parameters taken by Huang et al. [27], extrudate swell simulation is executed with our computer code. As shown in Table 3, when the shearing rate is low, the simulation results obtained by Huang et al. were similar with their experimental results, but when the shearing rate is high, the deviation became obvious. It is considered that the wall slip was neglected at high shearing rate, which led to higher simulation results. The tendency of simulation results obtained in our work is the same as Huang’s simulation results. As the shearing rate increases, our simulation results are slightly lower than Huang’s. Table 3

Figure 11 Effect of İ on swell ratio (volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s)

5.8

Effect of the shape of die

Considering the effect of the die shape, annulus and elliptical ring dies with the same cross area are adopted to calculate extrudate swell ratio. It can be seen from Fig. 12 that the extrudate swell ratio of elliptical ring die is higher than that of annulus die taking the same parameters except for the geometry. For elliptical die, on the ends of longer axis the curvature is relatively big which goes against forming fully developed flow. Thus, more elastic deformation is produced inside the die and the extrudate swell ratio is higher as a result.

Figure 12 Effect of the shape of the die on swell ratio (volumetric flow rate: 1.0×105 m3·s1; zero-shear viscosity: 1.0×104 Pa·s) die shape:ƶannulus;ƻelliptical ring

6

Comparison of swell ratio with other work [27]

Shearing rate/s1

Experimental data [27]

Simulation results [27]

Present simulation

0.1

1.375

1.300

1.3084

1

1.640

1.595

1.5103

10

1.750

1.965

1.8597

CONCLUSIONS

In this study, a mathematical model of 3D viscoelastic flow through elliptical ring die for polymer extrusion swell was investigated. Location of free surface was updated with streamline equation in a decoupled method which is simpler and needs less memory space than Gandharv Bhatara’s method adopted by Guo et al. Numerical simulation was executed with different technological and material parameters. Extrudate swell ratio reduces with the increase of zero-shear viscosity and the material parameter İ in the PTT constitutive equation. But the influence of zero-shear viscosity and İ on extrudate swell ratio is slight. Higher volumetric flow rate leads to higher extrudate swell ratio. As the relax time increases, elasticity of polymers increases and higher extrudate swell ratio is obtained. At the same average flow rate, it is found that more elasticity deformation is generated in elliptical die than annulus die and the extrudate swell ratio is higher as a result. Finally, we compared our numerical results with that of others using the same parameters, and the results agreed well. The results show that the method for free surface update adopted was suitable to obtain accurate numerical result. And further more, calculation time was reduced as less memory space was needed. In this paper, we only compared our numerical results with other’s experimental results due to the limitation of present experimental situation. In the future, more experiments will be carried out in succession to verify the feasibility and effectiveness of our numerical simulation method.

Chin. J. Chem. Eng., Vol. 19, No. 1, February 2011

NOMENCLATURE B hȗ hȘ hȟ I n p Qv S S0 t u ux uy uz uȘ uȗ uȟ İ Șr Șs Ș0 Ȝ Ȝp ȟ, Ș, ȗ ȟp IJ ȍe

’

swell ratio element characteristic lengths in the ȗ direction, m element characteristic lengths in the Ș direction, m element characteristic lengths in the ȟ direction, m unit tensor unit normal vector to the surface of fluid pressure, Pa volumetric flow rate, m3·s1 cross-sectional area of the extrudate, m2 cross-sectional area of the die, m2 unit tangential vector to the surface of fluid velocity vector, m·s1 velocity components in the x direction, m·s1 velocity components in the y direction, m·s1 velocity components in the z direction, m·s1 velocity components at the element center in the Ș direction, m·s1 velocity components at the element center in the ȗ direction, m·s1 velocity components at the element center in the ȟ direction, m·s1 parameter limiting the elongational viscosity of the fluid reference viscosity, Pa·s solvent contributions to the total viscosity of the liquid, Pa·s zero-shear viscosity, Pa·s relaxation time, s penalty constant local coordinates slip parameter which determines the shear behavior of the model extra stress tensor, Pa element region gradient operator

REFERENCES 1

2 3 4

5 6 7

Sombatsompop, N., O-Charoen, N., “Extrudate swell behavior of PS and LLDPE melts in a dual die with mixed circular/slit flow channels in an extrusion rheometer”, Polym. Adv. Technol., 14, 699710 (2003). Sombatsompop, N., Intawong, N., “Extrudate swell and flow analysis of polystyrene melt flowing in an electro-magnetized die in a single screw extruder”, Polym. Adv. Technol., 16, 505514 (2005). Carneiro de Araujo, J.H., Ruas, V., “A stable finite element method for the axisymmetric three-field Stokes system”, Comput. Methods Appl. Mech. Eng., 164, 267286 (1998). Tome, M.F., Grossi, L., Castelo, A., Cuminato, J.A., McKee, S., Walters, K., “Die-swell, splashing drop and a numerical technique for solving the Oldroyd-B model for axisymmetric free surface flows”, J. Non-Newt. Fluid Mech., 141,148166 (2007). Tome, M.F., Castelo, A., Ferreira, V.G., McKee, S., “A finite difference technique for solving the Oldroyd-B model for 3D-unsteady free surface flows”, J. Non-Newt. Fluid Mech., 154, 179206 (2008). Huang, X., Phan-Thien, N., Tanner, R.I., “Viscoelastic flow between eccentric rotating cylinders: Unstructured control volume method”, J. Non-Newt. Fluid Mech., 64, 7192 (1996). Eleni, T., Georgios, C.G., Evan, M., “Numerical simulation of the extrusion of strongly compressible Newtonian liquids”, Rheol. Acta,

8 9 10 11 12 13 14

15

16

17 18 19 20 21

22 23 24

25

26 27

17

47, 4962 (2008). Beris, A.N., Armstrong, R.C., Brown, R.A., “Spectral finite-element calculations of the flow of a Maxwell fluid between excentric rotating cylinders”, J. Non-Newt. Fluid Mech., 22, 129167 (1987). Fan, X.J., Phan-Thien, N., Zheng, R., “A direct simulation of fibre suspensions”, J. Non-Newt. Fluid Mech., 74, 113135 (1998). George, K., John, T., “Steady extrusion of viscoelastic materials from an annular die”, J. Non-Newt. Fluid Mech., 154, 136152 (2008). Mitsoulis, E., Heng, F.L., “Extrudate swell of Newtonian fluids from converging and diverging annular dies”, Rheol. Acta, 26, 414417 (1987). Garcia-Rejon, A., DiRaddo, R.W., Ryan, M.E., “Effect of die geometry and flow characteristics on viscoelastic annular swell”, J. Non-Newt. Fluid Mech., 60, 107128 (1995). Mitsoulis, E., “Annular extrudate swell of pseudoplastic and viscoplastic fluids”, J. Non-Newt. Fluid Mech., 141, 138147 (2007). Ganvir, V., Lele, A., Thaokar, R., Gautham, B.P., “Prediction of extrudate swell in polymer melt extrusion using an Arbitrary Lagrangian Eulerian (ALE) based finite element method”, J. Non-Newt. Fluid Mech., 156, 2128 (2009). Intawong, N-T., Sombatsompop, N., “Experimental studies on radial extrudate swell and velocity profiles of flowing PS melt in an electro-magnetized die of an extrusion rheometer”, Polym. Eng. Sci., 44 (12), 22982307 (2004). Guo, J.L., Zhou, G.F., Zhou, Y.F., Yan, L., “The full three dimensional isothermal viscoelastic numerical simulation on the extrusion die swell of polymer profile”, Journal of Nanchang University (Engineering & Technology), 29 (1), 48 (2007). (in Chinese) Siddiqui, A.M., Mahmood, R., Ghori, Q.K., “Some exact solutions for the thin film flow of a PTT fluid”, Phys. Lett. A, 356, 353356 (2006). Takehiro, Y., Masakazu, I., Masaki, N., Kiyoji, N., Noriyasu, M., “Three-dimensional viscoelastic flows through a rectangular channel with a cavity”, J. Non-Newt. Fluid Mech., 114, 1331 (2003). Paulo, G.S.de, Tome, M.F., McKee, S., “A marker-and-cell approach to viscoelastic free surface flows using the PTT model”, J. Non-Newt. Fluid Mech., 147, 149174 (2007). Mitsoulis, E., “Annular extrudate swell of pseudoplastic and viscoplastic fluids”, Journal of Non-Newtonian Fluid Mechanics, 141, 138147 (2007). Ki, B.S., Seung, J.P., Seong, J.L., Kyung, H.A., Seung, J.L., “Numerical simulation of three-dimensional viscoelastic flow using the open boundary condition method in coextrusion process”, Journal of Non-Newtonian Fluid Mechanics, 99, 125144 (2001). Beverly, C.R., Tanner, R.I., “Numerical analysis of three-dimensional Newtonian extrudate swell”, Rheol. Acta, 30, 341356 (1991). Gifford, W.A., “A three-dimensional analysis of coextrusion”, Polym. Eng. Sci., 37 (2), 315320 (1997). Qin, S.X., Zhao, G.Q., Mu, Y., Xu, X.M., “Numerical simulation of driven and pressure flow in compound shaped part of co-extrusion process of polymer with metal insert”, Key Eng. Mater., 392-394, 299303 (2009). Sun, J., Smith, M.D., Armstrong, R.C., Brown, R.A., “Finite element method for viscoelastic flows based on the discrete adaptive viscoelastic stress splitting and the discontinuous Galerkin method: DAVSS-G/DG”, J. Non-Newt. Fluid Mech., 86, 281307 (1999). Sun, J.S., Phan-Thien, N., Tanner, R.I., “An adaptive viscoelastic stress splitting scheme and its applications: AVSS/SI and AVSS/SUPG”, J. Non-Newt. Fluid Mech., 65, 7591 (1996). Huang, S.H., Jiang, T.Q., Lu, C.J., Huang, J., “Experiments and numerical simulation of die swell for LDPE melt”, Journal of Shanghai Jiaotong University, 37 (4), 535640 (2003). (in Chinese)